Journal of Modern Power Systems and Clean Energy

ISSN 2196-5625 CN 32-1884/TK

网刊加载中。。。

使用Chrome浏览器效果最佳,继续浏览,你可能不会看到最佳的展示效果,

确定继续浏览么?

复制成功,请在其他浏览器进行阅读

Detection of Nonlinear Behavior Induced by Hard Limiting in Voltage Source Converters in Wind Farms Based on Higher-order Spectral Analysis  PDF

  • Zetian Zheng 1
  • Shaowei Huang 1
  • Qiangsheng Bu 2
  • Chen Shen 1
  • Jun Yan 1
1. the State Key Laboratory of Power Systems, Department of Electrical Engineering, Tsinghua University, Beijing 100084, China; 2. State Grid Jiangsu Electric Power Co., Ltd. Research Institute, Nanjing, China

Updated:2024-03-26

DOI:10.35833/MPCE.2022.000784

  • Full Text
  • Figs & Tabs
  • References
  • Authors
  • About
CITE
OUTLINE

Abstract

In recent years, sub-synchronous oscillation accidents caused by wind power integration have received extensive attention. The recorded constant-amplitude waveforms can be induced by either linear or nonlinear oscillation mechanisms. Hence, the nonlinear behavior needs to be distinguished prior to choosing the analysis method. Since the 1960s, the higher-order statistics (HOS) theory has become a powerful tool for the detection of nonlinear behavior (DNB) in production quality control wherein it has mainly been applied to mechanical condition monitoring and fault diagnosis. This study focuses on the hard limiters of the voltage source converter (VSC) control systems in the wind farms and attempts to detect the nonlinear behavior caused by bi- or uni-lateral saturation hard limiting using the HOS analysis. First, the conventional describing function is extended to obtain the detailed frequency domain information on the bi- and uni-lateral saturation hard limiting. Furthermore, the bi- and tri-spectra are introduced as the HOS, which are extended into bi- and tri-coherence spectra to eliminate the effects of the linear parts on the harmonic characteristics of hard limiting in the VSC control system, respectively. The effectiveness of the HOS in the DNB and the classification of the hard-limiting types is proven, and its detailed derivation and estimation procedure is presented. Finally, the quadratic and cubic phase coupling in the signals is illustrated, and the performance of the proposed method is evaluated and discussed.

I. Introduction

WITH the implementation of the energy development strategy in China, the installation capacity of wind power has increased [

1]. However, the stability problems caused by wind power integration are more serious, with sub-synchronous oscillation (SSO) being the most severe [2]. After SSO accidents, the amplitudes of the recorded accident waveforms tend to be constant because the transition process is very fast. Many studies have analyzed SSO based on the weak-damping assumption. Here, the researchers linearized the power-system dynamic equations and conducted studies using linear analysis methods, including the eigenvalue [3], dynamic equivalent [4], [5], impedance [6], [7], and complex torque coefficient [8], [9] methods. However, nonlinear effects such as hard limiting in the control system of the wind generators, while the induced self-sustained oscillation also exhibits a constant amplitude.

Current research on nonlinear oscillation analysis can be broadly classified into two categories.

1) Under small perturbations, constant-amplitude oscillations induced by the negative damping of the voltage source converter (VSC) generally begin with divergent oscillations near the equilibrium point of the system. Owing to the influence of the nonlinear factors in the VSC, such as hard limiting and pulse-width modulation (PWM) saturation, the oscillations will not continue to diverge, thereby resulting in constant-amplitude self-sustained oscillations [

10]. To address the constant-amplitude oscillations caused by hard limiting, the describing function (DF) can be used to model the nonlinear parts of the system, combined with a small-signal model of the system to analyze the characteristics of the oscillations. For example, [11] established an impedance model of a VSC considering the PWM saturation and analyzed the mechanism of the constant-amplitude oscillations. References [12] and [13] considered the bi-lateral hard limiting of the current and DC-voltage loops of a grid-connected VSC, respectively, and used the Nyquist criterion to obtain the approximate calculation formula for the oscillation frequency and amplitude.

2) Under large disturbances such as faults, self-sustained oscillations can be caused by the effects of nonlinear parts such as hard limiting. Meanwhile, if the corresponding hard-limiting parts are removed, the system can return to its original equilibrium point after the large disturbance. Current research on the mechanism of oscillations under large disturbances is still in the preliminary stage. Reference [

14] studied the self-sustained oscillations caused by the hard limiting of static var generators (SVGs) in double-fed induction generators after a fault and analyzed the influence of the system parameters on the non-smooth bifurcations. Reference [15] analyzed the instable oscillations of the grid-connected VSCs at the positive-damping equilibrium points after they experienced large disturbances. When the d-axis current reaches its hard limit, the outer control loop of the DC voltage loses its effectiveness, thereby resulting in electromagnetic oscillations in the system, accompanied by a continuous increase in the DC voltage.

The mechanism and corresponding analysis method of the nonlinear oscillations are different from those of linear oscillations. The direct analysis of the amplitude and frequency of the SSO using a linearized method and adoption of the corresponding measures to suppress the oscillation are unacceptable. Hence, before choosing the analysis method, it is important to determine the type of oscillation (linear or nonlinear) from the waveform records.

Since its emergence in the early 1960s, the higher-order statistics (HOS) theory has become a powerful tool for the condition monitoring and fault diagnosis of the mechanical equipment. Its applications include harmonic retrieval [

16], system identification [17], and feature extraction [18]. Furthermore, it is used for the detection of nonlinear behavior (DNB) [19]. Based on the bi-spectrum, the definitions of the bi-coherence spectrum and inverted bi-spectrum are presented in [20] and [21], respectively, and different statistical indicators based on HOS for DNB have been introduced.

Inspired by the aforementioned ideas in 1995, the researchers of power systems applied the bi-spectral analysis to the fault diagnosis and condition monitoring of the three-phase induction motors to analyze and identify motor asymmetric faults and stator winding failure [

22]. Subsequently, studies on the fault diagnosis using HOS have achieved considerable results. Reference [23] detected and identified asymmetric faults in induction motors by measuring the vibration data and analyzing the motor nonlinearity using the bi-coherence spectrum. In [24], considering the Gaussian and non-Gaussian noises in the mechanical signals, a new rolling bearing detection method was proposed, which integrated the bi-spectral analysis with the improved ensemble empirical mode decomposition.

In wind farms, [

25] used a modulated signal bi-spectrum detector to diagnose bearing faults in doubly-fed induction generators. Reference [26] used a bi-spectral analysis to identify single-point defects in rolling bearings. Reference [27] proposed an improved signal separation method based on the Vold-Kalman filter and HOS analysis for rotating mechanical systems under strong background noise.

However, the current applications of the bi-spectral analysis focus on the mechanical defects in the power systems, and few studies have analyzed the nonlinearity of the control modules in these systems, particularly in wind farms.

In this study, we focus on the hard limiters in the VSC, which is a crucial component of wind turbines and SVGs. To apply the HOS to the hard-limiting DNB in the control system of VSCs in wind farms, four aspects must be considered: ① suitability of the HOS for characterizing the nonlinearity of the hard limiters in the VSC control system and selection of an appropriate HOS; ② approach to determine the characteristics using only the waveforms collected at the terminal of the VSC; ③ method to detect the nonlinear behavior caused by the bi- or uni-lateral saturation hard limiting, which is considered as the source of the nonlinearity in this study; and ④ approach to improve the effectiveness and quality of the spectrum.

By fully addressing the aforementioned aspects, this study seeks to analytically prove the effectiveness of the HOS applied to the DNB induced by hard limiting in the VSC control system.

The rest of this paper is organized as follows. Section II extends the conventional DF and analyzes the two types of hard limiting. In Section III, the HOS is introduced. In Section IV, the VSC control system is modeled while the HOS-based hard-limiting DNB is analyzed and confirmed. The detailed procedures for the hard-limiting DNB in the VSC control system are described in Section V. Its effectiveness is demonstrated through case studies in Section VI. Finally, Section VII draws the conclusions.

