A Periodogram-Based Method to Identify Forced and Natural Oscillations Using PMUs
Qinye Tang Xiaozhe Wang






A Periodogram-Based Method to Identify Forced and Natural Oscillations Using PMUs
This work is supported by the Fonds de Recherche du Québec - Nature et technologies under Grant FRQ-NT NC-253053 and by Mitacs Globalink Research Internship program.
    Qinye Tang

Zhejiang University

Hangzhou, China

Email:qtang@zju.edu.cn 
Xiaozhe Wang
McGill University

Montreal, QC Canada

Email: xiaozhe.wang2@mcgill.ca

    December 30, 2023
=======================================================================================================================================================================================================================================================




Sustained oscillations in power systems are dangerous. There are various mechanisms, for instance, limit cycle and forced oscillation, that may lead to sustained oscillations, which nevertheless are hard to differentiate. In this paper, a novel periodogram-based method to distinguish different oscillation mechanism is proposed, which can quantitatively extract essential signatures of different mechanisms from power spectral density. Numerical study shows that the proposed method can accurately distinguish different oscillation mechanisms even when the forced oscillation frequency is close to the natural frequency. 





rect=[draw,rectangle,fill=red!2,text width=1.75cm,text centered,minimum height=0.8cm]
rect_for_op=[draw,rectangle,fill=blue!10,text width=1.75cm,text centered,minimum height=0.8cm]
rect_for_decision=[draw,rectangle,fill=yellow!10,text width=1.75cm,text centered,minimum height=1cm]
diam=[draw,diamond,fill=green!5,text width=1.25cm,text badly centered,inner sep=0pt,minimum height = 0.5cm]
ellipse=[draw,Ellipse,fill=blue!30,text width=3cm,text centered,minimum height=1cm]


