Abstract
Oscillation has become one of the important problems faced by modern power grids. Multitypes of oscillations may occur simultaneously in the power system and the oscillation frequency span is extremely large. For signals with wideband oscillation modes, the signals in different frequency bands are first separated by a bandpass filter, and then the Improved Variational Mode Decomposition (IVMD) method with high noise robustness is used to extract each oscillating mode signal. Finally, the combinations of Hankel total least squares (HTLS) and adaptive neural network algorithm (Adaline ANN) is used to estimate the frequency, attenuation factor, amplitude and phase of lowfrequency oscillations. Furthermore, the introduction of Adaline neural network solves the problem that the mode amplitude and phase are difficult to determine after IVMD processing, so that the detection accuracy is improved. Simulation and case analysis show that this method can effectively distinguish and extract different types of oscillation modes in the signal, and accurately identify the information of each mode. The IVMDHTLSAdaline method can effectively identify signals that have experienced severe oscillations or noiselike signals with potential oscillations.
1. Introduction
With the largescale access of renewable energy, the widespread application of power electronic equipment, and the largescale interconnection of AC and DC in modern power grids, the resource allocation is optimized and the system reliability is improved; the weak links in the system increase and the antiinterference performance decreases. In recent years, new types of faults have continued to emerge, and security and stability issues have become increasingly prominent [13]. Oscillation is one of the issues that threatens the stable operation of the power system. Typical oscillations in current power systems include Lowfrequency oscillations (LFO) [46] caused by regional or interval weak damping and subsynchronous Oscillation (SSO) caused by series compensation capacitors or energy interaction between the power electronics equipment and generator set [710]. In addition, Supersynchronous Oscillation (SurSO) occasionally occurs with subsynchronous oscillation [11]. Normally, oscillation can be tested in four ways: mechanism, acoustics, electrics and electromagnetic [1219]. Here, we will focus on the electrical signal waves.
According to different damping, the oscillation signals of power system can be divided into two categories: (1) Maintain or diverge the oscillation signal (when the system is weakly damped or negatively damped); (2) Damped oscillation signal (when the system is in positive damping). The former occurs less frequently, but once it occurs, it will cause great harm to the system. The latter occurs more frequently and is easily masked by noise, also known as noiselike oscillation signals. The rapid and realtime extraction of modal information from the oscillating signal is helpful to the adjustment of power grid dispatching and control strategies. Before the oscillation, the mode identification can be found from the noiselike oscillation signal, and the potential oscillation mode of the system can also be found to provide an early warning for the system. Oscillation mode identification is the key of realtime efficient control and risk early warning of power systems.
Frequency, damping ratio, oscillation amplitude and phase are the key information of oscillation mode, as well as the key parameters to realize oscillation monitoring, early warning, control and protection. At present, the research methods of power system oscillation mainly include model analysis and measurement analysis. The model analysis method depends on the precise model of the system. Largescale systems and nonlinear power electronic devices will affect the accuracy of the model, leading to analysis errors [20]. Measurement analysis analyzes the actual measurement data of the system, and the results are close to the operating conditions. Combining measured data with signal processing technology is a common method for power system oscillation identification.
The Prony algorithm can directly extract the characteristic quantity of the signal to identify more accurate oscillation modal information [21, 22], but it is highly sensitive to noise. In actual engineering, the filtering algorithm is often combined with the Prony algorithm, and the dimensionality of the signal is reduced through the filtering algorithm, or the modal signal is extracted from it. Mean filtering (Average Filter, AF), empirical mode decomposition, autoregressive moving average algorithm, wavelet method and singular value decomposition are commonly used filtering algorithms, which have been applied to the identification of oscillation modes [2329]. However, the mean filtering can only reduce the noise interference and cannot remove the noise. The essence of wavelet denoising is to fit the signal in a specific frequency band. If there are multiple oscillation modes with similar frequencies, the wavelet method cannot distinguish them. Methods such as Hilbert transform and singular value decomposition may produce large errors when processing low signaltonoise ratio signals. In addition, in the face of multitypes of coexisting panband oscillations with large oscillation frequency spans, the above methods often identify high amplitude oscillations during processing, and other oscillations may be considered as noise and ignored. How to carry out unified modal identification of panband oscillation signals and extract modal information of multiple types of oscillations is an important issue for modal identification of power systems.
Variational Mode Decomposition (VMD) algorithm can effectively separate modal signals and is not sensitive to noise [30]. Based on the Improved VMD (Improved VMD, IVMD) method, a method for identifying panband oscillations in complex power systems is proposed. First, different types of bandpass filters are used for separation. Secondly, the modal signals in each bandpass filtered signal are extracted using IVMD. Finally, the Hankel total least squares (HTLS) algorithm and the adaptive artificial neural network (Adaline Artificial Neural Network, Adaline ANN) are used to estimate the frequency, attenuation factor, amplitude, and phase of lowfrequency oscillations respectively, so as to achieve unified identification of panband oscillations. Meanwhile, the performance of the proposed algorithm is verified by simulation and measurement examples.
2. Modal identification framework for BroadBand oscillation
2.1. BroadBand oscillation signal
The oscillation may include multiple oscillation modal signals in different frequency bands and noise signals generated by measurement or the system itself. The original measurement signal with broadband oscillation can be expressed as ${y}_{0}\left(t\right)$:
In Eq. (1): ${y}_{Noise}\left(t\right)$ is the noisy signal, and ${y}_{n}\left(t\right)(n=\mathrm{1,2},\cdots ,N)$is the modal signals for $N$ different frequency bands, and:
In Eq. (2): ${A}_{n}\left(t\right)$ is the amplitude of the $n$th oscillation mode at time $t$; ${\theta}_{n}$ is the initial phase; ${D}_{n}$ and ${f}_{n}$ are the damping ratio and frequency of oscillation respectively.
2.2. Modal identification framework
The framework for modal oscillations identification in a broadband is shown in Fig. 1. The identification process is divided into the following three steps.
Step1. Multiple Band Pass Filters (BPF) is used to decompose the original signal into multiple filtered signals in different frequency bands, so as to separate the oscillation signals.
Step2. The IVMD algorithm is used to extract the oscillation modal signal from each BPF filtered signal.
Step3. The HTLSAdaline algorithm is used to identify the oscillation mode signals provided by IVMD, and obtain the information of each oscillation mode. The original signal is filtered through bandpass filtering, IVMD modal signal extraction and HTLSAdaline identification to obtain modal parameters.
Fig. 1Framework of the mode identification for broadband oscillation
3. Broadband oscillation mode identification base on VMDHLSAdaline
3.1. Bandpass filter design
There are mainly three types of oscillations that cause power systems accidents: lowfrequency oscillation, subsynchronous oscillation and supersynchronous oscillation. The lowfrequency oscillation is further divided into local mode and interval mode. The frequencies of each type of oscillations are in different ranges. Thus, this article designs four BPFs. A proper signal length selection can improve the accuracy of pattern recognition and ensure the rapidity of recognition. Because each oscillation frequency is different, it is obviously inappropriate to use the same length of data for identification. Table 1 shows the BPF filter bandwidth and identification sampling time corresponding to various oscillations in this paper.
Table 1Time lengths and frequency bands for different BPFs
Bandpass filter  Oscillation type  Sampling time / s  Oscillation frequency / Hz 
BPF1  LFO (interval modal)  10  0.21 
BPF2  LFO (local mode)  2  15 
BPF3  SSO  0.4  1050 
BPF4  SurSO  0.4  70110 
According to the analysis of the modal identification algorithm [31], if it is necessary to obtain the oscillation frequency and amplitude more accurately, the sampling time of the identification signal must be greater than 1 oscillation period. In order to obtain accurate damping ratio information, the sampling time needs to be longer than two oscillation periods. Therefore, for lowfrequency oscillations, in order to ensure rapid identification, the minimum bandpass frequency of the bandpass filter is 0.2 Hz, and the sampling time is selected to be twice the oscillation period, that is 10 s. For subsynchronous oscillation, the identification speed is considered comprehensively. The minimum oscillation frequency is 10 Hz as the benchmark and the sampling time is 4 times. The oscillation period is 0.4 s. The mechanism of supersynchronous oscillation determines that it appears in pairs with subsynchronous oscillation and the frequencies are complementary. The sampling time of the synchronous oscillation is the same as that of the subsynchronous oscillation, and is also 0.4 s.
3.2. Modal signal extraction based on improved VMD algorithm
3.2.1. VMD algorithm
The VMD algorithm is a new adaptive signal decomposition method proposed by Dragomiretskiy et al. in 2014. The target modal is solved by the inherent modal function [30]. In view of the accuracy and noise robustness of the VMD algorithm, this paper uses the VMD method to extract and separate modal signals. The variational problem corresponding to the VMD algorithm is to find the KIMF components with the smallest sum of the estimated bandwidth. The variational problem is transformed into an augmented Lagrange equation and the equation is solved by the Alternating Direction Method of Multipliers (ADMM) to obtain the solution of the modal function ${u}_{k}^{n+1}\left(t\right)$:
Similarly, the solution of the center frequency value of the modal component is:
where, $\left\{{\omega}_{k}\right\}=\{{\omega}_{1},{\omega}_{2},\cdots ,{\omega}_{K}\}$ and $\left\{{\omega}_{k}\right\}=\{{\omega}_{1},{\omega}_{2},\cdots ,{\omega}_{K}\}$ is the component and its center frequency respectively. The flow of VMD algorithm is as follows:
Step 1. Initialize $\left\{{\omega}_{k}\right\}=\{{\omega}_{1},{\omega}_{2},\cdots ,{\omega}_{K}\}$ and $n=0$.
Step 2. $n\leftarrow n+1$, update ${u}_{k}$ and ${\omega}_{k}$ according to Eq. (3) and Eq. (4).
Step 3. Update $\lambda $:
Step 4. Repeat steps 2 and 3 until the iteration stop condition is satisfied with $\sum _{k}{\Vert {\widehat{u}}_{k}^{n+1}{\widehat{u}}_{k}^{n}\Vert}_{2}^{2}/{\Vert {\widehat{u}}_{k}^{n}\Vert}_{2}^{2}<\epsilon $. End the loop and output the results to get $K$ modal components and their center frequencies.
3.2.2. Determination of VMD penalty factor based on PSO algorithm optimization
The penalty factor $\alpha $ in the VMD algorithm has a great impact on the decomposition results. The study found that the smaller the penalty parameter $\alpha $, the larger the bandwidth of each IMF (Intrinsic Mode Function) component, and vice versa [32]. Therefore, when using VMD to decompose the oscillation signal of the power system, it is very important to choose the appropriate penalty factor parameter $\alpha $. In this paper, genetic mutation particle swarm optimization is used to optimize the penalty parameters to obtain the optimal $\alpha $.
Particle swarm optimization is a global optimization algorithm proposed by Eberh and Kennedy et al. in 1995. This method is a swarm intelligent optimization algorithm, and it has the advantages of few parameters, easy adjustment, and easy to fall into a local optimum. In order to obtain the global optimal approximate solution [33], this paper introduces the idea of genetic algorithm mutation in particle swarm algorithm to construct a genetic mutation particle swarm algorithm.
Definition of genetic mutation particle swarm algorithm: In an $D$dimensional search space, label population $X$, and it is composed of m particles, $\mathbf{X}=[{\mathbf{x}}_{1},{\mathbf{x}}_{2},\cdots ,{\mathbf{x}}_{m}]$. The position of each particle in the search space can be represented by a $D$dimensional vector, that is, ${\mathbf{x}}_{i}=[{x}_{i1},{x}_{i2},\cdots ,{x}_{iD}]$. $D$ is the number of parameters to be optimized, and the moving speed of $i$th particle is ${\mathbf{v}}_{i}=[{v}_{i1},{v}_{i2},\cdots ,{v}_{iD}]$. The local extremum of particles is ${\mathbf{p}}_{i}=[{p}_{i1},{p}_{i2},\cdots ,{p}_{iD}]$; the global extremum of the population is ${\mathbf{G}}_{1}=[{g}_{1},{g}_{2},\cdots ,{g}_{D}]$, and the subglobal optimal value is ${\mathbf{G}}_{2}=[{g\text{'}}_{1},{g\text{'}}_{2},\cdots ,{g\text{'}}_{D}]$. The maximum individual optimal algebra is $\mathrm{m}\mathrm{a}\mathrm{x}Age$, and the mutation probability is$q$. In order to prevent particles from falling into the local optimum, it is necessary to record the individual optimal maintaining algebra of the particles during the iteration process. When the individual optimal maintaining algebra does not reach $\mathrm{m}\mathrm{a}\mathrm{x}Age$, each particle updates the position of the next generation by individual local extreme value and global extreme value. The formula is updated as:
In the formula, $\omega $ is the inertia weight; $\eta $ is a random number between [0, 1]; ${c}_{1}$ and ${c}_{2}$ are the learning factors respectively that represent the local search ability and the global search ability; $n$ is the number of iterations, and ${\mathbf{v}}_{i}$, ${\mathbf{p}}_{i}$, ${\mathbf{G}}_{1}$ and ${\mathbf{x}}_{i}$ are $D$dimensional vectors respectively. The determination of the inertia weight $\omega $ of the current number of iterations adopts the linear decreasing weight method proposed by Shi [34], and the formula is as follows:
In the formula, ${\omega}_{\mathrm{m}\mathrm{a}\mathrm{x}}$ and ${\omega}_{\mathrm{m}\mathrm{i}\mathrm{n}}$ are the maximum and minimum inertia weights; $n$ is the current number of iterations, and ${n}_{\mathrm{m}\mathrm{a}\mathrm{x}}$ is the defined maximum number of iterations. When the individual optimal retention algebra reaches $\mathrm{m}\mathrm{a}\mathrm{x}Age$, the genetic mutation operation is used to update the position and velocity of the particle to make it jump out of the local optimal. Selection of fitness function for genetic mutation particle swarm optimization algorithm: In the parameter optimization, the evaluation criterion of the decomposition effect of VMD method uses the concept of envelope entropy ${E}_{p}$ proposed by Tang Guiji et al. [35]. The envelope entropy of $x\left(j\right)$ time signal $c$ of length $N$ is defined as:
In the formula, $a\left(j\right)$ is the envelope signal of $x\left(j\right)$ after Hilbert demodulation, and $j=\mathrm{1,2},\cdots ,N$. ${p}_{j}$ is the result of normalizing $a\left(j\right)$. Normalization not only avoids the influence of different envelope amplitudes of IMF components, but also reduces the interference of weak noise. ${E}_{p}$ is obtained according to the information entropy calculation rules. This article measures the decomposition effect of VMD according to ${E}_{p}$.
The VMD method is used to decompose the BPF filtered oscillation signal. When the component contains more noise, the sparseness of the component signal is weak and the envelope entropy is large. On the contrary, when a regular oscillation signal appears in the component, the signal will have strong sparseness, and the calculated envelope entropy is small at this time. Therefore, under the influence of parameter $\alpha $, the minimum entropy ${E}_{p}$ of the $K$ components is selected as the local minimum entropy $\mathrm{m}\mathrm{i}\mathrm{n}{E}_{p}$. The component corresponding to the minimum entropy value contains rich feature information. The local minimum entropy is used as part of the fitness function of the entire search process to find the parameter ${\alpha}_{0}$ corresponding to the global optimal component. Through the above analysis of parameter $\alpha $, a proper $\alpha $ will reduce the iterations of VMD, that is, the VMD method has a high decomposition efficiency. Therefore, it is necessary to achieve the highest decomposition efficiency in the case of the best decomposition effect. This article builds the fitness function based on $\mathrm{m}\mathrm{i}\mathrm{n}{E}_{p}$, add $time$ (iterations) as follows $\mathrm{m}\mathrm{i}\mathrm{n}F=\mathrm{m}\mathrm{i}\mathrm{n}{E}_{p}+\beta \cdot time$, where $\beta $ is the quantization factor of the fitness function.
In this paper, the number of modal components $K$ and penalty factor $A$ in VMD decomposition are set as model hyper parameters. In the process of VMD optimization by PSO, the number of particle is set as 30, and the maximum number of iterations is set as 500.The learning factor is set as ${c}_{1}={c}_{2}=$2; the velocity inertia factor is set as $w=$ 0.8, and the velocity coefficient is setas $\lambda =$ 1. The maximum and minimum inertia weights are set as ${\omega}_{\mathrm{m}\mathrm{a}\mathrm{x}}=$ 0.9 and ${\omega}_{\mathrm{m}\mathrm{i}\mathrm{n}}=$ 0.1 respectively, and the quantization factor of the fitness function is set as $\beta =$1/1000. The parameter optimization process based on genetic mutation particle swarm optimization algorithm is shown in Fig. 2.
Fig. 2Parameter optimization process based on genetic mulation particle swarm optimization
The BPF signal is decomposed by improved VMD to obtain the oscillation mode ${u}_{n}\left(t\right)$ extracted by VMD. When the accuracy of BPF and VMD is high enough, for modal signal ${u}_{n}\left(t\right)$, there is ${u}_{n}\left(t\right)\approx {y}_{n}\left(t\right)$.
Therefore, through BPF and IVMD, all oscillation modal signals can be separated from the original signal. The convergence process of VMD optimization using PSO is shown in Fig. 3.
Fig. 3Convergence process of fitness function in VMD optimization by PSO
3.3. VMDHTLS algorithm for frequency and attenuation factor
After the modal signals ${u}_{n}\left(t\right)$ are obtained by the IVMD decomposition, the HTLS algorithm is used to identify the modal parameters, such as the oscillation frequency and attenuation factor. HTLS algorithm is a subspace rotation invariant method, which has high computing efficiency and strong antinoise ability. Its calculation steps are described in [29]. The main idea is to construct a Hankel matrix using the sampled signal, and perform Vander Mang decomposition on it. Using the translationinvariant characteristic of the van der Mun matrix to construct the equivalent relationship, and the eigenvalues of the oscillating modes are obtained. The main idea of IVMDHTLS is: IVMD is used to decompose the BPF filtered sequence, and then the HTLS algorithm is used to calculate the oscillation frequency and attenuation factor for each component after decomposition. Because the FOMCHTLS algorithm cannot give the amplitude and phase of the original signal $y\left(n\right)$, and the information of each mode is incomplete. It is not conducive to the reconstruction of the signal and the quantitative evaluation of the algorithm. Therefore, this paper introduces the Adaline God network to oscillate modal information (Amplitude and Phase).
3.4. Adaline neural network solves amplitude and phase
3.4.1. Adaptive linear neural network
Adaptive linear (Adaline) neural network is a neuron model proposed originally by Widrow and Hoff [36]. It is widely used in signal processing and other fields.
Fig. 4The principle of adaptive linear neural network
In Fig. 4, ${x}_{1k}$, ${x}_{2k}$,⋯, ${x}_{nk}$ are the $n$ input signals of the adaptive linear neural network at time$k$.The input signal vector form is expressed as ${\mathbf{X}}_{ik}=[{x}_{1k},{x}_{2k},\cdots ,{x}_{nk}{]}^{T}$, and this is often called the input pattern vector of Adaline neural network. The weight vector corresponding to each group of input signals is ${\mathbf{W}}_{ik}=[{w}_{1k},{w}_{2k},\cdots ,{w}_{nk}{]}^{T}$. The Adaline neural network output is:
Let $A$ be the ideal response signal, and define the error function as:
The working process of Adaline neural network is as follows [36]: the ideal response signal $y\left(k\right)$ is compared with the output signal $\widehat{y}\left(k\right)$ of the neural network to obtain different $e\left(k\right)$. Feed$e\left(k\right)$into the learning rules, and adjust the weight vector according to the learning rules to make $\widehat{y}\left(k\right)$ and $y\left(k\right)$ consistent.
The learning rule of Adaline neural network is WidrowHoff rule, which is the least square error algorithm (LMS). The rule weight vector adjustment expression is:
In the formula, $\eta $ is the learning rate of the Adaline neural network, $\eta \in \left(\mathrm{0,1}\right)$. Its value directly affects the weight vector adjustment accuracy and the convergence velocity.
3.4.2. Solution of oscillation modes by Adaline neural network
The specific steps of Adaline neural network to solve the amplitude and phase are as follows, and a known lowfrequency oscillation discrete sampling signal model is established:
When the attenuation factor and frequency are known, Eq. (12) can be written as:
$=\sum _{i=1}^{M}\left[{A}_{i}\mathrm{c}\mathrm{o}\mathrm{s}{\theta}_{i}{e}^{{\alpha}_{i}n\mathrm{\Delta}t}\mathrm{c}\mathrm{o}\mathrm{s}\right(2\pi {f}_{i}n\mathrm{\Delta}t){A}_{i}\mathrm{s}\mathrm{i}\mathrm{n}{\theta}_{i}{e}^{{\alpha}_{i}n\mathrm{\Delta}t}\mathrm{s}\mathrm{i}\mathrm{n}\left(2\pi {f}_{i}n\mathrm{\Delta}t\right)]$
$\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}=\sum _{i=1}^{M}[{p}_{i}{e}^{{\alpha}_{i}n\mathrm{\Delta}t}\mathrm{c}\mathrm{o}\mathrm{s}(2\pi {f}_{i}n\mathrm{\Delta}t){q}_{i}{e}^{{\alpha}_{i}n\mathrm{\Delta}t}\mathrm{s}\mathrm{i}\mathrm{n}\left(2\pi {f}_{i}n\mathrm{\Delta}t\right)].$
In the formula, ${p}_{i}={A}_{i}\mathrm{c}\mathrm{o}\mathrm{s}{\theta}_{i}$ and ${q}_{i}={A}_{i}\mathrm{s}\mathrm{i}\mathrm{n}{\theta}_{i}$. The matrix expression of Eq. (13) is:
Among them:
$\mathbf{S}=\left[\begin{array}{cccc}\mathrm{s}\mathrm{i}\mathrm{n}\left(2\pi {f}_{1}\mathrm{\Delta}t\right)& \mathrm{s}\mathrm{i}\mathrm{n}\left(2\pi {f}_{1}2\mathrm{\Delta}t\right)& \cdots & \mathrm{s}\mathrm{i}\mathrm{n}\left(2\pi {f}_{1}N\mathrm{\Delta}t\right)\\ \mathrm{s}\mathrm{i}\mathrm{n}\left(2\pi {f}_{2}\mathrm{\Delta}t\right)& \mathrm{s}\mathrm{i}\mathrm{n}\left(2\pi {f}_{2}2\mathrm{\Delta}t\right)& \cdots & \mathrm{s}\mathrm{i}\mathrm{n}\left(2\pi {f}_{2}N\mathrm{\Delta}t\right)\\ \vdots & \vdots & \vdots & \vdots \\ \mathrm{s}\mathrm{i}\mathrm{n}\left(2\pi {f}_{Q}\mathrm{\Delta}t\right)& \mathrm{s}\mathrm{i}\mathrm{n}\left(2\pi {f}_{Q}2\mathrm{\Delta}t\right)& \cdots & \mathrm{s}\mathrm{i}\mathrm{n}\left(2\pi {f}_{Q}N\mathrm{\Delta}t\right)\end{array}\right].$
Similarly, define the error function:
In the formula, $x\left(n\right)$ is the actual sample, and $\widehat{x}\left(n\right)$ is the output of the neural network. Define the performance indicators as:
Because the attenuation factor and frequency are known, $\mathbf{p}$ and $\mathbf{q}$ are unknown in Eq. (14), $\mathbf{C}$ and $\mathbf{S}$ are the input vectors of the neural network. According to the training principle of the steepest descent method, the weight vectors p and q are adjusted to:
When the Adaline neural network training is completed, the amplitude and phase of the oscillation mode are solved from the obtained weight vector and Eq. (19):
In this paper, the number of neurons in Adaline neural network is equal to the number $K$ of modal components after VMD decomposition. The activation function of neurons is a constant function, that is, $f=$ 1. The learning rule is the minimum mean square error (LMS) criterion, and the learning rate $\eta =$ 0.0015. The maximum number of network iteration is 5000, and the network iteration stops when the error $\delta $ satisfies $\delta <$ 0.0001.
4. Simulation and analysis of actual examples
4.1. Simulation analysis
In order to verify the effectiveness of the method in this paper, an oscillating signal ${y}_{test}$ is constructed by simulation:
In the formula: ${y}_{LFO1}$ is the lowfrequency oscillation interval mode;${y}_{LFO2}$ is the lowfrequency oscillation local mode; ${y}_{SSO1}$ is the subsynchronous oscillation and paired with the supersynchronous oscillation signal ${y}_{SurSO}$; ${y}_{SSO2}$ is another independent subsynchronous oscillation, and ${y}_{Noise}$ is the white noise. In line with the real situation, the test signal satisfies the following conditions and assumptions: (1) All oscillation frequencies are randomly selected within the frequency band of the oscillation type; (2) Low frequency oscillation is the main mode of oscillation; the amplitude is higher than the subsynchronous oscillation, and the local mode of low frequency oscillation is higher than the interval mode; (3) In the pair of subsynchronous and supersynchronous oscillations, the amplitude of subsynchronous oscillation mode is slightly higher than the supersynchronous mode. Based on the above assumptions, the specific parameters of each mode of the test signal are finally selected as:
That is, the test signal ${y}_{test}$ contains 5 oscillating signals with different frequencies and a white noise signal with amplitude of 0.16. The frequency band of the test signal oscillation ${y}_{test}$ is 0.6393.42 Hz. Each test signal is equal amplitude oscillation. The time domain form of the test signal ${y}_{test}$ is shown in Fig. 5. Using the proposed method after the constructed test signal ${y}_{test}$ is extracted by bandpass filtering and IVMD, the oscillation modal signals in each frequency band are shown in Fig. 6. It can be seen from Fig. 6 that the method proposed in this paper can accurately distinguish the modalities of different frequency bands and effectively extract all corresponding modal signals.
Fig. 5Test signal ytest
Fig. 6Mode signals extracted by IVMD
The modal signals are extracted for each IVMD obtained in Fig. 6, and the modal identification is performed through the HTLSAdaline algorithm to obtain the parameter information of the corresponding modal. The identification results of this method compared with those of other methods is shown in Table 2. The selected comparison methods are: (1) the classic EMDbased HTLS algorithm (EMDHTLS); (2) the VMDbased HTLS algorithm (VMDHTLS); (3) the proposed method is based on the algorithm of IVMDHTLS and Adaline (IVMDHTLSAdaline).
Table 2Mode identification of test signal
Test signal  EMDHTLS  VMDHTLS  Proposed method  
Modal  Freq / Hz  Amp  Phase / rad  Fre  Amp  Phase / rand  Freq / Hz  Amp  Phase / rad  Freq / Hz  Amp  Phase / rad 
error  
${y}_{LPO1}$  0.65  2.35  1.0472  0.6648  2.2517  1.1815  0.6619  2.3841  0.9122  0.6512  2.3492  1.0488 
2.029%  15.637%  12.82%  1.044%  5.199%  12.89%  0.264%  0.945%  0.153%  
${y}_{LPO2}$  2.05  3.00  0.6283  2.0641  3.1125  0.4868  2.4833  2.0855  0.7318  2.0525  3.0263  0.6338 
4.321%  5.670%  22.52%  2.762%  4.775%  16.473%  0.233%  1.315%  0.875%  
${y}_{SSO1}$  25.56  0.58  0.5236  26.992  0.6233  0.6670  25.173  0.5922  0.5961  25.612  0.5842  0.5168 
5.384%  23.122%  27.38%  2.231%  8.612%  13.846%  0.158%  4.122%  1.299%  
${y}_{SSO2}$  22.84  0.25  0.7854  25.237  0.2747  1.0141  22.157  0.2419  0.7741  22.774  0.2484  0.7880 
10.13%  15.772%  29.11%  6.002%  5.409%  1.438%  0.487%  0.727%  0.331%  
${y}_{SurSO}$  92.32  0.42  1.0472  –  –  –  93.451  0.4438  1.1187  92.259  0.4291  1.0482 
1.104%  7.437%  6.828%  0.215%  2.843%  0.096% 
It can be seen from Table 2 that when the test signal oscillates in multiple frequency bands, the proposed method can effectively identify all oscillation modes. The maximum error of the oscillation frequency identification result is 1.55 %, with an average error of 0.67 %, and that of the oscillation amplitude identification result is 9.09 %, with an average error of 2.54 %. In comparison, the EMDHTLS method has the worst frequency identification and amplitude identification results, and the average error is also the highest. The VMDHTLS method has a good identification effect on the dominant oscillation mode (${y}_{LFO2}$) with the highest amplitude. The frequency and amplitude identification results are close to those of this method, but the phase identification results are not as good as those of this method. For the submodes ${y}_{LFO1}$ and ${y}_{SSO1}$, the VMDHTLS frequency identification effect is close to the proposed method, but the amplitude and phase identification results are far worse than it. In particular, for mode ${y}_{SSO2}$ with low amplitude and high frequency, the recognition results of EEMDHTLS and VMDHTLS are poor. For the supersynchronous oscillation mode ${y}_{SurSO}$ with a frequency of 93.42, the EEMDHTLS method failed to identify it, and the result of VMDHTLS is inconsistent with the actual one. Therefore, the method proposed in this paper has obvious advantages in identifying multimode coexisting panband oscillations.
4.2. Actual study data
In order to prove the effectiveness of the proposed method, the actual oscillation data of Hunan Power Grid was selected and the oscillation mode identification analysis was carried out. The oscillation event is as follows: On June 24, 2018, a lowfrequency oscillation occurred in a thermal power plant in Hunan Power Grid. The system started to oscillate at low frequency in the 120th second. After 30 s, the system quickly started to emit an alarm. The system maintained an alarm for about 180 s. In 60120 s before the oscillation, the system has a noiselike oscillation, and the system did not issue an alarm at this time. The range of the oscillation signal and noiselike signal has been marked in the figure. The system sampling frequency is 100 Hz. The oscillation signal at 165195 s is selected, and the modal identification of this method is adopted. The original signal used for identification is shown in Fig. 7.
It can be seen from Fig. 6 that the amplitude of the oscillation signal is high, and the highest frequency fluctuation exceeds 100 MW, which is higher than 30 % of the system output power. At this time, the signal is affected by a certain degree of noise. For the original oscillation signal data shown in Fig. 7, the modal signals are extracted by the proposed algorithm as shown in Fig. 8. This method decomposes three modes from the oscillation signal data, which are the main modes of the system oscillation. Among the three modes, the first mode has the highest amplitude and is the dominant oscillation mode. The extracted modal signals are linearly superimposed, and the reconstructed signals are shown in Fig. 8.
Fig. 7Original oscillation data used for mode identification
Fig. 8Mode signals extracted from the original oscillation data
Fig. 9Reconstructed signal based on mode signals from IVMD
Comparing Fig. 9 with Fig. 7, it is not difficult to find that the reconstructed signal is basically the same as the original signal waveform. That is, the three modal signals separated cover the main oscillations of the system. The method in this paper can effectively extract all the main oscillation modal signals. Using the Prony algorithm to perform parameter identification on the modal signals in Fig. 6, the information of the main oscillation modes can be obtained. Table 3 compares the results of this method with the MFProny and EMDProny parameter identification results. It can be seen from the results that for the oscillation signal, both the proposed method MFProny or EMDProny can identify the mode with the largest amplitude. In the three methods, the MFProny frequency identification error is the highest, and the results of this method are similar to the EMDProny method. During the severe oscillation, the dominant modal frequency is about 1.48 Hz, and the minor dominant modal frequency is about 2.02 Hz.
Table 3Mode identification result and comparison of original oscillation data
Method  Modal  Frequency / Hz  Amplitude / MW  Damping  Phase / rad 
EEMDHTLS  1  1.44  28.2  –1.18  2.6539 
2  1.61  13.7  –1.59  0.2355  
3  2.32  8.5  –2.31  2.9813  
VMDHTLS  1  1.48  30.6  –0.52  3.1348 
2  1.51  16.3  –1.01  0.1653  
3  1.08  8.7  –1.03  2.1981  
IVMDHTLSAdaline  1  1.48  32.7  –0.94  3.0273 
2  2.02  17.9  –1.40  0.0985  
3  0.49  3.2  0.18  2.3499 
5. Conclusions
A VMDbased method for power system panband oscillation signal extraction and modal identification is proposed. This method can effectively extract lowfrequency oscillation, subsynchronous oscillation, and supersynchronous oscillation of multitype and panband oscillations from power system operating data. Based on the oscillation discrimination, modal extraction and parameter identification, the identification results plays an important role in the analysis and control of system dynamic stability, as well as the identification and location of oscillation sources. It is helpful for early warning and timely suppression of system oscillation.
References

Kang Jie, Liu Li, Shao Yupei, Ma Qinggang Nonstationary signal decomposition approach for harmonic responses detection in operational modal analysis. Computers and Structures, Vol. 242, Issue 1, 2021, p. 106377.

David G. L., Jose R. R. H., Vicente V. H., Juan O. G., et al. Harmonic PMU and fuzzy logic for online detection of shortcircuited turns in transformers. Electric Power Systems Research, Vol. 190, Issue 1, 2021, p. 106862.

Duan Guizhong, Qin Wenping, Lu Ruipeng, et al. Static voltage stability analysis considering the wind power and uncertainty of load. Power System Protection and Control, Vol. 46, Issue 12, 2018, p. 108114.

Zhang Xinran, Lu Chao, Liu Shichao, et al. A review on widearea damping control to restrain interarea low frequency oscillation for largescale power systems with increasing renewable generation. Renewable and Sustainable Energy Reviews, Vol. 57, 2016, p. 4558.

Li Shenghu, Sun Qi, ShiXuemei, et al. Suppression of weakly damped lowfrequency modes of wind power system based on regional pole placement. Power System Protection and Control, Vol. 45, Issue 20, 2017, p. 1420.

Gao Haixiang, WuShuangxi, Miao Lu, et al. Overview of reasons for generatorinduced power oscillations& its suppression measures. Smart Power, Vol. 46, Issue 7, 2018, p. 4955.

Suriyaarachchi D. H. R., Annakkage U. D., Karawita C., et al. A Procedure to study subsynchronous interactions in wind integrated power systems. IEEE Transactions on Power Systems, Vol. 28, Issue 1, 2013, p. 377384.

Karaagac U., Faried S. O., Mahseredjian J. Coordinated control of wind energy conversion systems for mitigating sub synchronous interaction in DFIGbased wind farms. IEEE Transactions on Smart Grid, Vol. 5, Issue 5, 2014, p. 24402449.

Cui Yahui, Zhang Junjie, Zhao Zongbin Sub synchronous oscillation torsional vibration model and simulation for a supercritical 600 MW unit with high voltage direct current transmission. Thermal Power Generation, Vol. 45, Issue 6, 2016, p. 106110.

Zhang Fan, Zhang Donghui, Liu Yongjun Risk assessment and suppression method study on subsynchronous oscillation resulted from Luxi mixed backtoback HVDC project. Smart Power, Vol. 45, Issue 4, 2017, p. 3033.

Xu Yanhui, Cao Yuping Research on mechanism of sub/supsynchronous oscillation caused by GSC controller of directdrive permanent magnetic synchronous generator. Power System Technology, Vol. 42, Issue 5, 2018, p. 15561564.

Kun Liang, Xiaogong Lin, Yu Chen, Juan Li, Fuguang Ding Adaptive sliding mode output feedback control for dynamic positioning ships with input saturation. Ocean Engineering, Vol. 206, Issue 15, 2020, p. 107245.

Alejandro Garces Smallsignal stability in island residential microgrids considering droop controls and multiple scenarios of generation. Electric Power Systems Research, Vol. 185, 2020, p. 106371.

Ayse Nihan Basmaci Characteristics of electromagnetic wave propagation in a segmented photonic waveguide. Journal of Optoelectronics and Advanced Materials, Vol. 22, Issues 910, 2020, p. 452460.

Ayse Nihan Basmaci Filiz Solution of Helmholtz equation using finite differences method in wires have different properties along Xaxis. Journal of Electrical Engineering, Vol. 6, 2018, p. 206211.

Zhineng Zhang, Ling Zheng 1D numerical study of nonlinear propagation of finite amplitude waves in traveling wave tubes with varying cross section. International Journal of Acoustics and Vibration, Vol. 25, Issue 1, 2020, p. 8895.

Eldad J. A., Neeshtha D. B., Giuseppe C. G., Touvia M. Sound scattering by an elastic spherical shell and its cancellation using a multipole approach. Archives of Acoustics, Vol. 42, Issue 4, 2017, p. 697705.

Mingsian Bai R., Jia Hong Lin, Kwan Liang Liu Optimized microphone deployment for nearfield acoustic holography: to be, or not to be random, that is the question. Journal of Sound and Vibration, Vol. 329, Issue 14, 2010, p. 28092824.

Mraa Abad, Ahmadi H., Moosavian A., Khazaee M., Mohammadi M. Discrete wavelet transform and artificial neural network for gearbox fault detection based on acoustic signals. Journal of Vibroengineering, Vol. 15, Issue 1, 2013, p. 459463.

Wang Henan, Zheng Chao, Ren Jie Review on mechanism and analysis methods of low frequency oscillation in power system. Advanced Materials Research, Vol. 986, 2014, p. 20102013.

Wadduwage D. P., Annakkage U. D., Narendrak K. Identification of dominant lowfrequency modes inringdown oscillations using multiple Prony models. IET Generation, Transmission and Distribution, Vol. 9, Issue 15, 2015, p. 22062214.

Jin Tao, Liu Dui Power grid low frequency oscillation recognition based on advanced Prony algorithm with improved denoising feature. Electric Machines and Control, Vol. 21, Issue 5, 2017, p. 3341.

Li Kuan, Li Xingyuan, Hu Nan ESPRIT analysis of sub synchronous oscillation based on the empirical mode decomposition selfadaptive filter. Power System Protection and Control, Vol. 40, Issue 13, 2012, p. 1823.

Zhang Yulin, Chen Hongwei Parameter identification of harmonics and interharmonics based on CEEMDWPT and Prony algorithm. Power System Protection and Control, Vol. 46, Issue 12, 2018, p. 115121.

Wang Shen, Huang Songling, Wang Qing Mode identification of broadband Lamb wave signal with squeezed wavelet transform. Applied Acoustics, Vol. 125, 2017, p. 91101.

Ren Zihui, Liu Haoyue, Xu Jinxia Power quality disturbance analysis based on wavelet transform and improved Prony method. Power System Protection and Control, Vol. 44, Issue 9, 2016, p. 122128.

Ye Y., Yuanzhang S., Lin C. Power system low frequency oscillation monitoring and analysis based on multisignal online identification. Science China (Technological Sciences), Vol. 53, Issue 9, 2010, p. 25892596.

Jiang Ping, Shi Hao, Wu Xi Localizing disturbance source of power system forced oscillation caused by wind power fluctuation. Electric Power Engineering Technology, Vol. 37, Issue 5, 2018, p. 16.

Xu Qi, Xu Jian, Shi Weil The forced oscillation analysis of wind integrated power systems based on singular value decomposition. Proceedings of the CSEE, Vol. 36, Issue 18, 2016, p. 48174827.

Dragomiretskiy K., Zosso D. Variational mode decomposition. IEEE Transactions on Signal Processing, Vol. 62, Issue 3, 2014, p. 531544.

Hauer J. F., Demeure C. J., Scharf L. L. Initial results in Prony analysis of power system response signals. IEEE Transactions on Power Systems, Vol. 5, Issue 1, 1990, p. 8089.

Tang Guiji, Wang Xiaolong Parameter optimized variational mode decomposition method with application to incipient fault diagnosis of rolling bearing. Journal of Xi’an Jiaotong University, Vol. 49, Issue 5, 2015, p. 7381.

Ghodousian A., Parvari M. R. A modified PSO algorithm for linear optimization problem subject to the generalized fuzzy relational inequalities with fuzzy constraints. Information Sciences, Vol. 418, 2017, p. 317345.

Liu Yanmin Research and Application of Particle Swarm Optimization. Ph.D. Dissertation, Shandong Normal University, China, 2011.

Tang Guiji, Wang Xiaolong Application of variational modal decomposition method for optimizing parameters inearly fault diagnosis of rolling bearings. Journal of Xi’an Jiaotong University, Vol. 49, Issue 5, 2015, p. 7381.

Papy J. M., De Lathauwer L., Van Huffel S. Common pole estimation in multichannel exponential data modeling. Signal Processing, Vol. 86, Issue 4, 2006, p. 846858.