II. Analysis of Harmonic Characteristics of Hard Limiting

The DF [

28] has been effectively used to analyze the characteristics of the self-sustained oscillations or limit cycles caused by nonlinearity. In the conventional DF, only the first-order Fourier series of the oscillation are maintained. In the DNB, however, the higher-order terms are required to distinguish the self-sustained oscillations induced by hard limiting from the negatively damped oscillations induced by incorrectly configured control parameters. Hence, the DF method is introduced and extended, and the harmonic characteristics of the bi- and uni-lateral saturation hard limiting are discussed.

A. Extended DF

As shown in Fig. 1, the input signal of the studied nonlinear part xt is assumed to be sinusoidal as expressed below:

xt=Asinωt (1)

Fig. 1  Typical hard limiter with a sinusoidal input.

where A and ω are the amplitude and frequency of the sinusoidal input signal, respectively.

The hard-limit output yt is a periodic nonsinusoidal signal, which can be expanded into a Fourier series as expressed below:

yt=A0+n=1Ancosnωt+Bnsinnωt (2)

where A0 is the magnitude of the DC component; and An and Bn are the coefficients of the cosine and sine parts of the nth Fourier harmonic, respectively.

In the conventional DF-based analysis method, the nonlinear part is considered oddly symmetrical and the linear part is considered low-pass. Consequently, A0=0, and yt can be approximated as:

yty1t=A1cosωt+B1sinωt=Y1ejφ1t (3)

where Y1=A12+B12 and φ1=arctanB1A1, and we have

A1=1π02πytcosωtdωtB1=1π02πytsinωtdωt (4)

The ratio of the first-order Fourier series of yt to the magnitude of the input signal is defined as the DF of the nonlinear part, i.e., NA.

NA=Y1Aejφ1 (5)

To extend the DF method, the higher-order harmonic coefficients An and Bn in (2) are calculated as:

An=1π02πxtcosnωtdωtBn=1π02πxtsinnωtdωt (6)

B. Bi-lateral Saturation Hard Limiting

As shown in Fig. 2(a), when a sine wave xt=Asinωt passes through a time-domain bi-lateral saturation hard limiter, its upper and lower parts that exceed the limits are set at the limit values, as expressed below:

yt=-a         xt<-axt    -axtaa           xt>a (7)

Fig. 2  Time-domain characteristics of bi- and uni-lateral saturation hard limiting. (a) Bi-lateral saturation. (b) Uni-lateral saturation.

where a a>0 is the upper hard limit.

Here, the bi-lateral saturation hard limiting is oddly symmetrical, and the output periodic signal is an odd function. Hence, the coefficients of the DC and cosine components in the Fourier series are 0, that is, An=0 (n=0,1,2,). According to (4), the fundamental Fourier coefficient of the output signal can be obtained as:

B1=1π02πxtsinωtdωt=2π0ϕAsin2ωtdωt+ϕπ-ϕasin ωtdωt+π-ϕπAsin2ωtdωt=2AπarcsinaA+aA1-aA2 (8)

The nth Fourier coefficient is calculated similarly, and the results show that

B2n=0n=1,2,3, (9)

The detailed calculation results are presented in Appendix A Table AI.

C. Uni-lateral Saturation Hard Limiting

As shown in Fig. 2(b), when a sine wave xt=Asinωt passes through a time-domain uni-lateral saturation hard limiter, the upper or lower (depending on the actual situation) part that exceeds the limit is set at the limit value. Without loss of generality, the corresponding equations of the upper hard limit are expressed as:

yt=xtxtA0+aA0+axt>A0+a (10)

where A0 is the offset of the sinusoidal input signal, and A0+a is the upper hard limit.

The nth Fourier coefficient is calculated in Section II-B, and the results show that:

A2n+1=0    n=0,1,2,B2n=0        n=0,1,2, (11)

The detailed calculation results are presented in Appendix A Table AII.

III. HOS Analysis

The definition of the HOS is introduced to analyze the harmonic characteristics of the nonlinear parts. The eigenfunction method is an important tool in statistical analysis since it can easily define the higher-order moments and cumulants.

The first joint eigenfunction of k continuous random variables x1,x2,,xk is defined as:

Φω1,ω2,,ωkEejω1x1+ω2x2++ωkxk=--fx1,x2,,xkejω1x1+ω2x2++ωkxkdx1dx2dxk (12)

where f· is the probability density function; and E· is the expectation function.

The kth-order moment and cumulant of k random variables are obtained as:

Ex1,x2,,xk=-jkkΦω1,ω2,,ωkω1ω2ωkω1=ω2==ωk=0 (13)
cumx1,x2,,xk=-jkklnΦω1,ω2,,ωkω1ω2ωkω1=ω2==ωk=0 (14)

For a stationary continuous random signal xt, we set x1=xt,x2=xt+τ1,,xk=xt+τk-1, which are later substituted into (14). Consequently, the kth-order cumulant of the random signal xt is expressed as:

ckxτ1,τ2,,τk-1=cumxt,xt+τ1,,xt+τk-1 (15)

The kth-order cumulant spectrum Skx is defined as the (k-1)-dimensional discrete Fourier transform of the kth-order cumulant, which is expressed as:

Skxω1,ω2,,ωk-1=τ1=-τk-1=-ckxτ1,τ2,,τk-1e-jω1τ1+ω2τ2++ωk-1τk-1 (16)

Generally, a higher-order cumulant spectrum is referred to as a higher-order spectrum. Particularly, the third-order spectrum S3xω1,ω2 is referred to as the bi-spectrum because it is a two-frequency energy spectrum, which is represented by Bω1,ω2. Similarly, the fourth-order spectrum S4xω1,ω2,ω3 is referred to as the tri-spectrum, expressed as Tω1,ω2,ω3.

IV. Hard-limiting DNB of VSC Control System

A. Modeling of VSC Control System

A typical direct-drive wind farm has 30-60 generators. Several direct-drive permanent magnet synchronous generators (PMSGs) are connected to form a string structure. The strings are further connected to a point of common coupling (PCC) and finally to the grid with a voltage level of 35 kV through a series of boosting transformers and transmission lines (equivalent to a set of impedances) [

29], as shown in Fig. 3.

Fig. 3  Model of wind farm connected to grid.

Figure 4 shows the structure of a typical direct-drive wind generator. Because the machine-side converter adopts the maximum power tracking control, its interaction with the grid is limited. As presented in [

30] and [31], the grid-related oscillation dynamics strongly depend on the DC capacitor and grid-side converter, and are slightly affected by the wind turbine, PMSG, and machine-side converter. Hence, the machine-side component, including the wind turbine, PMSG, and machine-side converter, is equivalent to a power source that outputs the power received by the wind turbine. Additionally, the grid-side converter is modeled as a VSC with its corresponding control system.

Fig. 4  Structure of a typical direct-drive wind generator.

Moreover, static reactive power compensation equipment, such as SVGs, is generally installed in wind farms. In this study, the SVG model adopts a double closed-loop control scheme. The d-axis control loop stabilizes the DC bus voltage, and the q-axis control loop varies based on the control mode. When the SVG operates in the constant-voltage control mode, the control target of the loop is the terminal voltage, and when it operates in the constant reactive power control mode, the control target of the loop is the output reactive power. Owing to its ability to produce high-quality output voltage waveforms with reduced harmonic distortion and low switching frequencies, many existing SVGs in wind farms are based on cascaded H-bridge topology [

32], [33]. The average-value model that ignores the PWM dynamics and models the cascaded H-bridge in the circuit as a controlled voltage source is adopted in this study. Therefore, the SVG is modeled as a VSC with its corresponding control system.