line=[draw,-latex']
A Periodogram-Based Method to Identify Forced and Natural Oscillations Using PMUs
This work is supported by the Fonds de Recherche du Québec - Nature et technologies under Grant FRQ-NT NC-253053 and by Mitacs Globalink Research Internship program.
    Qinye Tang

Zhejiang University

Hangzhou, China

Email:qtang@zju.edu.cn 
Xiaozhe Wang
McGill University

Montreal, QC Canada

Email: xiaozhe.wang2@mcgill.ca

    December 30, 2023
=======================================================================================================================================================================================================================================================





§ INTRODUCTION

One of the major threats to the security and stability of power systems is sustained oscillations <cit.><cit.>, 
which can damage equipment, degrade the power quality and increase the risks of cascading power failure. Oscillations in power systems can be categorized into natural oscillations and forced oscillations. Natural oscillations often result from inherent interaction among dynamic devices, which typically can be mitigated by power system stabilizer (PSS), intertie line controls, etc.  In contrast, forced oscillations refer to system responses to an external periodic perturbation <cit.>, which cannot be depressed by PSS. The only effective measure is to locate and remove the external driving source.


In previous literature, various methods have been proposed to identify natural oscillations including prony analysis <cit.>, Frequency Domain Decomposition Analysis <cit.>,  Subspace Identification <cit.>, Robust RLS method <cit.>, etc. However, these methods may not effectively identify forced oscillations. To distinguish different oscillation mechanisms based on PMU data, Liu et al. in <cit.> used support vector machines (SVM) to analyze different features of the two oscillation mechanisms. However, this method's accuracy relies directly on the envelop size of oscillation. Xie et al. in <cit.> proposed a distinguishing method based on spectrum, yet its performance may degrade when the external force frequency is close to system's natural frequency. Ghorbaniparvar et al. in <cit.> proposed a residual spectral analysis method leveraging on different residual spectral properties of the two oscillation mechanisms. However, a group of forecasting models needs to be calculated, which may hamper its real-time implementation. Wang et al. in <cit.>, used a periodogram-based method to distinguish limit cycle oscillation from forced oscillation, but its performance is not guaranteed if the bandwidth of the power spectral density (PSD) is blurred by strong noise or weak external force.  

In this paper,  a novel periodogram-based method to differentiate forced oscillation from natural oscillation is proposed. It will be shown that, mathematically, the limit cycle and the forced oscillation can be described by the same governing stochastic differential equation yet with different parameter values, based on which a PSD-based classifier can be used to distinguish different mechanisms. The proposed method can quantitatively extract signatures of different mechanisms from PSD and provide accurate classification results even when the forced oscillation frequency is close to the natural oscillation frequency.




§ MATHEMATICAL MODELS FOR SUSTAINED OSCILLATIONS IN POWER SYSTEMS

The power system dynamic model can be described as:

    ẋ   =   f(x, y)
    0   =   g(x, y, u)

where x∈ℝ^n_x is the corresponding state variables (generator rotor angles, rotor speeds, etc), y∈ℝ^n_y is
algebraic variables (bus voltages, bus angles, etc), and u is the vector describing stochastic behaviors in real-world power systems. 

Our main interest lies in the stochastic perturbations like load variations and renewable generations, which can be modeled as a vectorized Ornstein-Uhlenbeck process u̇ = -Cu + σξ, which is stationary, Gaussian and Markovian <cit.>.


C is a diagonal matrix related to the reversing times of the process; ξ is a vector of independent standard Gaussian white noise; σ the intensity of noise.

Replacing y by x and u using implicit function theorem and linearizing (<ref>) around its steady state, we have <cit.>:

    q̇ = A q + σ Bξ

where q = [ δ x,δ u ] , A = [ f_x-f_y g_y^-1g_x    -f_y g_y^-1g_u;                 0                -C ],
B = [ 0,I_n_u ]^T.
In the rest of the paper, we focus on the stochastic dynamic model described above, where q is a vector Ornstein-Uhlenbeck
process.

Sustained oscillations of power systems can be caused by different mechanisms and thus are represented by different mathematical models. In this work, we are interested in identifying the forced oscillation caused by external force and the limit cycle governed by system inherent dynamics.




 §.§ Forced Oscillation



 The normal form of the forced oscillation can be described as:
 
    q̇=[    γ -ω_0;  ω_0    γ ]q+F[ cosΩ t; sinΩ t ]+σ[ ξ_1; ξ_2 ]

 where  q∈ℝ^2, γ<0, and ξ_1(t), ξ_2(t) are independent standard Gaussian white noises. F denotes the external driving force acted to the system with a frequency Ω.

We want to find the closed-form solution of q(t) (particularly the first component), from which a nice structure of forced oscillation can be derived. Denote
M≜[    γ -ω_0;  ω_0    γ ],
then the solution q(t)= q_c(t) + q_ou(t)=[ x_c(t)+x_ou(t),y_c(t)+y_ou(t) ]^T can be found from the control process: 
q̇_̇ċ=Mq_c+F[ cosΩ t,sinΩ t ]^T, 

and the Ornstein-Uhlenbeck process:

q̇_ou=Mq_ou+σ[ ξ_1,ξ_2 ]^T. 







The control process 
can be easily solved by evaluating the matrix exponential:

    x_c=[ 1 0 ]∫_0^t e^(t-s)MF[ cosΩ s; sinΩ s ]ds=Fcos(Ω t)/√(γ^2+(Ω-ω_0)^2)

Similarly, x_ou can be solved by evaluating the Itô integral:

 



    x_ou=[ 1 0 ]∫_0^t e^(t-s)Mσ[ dB_s^(1); dB_s^(2) ]ds=σ∫_0^te^-γ(t-s)dB_s


where the Brownian terms dB_s^(1) and dB_s^(2) are independent. 
 
 
 
Therefore, the first component of q is of the form:

x(t)=x_c(t)+x_ou(t). 


By taking time derivatives on both sides and collecting terms, we obtain:

    ẋ(t)=-γ(x(t)-F/γ√(γ^2+Ω^2/γ^2+(Ω-ω_0)^2)cos(Ω t))+σξ(t)

where ξ(t) is standard Gaussian white noise (see detailed derivation in Appendix <ref>). 

In other words,  (<ref>) can be represented as:

    ẋ(t)=-β(x(t)-Q(t))
    +σξ(t)

where β=γ and Q(t)
=F/γ√(γ^2+Ω^2/γ^2+(Ω-ω_0)^2)cos(Ω t). 

In the following section, we'll show that incidentally the limit cycle can also be represented in the form of (<ref>), yet with the parameters λ, α, σ_1 locating in different ranges. Leveraging on this difference, we can develop a data-driven diagnostic to tell the two oscillations apart from each other.  



 §.§ Limit Cycle

 Stable limit cycle emerging from supercritical Hopf bifurcation is another common power system oscillation mechanism. It's usually regarded as an early warning sign of voltage collapse, since the Hopf bifurcation usually precedes the saddle-node bifurcation which leads to final voltage collapse <cit.>.
 
 
The normal form of the oscillation stochastic system near Hopf bifurcation can be described as:
 
    q̇ = [    γ -ω_h;  ω_h    γ ]q-|q|^2q+σ[ ξ_1; ξ_2 ]

where  q∈ℝ^2, γ>0, ξ_1(t), ξ_2(t) are independent standard Gaussian white noises.

As discussed in <cit.>, near the Hopf bifurcating point, the amplitude of the limit cycle grows with √(γ), and the angular frequency is approximately ω_h. The solution can be approximated as:

    q=([ x(t); y(t) ])≈([ (√(γ)+p(t))cos(ϕ(t)); (√(γ)+p(t))sin(ϕ(t)) ])≜([   N(t); N^'(t) ])

where ϕ̇=ω_h +σ/√(γ)ξ_ϕ and ṗ=-2γ p+σξ_p is an Ornstein-Uhlenbeck process independent of ϕ, which is a Brownian motion with deterministic drift <cit.>. 




We claim that x(t) can be described by the following dynamic equation:

    ẋ(t)=-β(x(t)-N(t))

where 
β→∞. To see this, taking Laplace transform on both sides of (<ref>), we have
x̂(s)=β/β+sN̂(s)+x(0)/β+s.
Hence by the uniqueness of inverse Laplace transform we obtain that

lim_β→∞x(t) = lim_β→∞ℒ^-1{x̂(s)}(t) = N(t).











 §.§ The Same Governing Equation Yet Different Parameters


Comparing (<ref>) and (<ref>), it is easy to see that they can be described by the same stochastic differential equation:

    ẋ(t)=-β(x(t)-Q(t;λ,α,σ_1))
    +σ_2ξ(t)

where 

    Q(t;λ,ω,σ_1)≈(λ+σ_1∫_0^te^-2λ^2(t-s)dW_s)cos(ω t+σ_1/λB_t)



However, the two oscillation mechanisms have different parameter values as summarized in Table <ref>. We intend to utilize this structure that is able to simultaneously describe the two oscillation mechanisms to develop a data-driven method to distinguish the two oscillation mechanisms.






§ PSD-BASED TECHNIQUE TO DISTINGUISH TWO OSCILLATION MECHANISMS


Inspired by the work <cit.>, we consider the power spectral density of the dynamical system (<ref>), from which a classifier can be readily proposed to distinguish the oscillation mechanisms. In addition, an optimization technique is exploited to to estimate parameters presented in Table <ref>.




 §.§ The Power Spectral Density


The power spectral density (PSD) of the linear dynamic system (<ref>) takes the following form <cit.>:



    PS[X](ω)=β^2/β^2+ω^2PS[Q](ω)+σ_2^2/β^2+ω^2

where the PSD of process Q(t;λ,α,σ_1) can be approximated by: 

    PS[Q](ω) ≈
    λ^2 F(
    σ_1^2/λ^2,α,ω)+
    σ_1^2/4λ^2F(
    σ_1^2/λ^2+4λ^2,α,ω)

The detailed expression of F is:

    F(s,α,ω)=2s4(α^2+ω^2)+s^2/[4(α-ω)^2+s^2][4(α+ω)^2+s^2]






 §.§ A PSD-Based Classifier

The general perception about the distinction between the limit cycle and the forced oscillation lies in the temporal decoherence <cit.><cit.>. The forced oscillation is characterized by a long-term temporal coherence because of the periodic driving, in contrast to intrinsic limit cycle, which has increasing spectral bandwidth due to temporal decoherence.

Particularly, the rightmost term in (<ref>) describes the effect of the non-decohering perturbation on the periodic solution Q(t;λ,α,σ_1). Ideally, the rightmost term in (<ref>) should be zero for limit cycle (σ_2=0) and far from zero for forced oscillation (σ_2=σ>0). In light of this, we propose the following classier: 


    V_PS:=
    σ_2^2/β^2/PS[x](0)
    =
    σ_2^2/β^2/PS[Q(t;λ,α,σ_1)](0)+σ_2^2/β^2

which relates the power due to non-decohering perturbations to the total power of DC term (i.e, the process power at zero frequency). 




It is straightforward from its form that  V_PS∈[0,1]. If the periodic time series data originates from a forced oscillation, σ_2>0 while PS[Q(t;λ,α,0)](0)=0, indicating that V_PS should be close to 1. In contrast, if the cyclic data results from a limit cycle, σ_2=0 and β=∞ ideally, thus V_PS should be very close to 0. 

It is worth mentioning that the value of classifier V_PS will not be dramatically influenced if the external force frequency Ω is close to system's mode frequency ω_0, since the only affected parameter is λ from Table <ref>, which nevertheless will not greatly affect V_PS. Such result is also validated in the numerical study in Section <ref>.





 §.§ A Data-Driven Method to Distinguish Different Oscillation Mechanisms



We have shown that the dynamics of two oscillation mechanisms, namely, the forced oscillation and the limit cycle, can be described with the same governing equation (<ref>) with different parameter values. In addition, there is a classifier V_ps that possesses a higher value (≈ 1) in the case of forced oscillation and a lower value (≈ 0) for limit cycle. 

We intend to leverage the above analytical results to develop a diagnostic method to distinguish the two oscillation mechanisms.  
Given PMU data with oscillations, we will first solve two constrained optimization problems to fit the parameters to each oscillation mechanism. The one with less loss function value will be selected for calculating the classifier V_PS, based on which the oscillation mechanism can be determined.   



More specifically, given a  pair of oscillator's time series {x_n}_n=1^N, {y_n}_n=1^N and sampling period Δ t, 
we estimate the sample PSD by:

    PS[x]_estimated(ω) = Δ t/N|∑_n=1^Nx_ne^-iω nΔ t|^2






Next, we estimate two sets of the five parameters listed in Table <ref> by fitting the estimated sampe PSD. Particularly, two constrained optimization problems are solved. First,  a limit cycle is assumed so that we estimate the parameters by solving:







    min_λ,α,σ_1,β,σ_2∫_0^ω_c[PS[x]_estimated(ω)   -PS[x(t;λ,α,σ_1,β,σ_2)](ω)]^2dω
    
        s.t.    σ_2   <ϵ,   β>M



where ω_c is the upper bound of the fitting interval; 
ϵ and 
M are two bounds set to refine the searching space of the parameters according to Table <ref>. 

Second, a forced oscillation is assumed so that we estimate the parameters by solving the second constrained optimization:



    min_λ,α,σ_1,β,σ_2∫_0^ω_c[PS[x]_estimated(ω)   -PS[x(t;λ,α,σ_1,β,σ_2)](ω)]^2dω
    
        s.t.    σ_1<ϵ,      σ_2/β=PSD[x]_estimated(0)

where 
ϵ is an upper bound for σ_1 based on Table <ref> and the second constraint is applied to guarantee the fitting performance of the DC term. 

It is worth mentioning that it is more efficient to solve two constrained optimization problems for each oscillation mechanism compared to solving one global optimization, since the dimension of the searching space is greatly reduced. 


Once we get two sets of the estimated parameters from the above constrained optimization problems,  
we choose the set of parameters which has less objective function value, i.e., the one that fits the empirical PSD better, to calculate the statistic V_PS. 
We classify the oscillation to be a forced oscillation if V_PS>0.5 or to be a limit cycle if V_PS≤ 0.5.

Note that we do not make the classification decision directly based on the optimization results to avoid potential incorrect results due to overfitting. V_PS on the other hand can effectively identify the different signatures of the two oscillation mechanisms around zero frequency and around the peak. 
Thus, combining the analysis on the whole interval [0 ω_c]  with a particular focus around zero frequency and peak frequency may provide more reliable and accurate identification results.  



To summarize, assuming that other possibilities like quasi-periodic chaos and weakly-damped oscillation have been ruled out (e.g., via kurtosis <cit.>), a flowchart describing the detailed steps of the proposed method to distinguish the limit cycle and the forced oscillation is presented in Fig. <ref>.








§ SIMULATION STUDIES
 
Fig. <ref> shows two different oscillation scenarios in IEEE 14-bus systems, from which the exact mechanisms can  
be hardly distinguished.
In this section, we test the performance of the  proposed method in identifying the oscillation mechanism of each scenario.  
The emulated PMU data is obtained from the time-domain simulation of the 14-bus system. Stochastic load fluctuations are considered, which follow Orstein-Uhlenbeck process with σ=0.01 in (<ref>). Besides, we use exponential recovery loads to model load dynamics, the model of which can be described as:

    P   =   kP_0(V/V_0)^α
    
    Q   =   kQ_0(V/V_0)^β

where k is a dimensionless demand variable, V_0 is the reference voltage, and α and β are active and reactive power exponents depending on the type of load. All simulations were conducted in PSAT-2.1.8<cit.>. The parameter values of the test systems are available on: http://github.com/xiaozhew/Test-Systems-. 



The estimated PSD using the PMU data between 1000 to 1500s are shown in Fig. <ref>, both of which have relatively blurred spikes, indicating that it is impossible to distinguish the oscillation mechanisms by eyeballing the width of the spike as suggested by <cit.>.

 





 §.§ Case A


The time series data shown in Fig. <ref> indicates that the systems evolves from the old steady state (0-300s) to the new steady state (700s-1500s) probably after some change of parameters. The estimated PSD of the new steady state clearly shows that there is an oscillation around 0.12 Hz, yet the PSD has certain bandwidth. We intend to apply the proposed method to figure out the exact mechanism of the oscillation.

Following the procedure proposed, we get two sets of parameter estimates by solving two constrained optimization problems under different assumptions. The fitted PSDs are presented in Fig. <ref>. The parameter values and the loss functions are listed in Table. <ref>. It can be seen the parameter set under the assumption of limit cycle has less loss function value, which is therefore used to compute the classifier V_PS. Since V_PS=9.82e-16, much less than 0.5, the oscillation mechanism is classified to be limit cycle.  







The physical cause of the oscillation has been detailed in <cit.>. Starting from 300s, the exponential recovery loads at Bus 5 and Bus 12 increase gradually. By
420s, both loads grow by 12% and stop increasing afterwards. 
The competing effect between load
dynamics and the over excitation limiters of the generators leads to the voltage instability. The power system passes the Hopf bifurcation point around
400s when a stable limit cycle is born.



 §.§ Case B

In case B, the estimated PSD of the new steady state shown in Fig. <ref> does not have a perfectly thin spike, making the identification based on the bandwidth around the peak <cit.> fails. Applying the proposed method, we fit the estimated PSD and obtain two sets of parameters by solving two optimization problems under different assumptions. The results are shown in Fig. <ref> and Table. <ref>. 

The loss function of the forced oscillation constraint is smaller,  leading to a classifier value V_PS=1. Therefore, the oscillation mechanism is identified to be forced oscillation. 





The actual situation is that the exponential recovery
load at Bus 5 increases by 5% between 300s and 350s, at which one cyclic load joins at Bus 5 with a forced frequency 0.15Hz. Therefore, the oscillation is indeed a forced oscillation. More importantly, although the forced oscillation frequency is very close to the natural frequency mode 0.12Hz, the proposed method is still able to provide accurate identification result. 



§ CONCLUSIONS AND PERSPECTIVES

This paper has proposed a periodogram-based method to distinguish natural oscillation and forced oscillation. It has been shown that the limit cycle and the forced oscillation can be described by the same stochastic differential equation yet with different parameter values, based on which a PSD-based classifier is proposed to distinguish different mechanisms. It has been shown via numerical simulation that the proposed method can provide accurate classification results even though the forced oscillation frequency is close to the natural one. Further investigation about locating the source of forced oscillation is needed.




§ APPENDIX



 §.§ Derivation of Equation <ref>


    x(t)=F/√(γ^2+(Ω-ω_0)^2)cos(Ω t)+σ∫_0^te^-γ(t-s)dB_s

Taking time derivative on both sides we obtain:

    ẋ=-ΩF/√(γ^2+(Ω-ω_0)^2)sin(Ω t)-γσ∫_0^te^-γ(t-s)dB_s+σξ(t)

Let ẋ=-γ(x(t)-Q(t))+σξ(t) and solve for Q(t): 

    Q(t)   =ẋ(t)+γ x(t)-σξ(t)/γ
       =F/√(γ^2+(Ω-ω_0)^2)cos(Ω t)-Ω/γF/√(γ^2+(Ω-ω_0)^2)sin(Ω t)
       =F/γ√(γ^2+Ω^2/γ^2+(Ω-ω_0)^2)sin(arctanγ/Ω - Ω t)

Note that the phase shift arctanγ/Ω depends on the initial condition. Without loss of generality, Q(t) takes the form:

    Q(t)=F/γ√(γ^2+Ω^2/γ^2+(Ω-ω_0)^2)cos( Ω t)



ieeetran