Thus, we obtain a unified VSC control system for the PMSGs and SVGs in the wind farms. The only difference is the choice of control targets in the d- and q-axis control loops, as shown in Fig. 5. In this study, four hard limiters are considered as the nonlinear parts: two in the inner control loop of the current and two in the outer control loop of the voltage, respectively.

Fig. 5  Structure of VSC control system.

In Fig. 5, vdc and vdcr* are the measured and reference DC bus voltages, respectively; C is the DC capacitor; vd (or Q) and vdr* (or Qr*) are the measured and reference d-axis terminal voltages (or the output reactive power), respectively; id and iq are the d- and q-axis components of the output current of the VSC, respectively; ed* and eq* are the d- and q-axis components of the output voltage reference of the VSC, respectively; v', v, vpcc, and vg are the voltages of the VSC output, terminal, PCC, and grid, respectively; θ is the output angle of the phase-locked loop (PLL); and R, L and R2, L2 are the equivalent resistances and inductances of the transformer and transmission line between the VSC output and terminal, respectively; Rl, Ll and Rl2, Ll2 are the equivalent resistances and inductances of the transmission line between the terminal and PCC, respectively; and Rg and Lg are the equivalent resistance and inductance of the transmission line between PCC and the grid, respectively.

B. Elimination of Effects from Linear Parts

While the nonlinear parts take effect inside the control system, the accident waveform record collected from the phasor measurement unit provides only voltage and current information at the terminal of the VSC. Hence, the relationship between the HOS of the hard-limit output and that of the terminal electrical quantities must be derived.

When the hard limiting in the d-axis inner control loop of the current takes effect, v¯d is produced by the nonlinear part after proportional-integral (PI) control, as shown in Fig. 5. Assuming that the PLL performs satisfactorily and the VSC is synchronized with the system, the d-axis frame model of the main circuit is expressed as:

Δvd'-Δvd=sL+RΔid-ω0LΔiq (17)

where Δ represents the small-signal perturbation; and ω0 is the fundamental angular frequency.

Without considering the dynamics of the PWM, we assume that the VSC output tracks the reference signal:

Δvd'=Δed* (18)
Δvq'=Δeq* (19)

From Fig. 5, the relationship between v¯d and vd' is obtained as:

Δvd'=Δed*=Δv¯d-ω0LΔiq+Δvd (20)

where v¯d is the hard-limit output in the d-axis inner control loop of the current.

Then, we obtain the relationship between id and v¯d by combining (17) and (20), as expressed below:

Δid=1sL+RΔv¯d (21)

Similarly, when the hard limiting in the q-axis inner control loop of the current takes effect, the relationship between iq and v¯q is expressed as:

Δiq=1sL+RΔv¯q (22)

where v¯q is the hard-limit output in the q-axis inner control loop of the current.

When the hard limit in the d-axis outer control loop of the voltage takes effect, the output i¯d is produced by the nonlinear part after PI control. Thus, (20) can be rewritten as:

Δvd'=GisΔi¯d-Δid-ω0LΔiq+Δvd (23)

where i¯d is the output of the hard limit in the d-axis outer control loop of the voltage; and Gi(s) is the transfer function of the inner control loop of the current.

Then, we obtain the relationship between id and i¯d by combining (17) and (23), as expressed below:

Δid=GissL+R+GisΔi¯d (24)

Similarly, when the hard limit in the q-axis outer control loop of the voltage takes effect, the relationship between iq and i¯q is expressed as:

Δiq=GissL+R+GisΔi¯q (25)

where i¯q is the output of the hard limit in the q-axis outer control loop of the voltage.

Hence, the output of the nonlinear part in the control system can always be obtained from the current measured at the terminal after passing through the linear part by combining (21), (22), (24), and (25). The system shown in Fig. 6 illustrates how to eliminate the effects of the linear parts on HOS of the terminal current, wherein en and yn are the input and output, respectively; and Hz represents the linear time-invariant part.

Fig. 6  Typical linear system.

The relationship between the HOS of en and that of yn can be obtained based on the definition and properties of the HOS [

34] as expressed below:

Skyω1,ω2,,ωk-1=Skeω1,ω2,,ωk-1Hω1Hωk-1H*ω1+ω2++ωk-1 (26)

where Skeω1,ω2,,ωk-1 and Skyω1,ω2,,ωk-1 are the kth-order cumulant spectra of en and yn, respectively; Hω is the continuous transfer function of Hz; and * is the conjugation operator. Let k=2,3 in (26). S2ω is denoted by Pω and S3ω1,ω2 is denoted by Bω. Consequently, we can obtain:

Pyω=PeωHω2 (27)
Byω1,ω2=Beω1,ω2Hω1Hω2H*ω1+ω2 (28)

Define bicω1,ω2 as the bi-coherence [

35] at the frequency pair ω1,ω2, which is calculated as:

bicω1,ω2=Bω1,ω2Pω1Pω2Pω1+ω2 (29)

By combining (27), (28), and (29), we can obtain:

bicyω1,ω2=biceω1,ω2 (30)

Further, we define the tri-coherence tricω1,ω2,ω3 [

35] as:

tricω1,ω2,ω3=Tω1,ω2,ω3Pω1+ω2+ω3Pω1Pω2Pω3 (31)

Similarly, it can be proven that tricyω1,ω2,ω3=triceω1,ω2,ω3.

Consequently, the linear part does not change the bi- or tri-coherence of the system. Obtaining the waveforms of id and iq at the terminal of the VSC and performing the HOS analysis enable the detection of the nonlinearity of the VSC system.

C. Bi-coherence Spectrum and Uni-lateral Saturation Hard-limiting Detection

Based on the results in Section II-C, when a self-sustained oscillation occurs, the output of the uni-lateral saturation hard limiter contains the second harmonic of the oscillation frequency with the same phase as that of the fundamental frequency. Without loss of generality, let its initial phase be zero, i.e.,

yt=B1sin2πft+B2sin2π×2ft (32)

where f denotes the oscillation frequency; and B1 and B2 are the coefficients of the fundamental and second Fourier harmonics, respectively.

A Fourier transform is performed twice on the second-order autocorrelation function of yt. Consequently, the bi-spectrum of yt can be expressed as:

Byω1,ω2=--f01/fytyt+τ1yt+τ2dte-j2πω1+ω2dτ1dτ2 (33)

The bi-spectrum has 12 symmetrical regions [

35], as shown in Fig. 7(a).

Fig. 7  Related figures of bi-spectrum and bi-coherence spectrum. (a) Symmetric region of bi-spectrum. (b) Peak of bi-coherence spectrum.

Therefore, to completely describe the entire bi-spectrum, only a symmetrical region in the resulting Byω1,ω2 is required for analysis. Considering ω1,ω2|0ω1ω2, Byω1,ω2 can be calculated as:

Byω1,ω2=π4jB12B2δω1-2πfδω2-2πf (34)

where δ is the Dirac delta function [

36], which is expressed as:

δx=0    x0-δxdx=1 (35)

where δx is a finite maximum at x=0 when the input signal is discretized. In (34), Byω1, ω2 is a finite maximum if and only if ω1=ω2=2πf; otherwise, it is zero. Therefore, a peak can be observed at the x-y coordinates 2πf, 2πf in the three-dimensional (3D) graph of ω1,ω2,Byω1,ω20ω1ω2. Furthermore, because 2πf,

2πf is on the symmetry axis ω1=ω2, when the area is extended to ω1,ω2ω10,ω20 (areas 1 and 2 in Fig. 7(a)), the peak at 2πf, 2πf still exists.

The power spectrum of yt is derived as:

Pyω=-f01/fytyt+τdte-j2πωdτ=12B12π2δω-2πf+12B22π2δω-4πf (36)

Based on (29), (34), and (36), the bi-coherence spectrum ω1,ω2ω10,ω20 can be calculated as:

bicyω1,ω2=Byω1,ω2Pyω1+ω2Pyω1Pyω2=       π4B12B2δω1-2πfδω2-2πf/     12B12π2δω1-2πf+12B22π2δω1-4πf/     12B12π2δω2-2πf+12B22π2δω2-4πf/   12B12π2δω1+ω2-2πf+12B22π2δω1+ω2-4πf                                                                                    ω10,ω20 (37)

Since F-1δω=12π, it can be proven that the bi-coherence spectrum of yt reaches its peak if and only if ω1=ω2=2πf, as shown in Fig. 7(b). It is further calculated as:

bicy2πf,2πf=π4B12B212π12π12B12π212π12B1212ππ212B2212ππ2=1 (38)

In (35), bicyω1,ω20. Therefore, the range of the corresponding bi-coherence value of each x-y coordinate in the bicoherence spectrum is 0,1. The larger the bi-coherence value, the stronger the nonlinear phase coupling between the two frequencies corresponding to the coordinate (the stronger the nonlinearity).

In (32), when yt extends to yt=n=1+Bnsin (2πnft), it can be proven that in the 3D graph of ω1,ω2, bicyω1,ω2ω10,ω20, there are peaks at the x-y coordinates i2πf,j2πf(i,j=1,2,3,), and their corresponding bi-coherence values are equal to one. Furthermore, it can be proven that the conclusion remains unchanged when yt=n=1+A2ncos2π×2nft+n=0+B2n+1sin2π2n+1ft, which accurately represents the output of the uni-lateral saturation hard limiter according to (11).

D. Tri-coherence Spectrum and Bi-lateral Saturation Hard-limiting Detection

Based on the results in Section II-B, when yt=n=0+B2n+1sin2π2n+1ft, which accurately represents the output of the bi-lateral saturation hard limiter according to (9), it can be proven that in the 4D graph of ω1,ω2,ω3,tricyω1,ω2,ω3ω10,ω20,ω30, peaks exist at the x-y-z coordinates 2i+1×2πf, 2j+1×2πf, 2k+1×2πfi,j,k=0,1,2,, with the peak values of 1. The range of the corresponding tri-coherence value of each x-y-z coordinate in the tri-coherence spectrum is 0,1. Furthermore, the larger the tri-coherence value, the stronger the nonlinear phase coupling among the three frequencies corresponding to the coordinate (the stronger the nonlinearity).

The detailed derivation process is presented in Appendix B.

E. Nonlinearity Detection and Classification of VSC Control Systems

Based on Section IV-B, the nonlinearity induced by hard limiting in the VSC control system can be detected by transforming the current waveform into id and iq at the terminal of the VSC and performing the HOS (bi-/tri-coherence) analysis. The nonlinearity of id represents the nonlinearity in the d-axis control loop of the VSC control system, whereas the nonlinearity of iq represents the nonlinearity in the q-axis.

Table I summarizes the results obtained in Sections IV-C and IV-D. By examining the relationship between the peak coordinates of the HOS and harmonic characteristics of the analyzed signal, the bi-coherence spectrum is used to detect the “phase coupling”. Thus, the harmonics with frequencies f1+f2=f3 and phases φ1+φ2=φ3 are simultaneously satisfied. Furthermore, the tri-coherence spectrum is used to detect whether f1+f2+f3=f4 and φ1+φ2+φ3=φ4 are both satisfied. It is noteworthy that the phase equation is sufficient but unnecessary in the phase coupling phenomenon, which is thoroughly explained in Section VI-A.

TABLE I  Characteristics of Bi-coherence and Tri-coherence Spectra
TypePhase couplingUni-lateral saturationBi-lateral saturation
Fourier series yt=n=1+A2ncos2π×2nft+n=0+B2n+1sin2π2n+1ft yt=n=0+B2n+1sin2π2n+1ft
Bi-coherence f1+f2=f3,φ1+φ2=φ3 Peaks No peaks
Tri-coherence f1+f2+f3=f4,φ1+φ2+φ3=φ4 Peaks Peaks

Table I shows that the output of the uni-lateral saturation hard limiter contains harmonics of all orders and satisfies both the quadratic and cubic phase coupling. Thus, peaks exist both in the bi- and tri-coherence spectra. Hence, the bi-coherence spectrum must be first examined, followed by the tri-coherence spectrum. Nonlinearity can be assessed and classified, as shown in Fig. 8.

Fig. 8  Flowchart of nonlinearity detection and classification.

In this study, when detection is performed at the VSC terminal, only hard limiting related nonlinearity exists in the control system (average-value model), and all the other function blocks in the control system can be represented by the transfer functions (linear blocks). Therefore, we assume that hard limiting is the only source of nonlinearity in Fig. 8.

V. Procedure for Hard-limiting DNB of VSC Control System

The specific procedure for applying the hard-limiting DNB to the VSC control system is based on the theoretical analysis described in Section IV. The corresponding steps are summarized in Fig. 9, and the detailed process is as follows.

Fig. 9  Key points of procedure for hard-limiting DNB.

To obtain the time series for the calculation, we first collect the accident waveform records of the current at the terminal of the VSC. The studied signals xak, xbk, and xck are sampled from the three-phase currents iat, ibt, and ict, respectively. We begin the sampling from the initial time t0 and set the sampling interval to Δt. The sampling point is denoted by k, and the sample length is L, i.e.,

xak=iat0+kΔtxbk=ibt0+kΔtxck=ict0+kΔtk=1,2,,L (39)

Step 1:   perform a dq transformation on the three-phase sampling signals xak, xbk, and xck. The initial phase θ0 can be obtained by applying a PLL algorithm as expressed below:

xdk     xqk    x0kT=23cosθk+θ0cosθk-23πcosθk+23π-sinθk+θ0-sinθk-23π-sinθk+23π121212xak    xbk    xckT (40)

where θk=ω0kΔt+θ0; and xdk, xqk, and x0k are the transformation results in the dq0 coordinate system.

Step 2:   set the transformation result xdk of the previous step as the subsequent signal processing object, i.e.,

xk=xdk (41)

Step 3:   divide xk into M segments, wherein each segment length is N (MN=L). Each segment is recorded as xil (i=1,2,,M;l=1,2,,N).

Step 4:   select an appropriate window function, such as a Hanning window, which is expressed as:

wl=121+cos2πlN-1 (42)

Multiply each segment of the signal by the window function, and use the obtained results x'il for further calculations to reduce leakage errors, as expressed below:

x'il=xilwl (43)

Step 5:   for each segment x'il, subtract its mean value as:

x^il=x'il-x¯il (44)

Step 6:   perform the fast Fourier transform (FFT) on each segment x^il:

Xki=1Nl=1Nx^ile-i2πkl/Nk=1,2,,N/2,i=1,2,,M (45)

Step 7:   when dealing with the FFT results, consider a small parameter σ (for instance, σ=0.001). Sweep i=1,2,,M, and for any k, if Xki<σmaxk=1,2,,N/2Xki, then let Xki=σ2maxk=1,2,,N/2Xki. This step can further increase the difference in the order of magnitudes between the white noise and peak value in the spectrum. Thus, the analysis and conclusion of the peak value is not affected by the appearance of the values close to 0/0 in the bi- and tri-coherence spectrum.

Step 8:   the estimated values of the power spectrum, bi-spectrum, and tri-spectrum of xk are expressed as:

P^m=1Mi=1MXmiXm*i (46)
B^m,n=1Mi=1MXmiXniXm+n*i (47)
T^m,n,o=1Mi=1MXmiXniXoiXm+n+o*i (48)

Step 9:   calculate the bi-coherence spectrum as:

b^icm,n=B^m,nP^m+nP^mP^n (49)

Step 10:   the obtained bi-coherence spectrum is a 3D graph. Its x-y coordinates are the m,n frequencies and its z coordinate is the corresponding bi-coherence value, whose theoretical value range is 0,1.

Step 11:   calculate the tri-coherence spectrum as:

t^ricm,n,o=T^m,n,oP^m+n+oP^mP^nP^o (50)

Step 12:   determination of DNB thresholds.

Peak value threshold: define σp as the nonlinear threshold (preferably 0.3). A peak in the bi- or tri-coherence spectrum whose value is greater than σp is considered to characterize the existence of the quadratic or cubic phase coupling, and the coordinates of the peak represent the corresponding frequencies. The conclusion and classification of the nonlinearity can be completed using the process provided in Fig. 8.

Flatness threshold: in addition to the index σp for checking the peak values, a nonlinear index μ can be defined based on the graph flatness [

37] for the bi-coherence spectrum:

μb^icmax2-b^ic¯2+2σb^ic2 (51)

where b^ic¯2 is the average of the estimated squared bi-coherence; and σb^ic2 is the standard deviation of b^ic2. Based on the power-quality standard [

38], the corresponding result of the index μ is 0.0180 for 35 kV networks obtained by simulation, while considering the measurement of the output currents with the background harmonics. Consequently, we set the threshold μp=0.10 to keep a sufficient margin. Here, the nonlinear behavior is considered existent if μ>μp.

For a single VSC with slight harmonic pollution, the peak value threshold is enough for the DNB. Contrarily, for the multi-VSCs or complex systems, the flatness threshold can be considered. However, only the existence of nonlinear behavior in the control system can be verified, and misjudgment may occur when using the flatness threshold to locate the nonlinear oscillatory source, which is illustrated and discussed in Section VI-C.

Step 13:   the above steps implement the DNB on the d-axis control loop of the VSC control system. To study the q-axis control loop, return to Step 2: let xk=xqk and repeat Steps 3-12.

Consequently, the nonlinear behavior caused by the bi- or uni-lateral saturation hard limiting in the d- or q-axis control loop can be detected. Some of these steps increase the resolution and effectiveness of the HOS. Step 4 reduces the spectrum leakage, Steps 6 and 8 eliminate the effect of the random phases, and Step 7 adds credibility to the presence of peaks.

VI. Case Studies

Three cases are presented and discussed to evaluate the effectiveness of the proposed method for the hard-limiting DNB in the VSC control system.

Case 1 considers an artificially constructed signal that is abstracted from the harmonic characteristics of the uni-lateral saturation hard limiter in Section II-C to further demonstrate whether φ1+φ2=φ3 is a necessary condition of phase coupling in HOS.

Case 2 sets up a grid-connected PMSG model to demonstrate the effectiveness of the proposed process in detecting the nonlinearity owing to the uni-lateral saturation hard limiting by collecting accident waveform records at the terminals of the VSCs. In this case, the nonlinearity owing to bi-lateral saturation hard limiting is also detected later using the tri-coherence spectrum.

Case 3 sets up an IEEE 9-bus system with three SVGs and one static VAR compensator (SVC), wherein the self-sustained oscillation is induced by two SVGs. This case shows that the HOS can only detect the presence of nonlinearity but not locate the source of the nonlinear oscillations.

A. Case 1

Here, the following signal is considered in MATLAB.

xn=cos2πf1n/fs+φ1+w1n+cos2πf2n/fs+φ2+w2n+cos2πf3n/fs+φ3+w3n (52)

where f1=0.6381 Hz;  f2=0.8345 Hz;  f3=f1+f2; fs is the sampling frequency; and win i=1,2,3 is a -20 dB Gaussian white noise.

In (52), let φ1=φ2=φ3=0, and record xn as x1n. Further, let φ1=φ2=0, φ3=π/2, and record xn as x2n.

Figure 10(a) shows the frequency spectrum of x1n, from which three frequency components f1, f2, and f3 can be found. However, the relationship between them cannot be determined. Figure 10(b) shows the bi-coherence spectrum of x1n as a contour map. Its peak appears at (f1, f2) and (f2, f1), with a peak value of 1.0. Hence, the power at f3 comes entirely from the quadratic phase coupling of f1 and f2, which confirms the conclusion in Section IV-C. Figure 10(c) and (d) shows the frequency and bi-coherence spectrum of x2n, which are almost the same as those in Fig. 10(a) and (b), respectively. However, in our initial setting, φ1+φ2=φ3 in x1n, whereas φ1+φ2φ3 in x2n.

Fig. 10  Detection of quadratic phase coupling. (a) Frequency spectrum of x1n. (b) Bi-coherence spectrum of x1n. (c) Frequency spectrum of x2n. (d) Bi-coherence spectrum of x2n.

Previous studies [

39]-[41] considered a quadratic phase coupling which is equivalent to f1+f2=f3 and φ1+φ2=φ3. In this case, when the condition extends to “φ1+φ2-φ3 is constant”, the phase coupling phenomenon and the peaks will also exist in the bi-coherence spectrum. Therefore, the phase equation is a sufficient but not unnecessary condition.

In this study, according to (11), the output of the uni-lateral saturation hard limiter can be expressed as:

yt=n=1+A2ncos2π×2nft+n=0+B2n+1sin2π2n+1ft=n=1+A2nsin2π×2nft+π2+n=0+B2n+1sin2π2n+1ft (53)

Therefore, the peaks with values greater than σp will appear in its bi-coherence spectrum similarly to a “chessboard” because any two integer multiples of the fundamental frequency have the property of the quadratic phase coupling. However, if the bi-coherence could detect only those satisfying f1+f2=f3 and φ1+φ2=φ3 simultaneously, no peaks in the bi-coherence spectrum will exist, which does not match the actual situation.

B. Case 2

A detailed grid-connected PMSG model is set up in PSCAD/EMTDC with the structure of the VSC control system shown in Fig. 11.

Fig. 11  PMSG model for DNB.

First, the parameters are adjusted to make the hard limit of PI in the d-axis outer control loop of the voltage take effect. Table II presents the parameter settings.

TABLE II  Parameters in Control Blocks of PMSG
SymbolDescriptionValue
Vg Grid voltage 0.69 kV
f0 Fundamental frequency 50 Hz
PN Rated capacity of PMSG 1.5 MW
HPLLs PLL 500+900/s
P Active power 0.34 MW
C DC capacitor 200 mF
R Connection resistance 0.001 Ω
L Connection inductance 0.35 mH
Rg Grid-side resistance 0.005 Ω
Lg Grid-side inductance 0.4 mH
Gdcs DC-voltage controller 9+500/s
Gqs Reactive power controller 0.3+50.28/s
Gis Inner-loop current controller 0.012+12.5/s

The simulation is implemented as follows.

t=0  s: use a voltage source to charge the DC capacitor. In the initial state, the PMSG is off-grid, and the active power and reactive power references are both zero.

t= 0.2 s: the DC capacitor side is switched to the power source, and the PMSG is connected to the grid.

t= 1.0 s: the active power is set to be 0.34 MW.

t= 2.0 s: Gis is set to be 0.2+20/s.

t= 4.0 s: Gis is set to be 0.012+12.5/s.

Figure 12 shows the d-axis current reference idref. According to Fig. 12, the oscillation can be divided into two stages.

Fig. 12  d-axis current reference.

Stage I: after t=4.0  s and before any hard limit takes effect, the system can be analyzed with small-signal linearization. Based on the parameters presented in Table II, the dominating mode of the system state space model is 6.32±j190.06, which is on the right side of the imaginary axis. It induces a divergent oscillation and activates the hard limits of PI.

Stage II: after the oscillation causes idref to reach the hard limit of PI in the d-axis outer control loop of the voltage, the hard limit becomes uni-laterally saturated and the system cannot be analyzed by the linearization method. At the PCC, a constant-amplitude self-sustained oscillation of 33.8 Hz is observed, as shown in Fig. 13. The oscillating frequency can be computed and examined using the DF-based Nyquist criterion [

13]. Since the transient process is short, the collected accident waveform record may only contain the constant-amplitude part, which superficially resembles a weakly-damped linear oscillation.

Fig. 13  Voltage and current at PCC in d-axis.

Figure 14 shows the bi-coherence and power spectra of idref and idpcc. The bi-coherence spectra are shown as contour maps whose x-y axis range is restricted to f1,f2|ω10,ω20, respectively. As shown in Fig. 14(a), the bi-coherence spectrum of idref resembles a “chessboard”. Hence, the quadratic phase coupling exists between any two integer multiples of the fundamental frequency. Generally, the transfer function of the VSC control system is low-pass, which can be observed by comparing Fig. 14(b) and (d). However, the bi-coherence spectrum of idpcc can still maintain the property of the quadratic phase coupling, as shown in Fig. 14(c). Here, the peaks of the higher harmonics disappear because their amplitudes are too small to maintain a distinction from the background noise. Hence, Fig. 14(c) confirms that nonlinearity owing to uni-lateral saturation hard limit can be detected by collecting accident waveforms at the terminal of VSCs based on the HOS analysis. In this case, nonlinearity exists in the d-axis control loop of the VSC system.

Fig. 14  Bi-coherence and power spectra of idref and idpcc. (a) Bi-coherence spectrum of idref. (b) Power spectrum of idref. (c) Bi-coherence spectrum of idpcc. (d) Power spectrum of idpcc.

Next, we adjust the parameters to examine the effectiveness of the tri-coherence spectrum. We set the upper and lower limits of the d-axis current reference idref to ± and set those of the d-axis voltage reference vdref to ±0.08. Then, the hard limit of vdref takes effect and induces a bi-lateral saturation self-sustained oscillation.

Figure 15(a) shows the bi-coherence spectrum of idpcc, with a global maximum value of only 0.002134, which is much smaller than the nonlinear threshold σp. Consequently, we consider that no peaks representing the nonlinearity exist in the contour map and no uni-lateral saturation exists.

Fig. 15  Bi- and tri-coherence spectra of idpcc. (a) Bi-coherence spectrum of idpcc. (b) Tri-coherence spectrum of idpcc.

Figure 15(b) shows the tri-coherence spectrum of idpcc, with a peak at the x-y-z coordinates (33.8 Hz, 33.8 Hz, 33.8 Hz). Here, the aliasing error expands the peak range, thereby making the maximum value larger than one; however, it still demonstrates the bi-lateral saturation in the system, which proves the effectiveness of the calculation and classification processes proposed in Section IV-E and Section V.

C. Case 3

Figure 16 illustrates a topology derived from the IEEE 9-bus system, wherein the network parameters of buses 1-9 correspond with the IEEE 9-bus system. Contrary to the IEEE benchmark system, an SVC is located at the medium-voltage bus 13 for reactive power adjustment, and two grid-connected wind farms and their corresponding SVGs are located at buses 14 and 15. The PMSGs and SVGs use the double-loop VSC control, as presented in Section IV-A. The control parameters are presented in Table III.

Fig. 16  Topology of case study system.

TABLE III  Parameters of Network and VSC Control for Case study System
ParameterValue
KPυd, KIυd, KPυq, KIυq (SVG voltage control loop) 2.5 p.u., 1000 p.u., 2 p.u., 20 p.u.
Vref1, Vref2, Vref3 (reference of terminal voltage control)

1.005 p.u., 1.005 p.u., 1.005 p.u. (t<2 s)

1.005 p.u., 1.000 p.u., 1.005 p.u. (t2 s)

KPi, KIi (current control loop) 40 p.u., 6250 p.u.
Xl1, Xl2, Xl3 (connection impedance) 0.0051 p.u., 0.0038 p.u., 0.0256 p.u.
R6-10, R8-11 (line resistance) 0.0017 p.u., 0.0054 p.u.
X6-10, X8-11 (line impedance) 0.0092 p.u., 0.0178 p.u.
XT1, XT2, XT3 (transformer impedance) 0.0586 p.u., 0.0586 p.u., 0.0576 p.u.

The case study system is simulated using PSCAD/EMTDC. In the case, a self-sustained oscillation is induced by mismatching the reference terminal voltages of SVG1 and SVG2, which are denoted by Vref1 and Vref2, respectively. Figure 17 shows the oscillation-related waveforms of the voltages and currents. In Fig. 17(a), when t<2 s, set Vref1=Vref2=1.005 p.u., thereby stabilizing the voltage amplitudes V14, Vpcc1, and Vpcc2. Contrarily, when t2.45 s, set Vref1=1.005 p.u. and Vref2=1.000 p.u.. Consequently, V14 deviates from the previous equilibrium point, which fluctuates in the range of 0.97-1.04 p.u.. Based on Fig. 17(b), i14 significantly increases when t2.45  s. The major harmonic frequencies of the instantaneous current i14 are 17.5  Hz and 82.5  Hz 50  Hz±32.5  Hz, and the sub-synchronous current flows from the VSCs to the networks. As shown in Fig. 17(c), the SVC stimulates a linear LC resonance, whose combination of capacitance and network inductance matches the sub-synchronous frequency of the VSCs. Additionally, the current waveform of the SVC shows that the SVC functions as a harmonic amplifier.

Fig. 17  Voltage and current waveforms. (a) Voltage amplitude. (b) Three-phase. (c) Three-phase currents of SVC.

Calculated as (51), Table IV shows that the nonlinear index μ>μp for the SVGs and μμp for the SVC, which indicates that the SVGs are the main contributors to the nonlinearity. Hence, the nonlinearity is detected in the system. However, the hard limiting takes effect in only two SVGs. Therefore, the proposed method can only detect the nonlinear behavior in the system, but not its location. Furthermore, the precise localization method needs to be developed in future studies.

TABLE IV  Comparison of μ for Harmonic Source Searching
Partμ (p.u.)Threshold μp (p.u.)
SVG1 0.3727 0.1000
SVG2 0.1598
SVG3 0.1043
SVC 0.0079

VII. Conclusion

This study proposes a method based on the HOS analysis for the hard-limiting DNB in the VSC control system in wind farms, wherein the PMSGs and VSGs are modeled using a unified VSC control model. The contributions are summarized as follows.

1) The effectiveness of the bi- and tri-spectra is proven when characterizing the nonlinear behavior induced by hard limiting in the VSC control system.

2) The effects of the linear parts in the VSC control system are eliminated using the bi- and tri-coherence, which facilitate the hard-limiting DNB with only the waveforms collecte at the terminal of the VSC. Additionally, the phase equation is proven to be a sufficient and unnecessary condition of the phase coupling phenomenon, which is not examined thoroughly in previous studies.

3) The detailed procedure is proposed to detect and classify the nonlinear behavior caused by the bi- or uni-lateral saturation hard limiting.

4) The data processing problems are solved in the DNB procedures based on the HOS, such as the “0/0” phenomenon and spectrum leakage, to improve the resolution and quality of the spectra.

Future studies include two aspects.

First, the HOS analysis is a measurement-based method. However, the extension of the proposed method to any equipment requires further examination. This study confirms the applicability of the HOS-based DNB to VSCs, which exist in energy storage equipment and flexible DC transmission systems, in addition to wind power systems as discussed in this study. Here, we do not investigate eliminating the effects of different parts of other devices (in the path from the hard limiter to the terminal).

Second, as discussed in Section VI-C, the DNB is not sufficient to constitute an effective control measure and generator tripping strategy when a self-sustained oscillation accident occurs. Therefore, a nonlinear oscillatory source localization method needs to be introduced.

Appendix

Appendix A

TABLE AI  The nth Component of Bi-lateral Saturation Hard Limit Output
nAnBn
0 0 0
1 0 2a1-a2A2+2AarcsinaAπ
2 0 0
3 0 4a1-a2A23/23π
4 0 0
5 0 4a1-a2A2(8a4-11a2A2+3A4)15A4π
6 0 0
7 0 48acos7arcsinaA+28Asin6arcsinaA-21Asin8arcsinaA84π
TABLE AII  The nth Component of Uni-lateral Saturation Hard Limit Output
nAnBn
0 a+2A0-2π1-a2A2A-2aπarcsinaA a+2A0-2π1-a2A2A-2aπarcsinaA
1 0 2a1-a2A2+Aπ+2AarcsinaA2π
2 21-a2A2-a2+A23Aπ 0
3 0 2a1-a2A23/23π
4 21-a2A26a4-7a2A2+A415A3π 0
5 0 2a3+8a4A4-11a2A21-a2A215π
6 21-a2A2105A5π-80a6+128a4A2-51a2A4+3A6 0
7 0 1168π48acos7arcsinaA+28Asin6arcsinaA-21Asin8arcsinaA

Appendix B

Based on the results in Section II-B, when self-sustained oscillation occurs, the output yt of the bi-lateral saturation hard limiter contains the 3rd and 5th harmonics of the oscillation frequency with the same phase as that of the fundamental frequency. Without loss of generality, let its initial phase be zero, that is,

yt=B1sin2πft+B3sin2π×3ft+B5sin2π×5ft (B1)

where B1, B3, and B5 are the Fourier coefficients of the 1st, 3rd, and 5th Fourier harmonics, respectively.

The Fourier transform is performed three times on the third-order autocorrelation function of yt. Thus, the tri-spectrum of yt can be derived as:

Rfτ1,τ2,τ3=f01/fytyt+τ1yt+τ2yt+τ3dtTyω1,ω2,ω3=---Rfτ1,τ2,τ3e-j2πω1+ω2+ω3dτ1dτ2dτ3 (B2)

Additionally, the tri-spectrum has 96 symmetrical regions [

35]. Therefore, a symmetrical region in the result of Tyω1,ω2,ω3 can be used for analysis to completely describe the entire trispectrum. Considering ω1,ω2,ω30ω1ω2ω3, Tyω1,ω2,ω3 is calculated as:

Tyω1,ω2,ω3=B13B3π3/242δω1-2πfδω2-2πfδω3-2πf+B12B3B5π3/242δω1-2πfδω2-2πfδω3-6πf (B3)

Tyω1,ω2,ω3 is a finite maximum if and only if ω1=ω2=ω3=2πf or ω1=ω2=2πf,ω3=6πf; otherwise, it is zero. Therefore, two peaks can be observed at the x-y-z coordinates 2πf,2πf,2πf and 2πf,2πf,6πf in the four-dimensional (4D) graph of ω1,ω2,ω3,Tyω1,ω2,ω30ω1ω2ω3. For an easy intuitive visualization of the graphs, when the area is extended to ω1,ω2,ω3ω10,ω20,ω30, four peaks at 2πf,2πf,2πf, 2πf,2πf,6πf, 2πf,6πf,2πf, and 6πf,2πf,2πf occur.

The power spectrum of yt is expressed as:

Pyω=-f01/fytyt+τdte-j2πωdτ=12B12π2δω-2πf+12B32π2δω-6πf+12B52π2δω-10πf (B4)

According to (B3) and (B4), the tri-coherence spectrum ω1,ω2,ω3ω10,ω20,ω30 can be expressed as:

tricyω1,ω2,ω3=Tyω1,ω2,ω3Pyω1+ω2+ω3Pyω1Pyω2Pyω3 (B5)

In (B5), the tri-coherence spectrum of yt reaches its peak if and only if ω1,ω2,ω3=2πf,2πf,2πf, 2πf,2πf,6πf, 2πf,6πf,2πf, or 6πf,2πf,2πf.

When ω1,ω2,ω3=2πf,2πf,2πf, (B5) is calculated as (B6). When ω1,ω2,ω3=2πf,2πf,6πf, 2πf,6πf,2πf, or 6πf,2πf,2πf, (B5) is calculated as (B7).

tricy2πf,2πf,2πf=B13B3π3/24212π12π12π12B32π212π12B1212ππ212B1212ππ212B1212ππ2=1 (B6)
tricy2πf,2πf,6πf=B12B3B5π3/24212π12π12π12B52π212π12B3212ππ212B1212ππ212B1212ππ2=1 (B7)

In (B5), tricyω1,ω2,ω30. Therefore, the range of the corresponding tri-coherence value of each x-y-z coordinate in the tri-coherence spectrum is 0,1. The larger the tri-coherence value, the stronger the nonlinear phase coupling among the three frequencies corresponding to the coordinate (the stronger the nonlinearity).

In (B1), when yt extends to yt=n=0+B2n+1sin2π2n+1ft, according to (9), which accurately represents the output of the bi-lateral saturation hard limiter, it can be proved that in the 4D graph of ω1,ω2,ω3,tricyω1,ω2,ω3ω10,ω20,ω30, peaks exist at the x-y-z coordinates 2i+12πf,2j+12πf,2k+12πf(i,j,k=0,1,2,).

References

1

S. Zhang, J. Wei, X. Chen et al., “China in global wind power development: role, status and impact,” Renewable and Sustainable Energy Reviews, vol. 127, p. 109881, Jul. 2020. [Baidu Scholar] 

2

J. Shair, X. Xie, L. Wang et al., “Overview of emerging subsynchronous oscillations in practical wind power systems,” Renewable and Sustainable Energy Reviews, vol. 99, pp. 159-168, Jan. 2019. [Baidu Scholar] 

3

E. Ebrahimzadeh, F. Blaabjerg, X. Wang et al., “Modeling and identification of harmonic instability problems in wind farms,” in Proceedings of the 8th Annual IEEE Energy Conversion Congress & Exposition (ECCE), Milwaukee, Wisconsin, USA, Sept. 2016, pp. 1-6. [Baidu Scholar] 

4

A. Shafiu, O. Anaya-Lara, G. Bathurst et al., “Aggregated wind turbine models for power system dynamic studies,” Wind Engineering, vol. 30, no. 3, pp. 171-185, May 2006. [Baidu Scholar] 

5

J. Bi, W. Du, and H. F. Wang, “An aggregation method of wind farms model for studying power system low frequency power oscillation,” in proceedings of 2017 2nd International Conference on Power and Renewable Energy (ICPRE), Chengdu, China, Sept. 2017, pp. 422-427. [Baidu Scholar] 

6

J. Sun, “Impedance-based stability criterion for grid-connected inverters,” IEEE Transactions on Power Electronics, vol. 26, no. 11, pp. 3075-3078, Apr. 2011. [Baidu Scholar] 

7

L. Fan and Z. Miao, “Nyquist-stability-criterion-based SSR explanation for Type-3 wind generators,” IEEE Transactions on Energy Conversion, vol. 27, no. 3, pp. 807-809, May 2012. [Baidu Scholar] 

8

F. Zhao, L. Wu, J. Zhang et al., “Suppression method of parallel-damping controller for DFIG,” Journal of Engineering, vol. 2017, no. 13, pp. 880-884, Oct. 2017. [Baidu Scholar] 

9

L. Sainz, M. Cheah-Mane, L. Monjo et al., “Positive-net-damping stability criterion in grid-connected VSC systems,” IEEE Journal of Emerging & Selected Topics in Power Electronics, vol. 5, no. 4, pp. 1499-1512, May 2017. [Baidu Scholar] 

10

X. Xie, H. Liu, J. He et al., “On new oscillation issues of power systems,” Proceedings of the CSEE, vol. 38, no. 10, pp. 2821-2828, May 2018. [Baidu Scholar] 

11

S. Shah, P. Koralewicz, V. Gevorgian et al., “Large-signal impedance-based modeling and mitigation of resonance of converter-grid systems,” IEEE Transactions on Sustainable Energy, vol. 10, no. 3, pp. 1439-1449, Jul. 2019. [Baidu Scholar] 

12

Y. Xu, Z. Gu, and K. Sun, “Characterization of subsynchronous oscillation with wind farms using describing function and generalized Nyquist criterion,” IEEE Transactions on Power Systems, vol. 35, no. 4, pp. 2783-2793, Jul. 2020. [Baidu Scholar] 

13

T. Wu, Q. Jiang, J. Shair et al., “Inclusion of current limiter nonlinearity in the characteristic analysis of sustained subsynchronous oscillations in grid-connected PMSGs,” IEEE Transactions on Energy Conversion, vol. 36, no. 3, pp. 2416-2426, Sep. 2021. [Baidu Scholar] 

14

A. Xue, Y. Wang, X. Fu et al., “Analysis of switched subsynchronous oscillation and its non-smooth bifurcation characteristics in doubly-fed system with SVG,” Power System Technology, vol. 45, no. 3, pp. 918-926, Mar. 2021. [Baidu Scholar] 

15

G. Xing, Y. Min, L. Chen et al., “Limit induced bifurcation of grid-connected VSC caused by current limit,” IEEE Transactions on Power Systems, vol. 36, no. 3, pp. 2717-2720, May 2021. [Baidu Scholar] 

16

A. Swami and J. M. Mendel, “Cumulant-based approach to the harmonic retrieval problem,” in Proceedings of 1988 International Conference on Acoustics, Speech, and Signal Processing (ICASSP 88), New York, USA, Apr. 1988, pp. 2264-2265. [Baidu Scholar] 

17

K. Lii and M. Rosenblatt, “Deconvolution and estimation of transfer function phase and coefficients for nonGaussian linear processes,” The Annals of Statistics, vol. 10, no. 4, pp. 1195-1208, Dec. 1982. [Baidu Scholar] 

18

G. Zhang, T. Shi, and S. Yang, “A method for extracting mechanical faults features based on higher-order statistics,” Journal of Huazhong University of Science and Technology, vol. 27, no. 3, pp. 6-8, Mar. 1999. [Baidu Scholar] 

19

M. Rosenblatt and J. W. V. Ness, “Estimation of the bispectrum,” The Annals of Mathematical Statistics, vol. 36, no. 4, pp. 1120-1136, Aug. 1965. [Baidu Scholar] 

20

D. R. Brillinger, “An introduction to polyspectra,” The Annals of Mathematical Statistics, vol. 36, no. 5, pp. 1351-1374, Oct. 1965. [Baidu Scholar] 

21

A. M. Tekalp and A. T. Erdem, “Higher-order spectrum factorization in one and two dimensions with applications in signal modeling and nonminimum phase system identification,” IEEE Transactions on Acoustics Speech & Signal Processing, vol. 37, no. 10, pp. 1537-1549, Oct. 1989. [Baidu Scholar] 

22

T. W. S. Chow and G. Fei, “Three phase induction machines asymmetrical faults identification using bispectrum,” IEEE Transactions on Energy Conversion, vol. 10, no. 4, pp. 688-693, Dec. 1995. [Baidu Scholar] 

23

B. Jang, C. Shin, E. J. Powers et al., “Machine fault detection using bicoherence spectra,” in Proceedings of the 21st IEEE Instrumentation and Measurement Technology Conference, Como, Italy, May 2004, pp. 1661-1666. [Baidu Scholar] 

24

Y. Jiang, C. Tang, X. Zhang et al., “A novel rolling bearing defect detection method based on bispectrum analysis and cloud model-improved EEMD,” IEEE Access, vol. 8, pp. 24323-24333, Jan. 2020. [Baidu Scholar] 

25

X. Chen, W. Xu, Y. Liu et al., “Bearing corrosion failure diagnosis of doubly fed induction generator in wind turbines based on stator current analysis,” IEEE Transactions on Industrial Electronics, vol. 67, no. 5, pp. 3419-3430, May 2020. [Baidu Scholar] 

26

J. R. Stack, R. G. Harley, and T. G. Habetler, “An amplitude modulation detector for fault diagnosis in rolling element bearings,” IEEE Transactions on Industrial Electronics, vol. 51, no. 5, pp. 1097-1102, Oct. 2004. [Baidu Scholar] 

27

Y. Li, Z. Han, and Z. Wang, “Research on a signal separation method based on Vold-Kalman filter of improved adaptive instantaneous frequency estimation,” IEEE Access, vol. 8, pp. 112170-112189, Jun. 2020. [Baidu Scholar] 

28

J. H. Taylor, Describing Functions. Hoboken: Wiley Encyclopedia of Electrical and Electronics Engineering, 2001. [Baidu Scholar] 

29

Z. An, C. Shen, Z. Zheng et al., “Assessment method for equivalent models of wind farms based on direct-driven wind generators considering randomness,” Proceedings of the CSEE, vol. 38, no. 22, pp. 6511-6519, Nov. 2018. [Baidu Scholar] 

30

K. M. Alawasa, Y. A. R. I. Mohamed, and W. Xu, “Active mitigation of subsynchronous interactions between PWM voltage-source converters and power networks,” IEEE Transactions on Power Electronics, vol. 29, no. 1, pp. 121-134, Mar. 2013. [Baidu Scholar] 

31

K. M. Alawasa, Y. A. R. I. Mohamed, and W. Xu, “Modeling, analysis, and suppression of the impact of full-scale wind-power converters on subsynchronous damping,” IEEE Systems Journal, vol. 7, no. 4, pp. 700-712, Jan. 2013. [Baidu Scholar] 

32

W. Song and A. Q. Huang, “Fault-tolerant design and control strategy for cascaded H-bridge multilevel converter-based STATCOM,” IEEE Transactions on Industrial Electronics, vol. 57, no. 8, pp. 2700-2708, Aug. 2010. [Baidu Scholar] 

33

X. Hou, Y. Sun, H. Han et al., “A general decentralized control scheme for medium-/high-voltage cascaded STATCOM,” IEEE Transactions on Power Systems, vol. 33, no. 6, pp. 7296-7300, Nov. 2018. [Baidu Scholar] 

34

C. L. Nikias and M. R. Raghuveer, “Bispectrum estimation: a digital signal processing framework,” Proceedings of the IEEE, vol. 75, no. 7, pp. 869-891, Jul. 1987. [Baidu Scholar] 

35

L. A. Pflug, G. E. Ioup, J. W. Ioup et al., “Properties of higher-order correlations and spectra for bandlimited, deterministic transients,” The Journal of the Acoustical Society of America, vol. 91, no. 2, pp. 975-988, Feb. 1992. [Baidu Scholar] 

36

P. A. M. Dirac, “The quantum theory of the emission and absorption of radiation,” Proceedings of the Royal Society of London, vol. 114, no. 767, pp. 243-265, Mar. 1927. [Baidu Scholar] 

37

M. A. A. S. Choudhury, S. L. Shah, and N. F. Thornhill, “Diagnosis of poor control-loop performance using higher-order statistics,” Automatica, vol. 40, no. 10, pp. 1719-1728, Oct. 2004. [Baidu Scholar] 

38

Quality of Electric Energy Supply – Harmonics in Public Supply Network, Chinese National Standard 14549-1993. [Baidu Scholar] 

39

K. L. Siu, J. M. Ahn, K. Ju et al., “Statistical approach to quantify the presence of phase coupling using the bispectrum,” IEEE Transactions on Biomedical Engineering, vol. 55, no. 5, pp. 1512-1520, Apr. 2008. [Baidu Scholar] 

40

Z. Yanbing, L. Yibing, A. Hongwen et al., “Fault recognition of large steam turbine based on higher order spectral features of vibration signals,” in Proceedings of 2011 IEEE International Conference on Mechatronics and Automation, Beijing, China, Aug. 2011, pp. 1572-1577. [Baidu Scholar] 

41

C. Qiu, Z. Xue, and P. Shao, “An identification method for weak motor fault eigenfrequencies based on multitaper bispectrum,” in Proceedings of 2010 International Conference on Intelligent Control and Information Processing, Dalian, China, Aug. 2010, pp. 236-239. [Baidu Scholar]