Performance Analysis of Data-Driven Techniques for Detection and Identification of Low Frequency Oscillations in Multimachine Power System

In power systems, identification and damping of low-frequency oscillations(LFO) is very crucial to maintain the small signal stability. Hence the computation of eigenvalues, eigenmode shapes, participation factors, and coherency of generators are essential parameters of critical LFO modes. The existing data-driven approaches explore either one or two aspects of modal parameters from the dynamic pattern of the measurement data. In the present work, two approaches i) Iterative Approach(IA) ii) Non-Iterative Subspace(SS) method of data-driven techniques are used to estimate the state-space model of the system under study from the measurement data in a holistic framework. Based on the estimated system model, the eigenvalues of LFO, eigenmode shapes, participation factor, and coherency of associated generators participating in electromechanical oscillations are computed. Finally, from the estimated participation factors for the Inter-area oscillation mode (IAM), the Static Synchronous Compensator (STATCOM) damping controller is designed and placed at the generator with the highest participation factor for damping of inter-area oscillation. The enhancement of damping ratio of inter area mode with STATCOM damping controller is estimated and verified using IA & SS data driven approaches for the first time. In this work, IA uses prediction-error minimization algorithm (PEM) & Parallel computing techniques and SS method uses Multivariable Output Error State Space (MOESP) algorithm for the estimation of Hankel matrix from the measured data. The effectiveness of data-driven approaches are demonstrated by the simulation of a IEEE 4-machine,10-bus power system using MATLAB/Simulink. IA & SS methods incorporating wavelet based denoising techniques are very effective in identifying the LFO modes even with noisy measurement. The efficacy of the denoising to suppress the effect of noise is demonstrated by comparing with noiseless environment. The results of data-driven approaches indicate their high degree of accuracy and efficiency in being consistent with Eigenvalue analysis (EA) performed on the system.


I. INTRODUCTION
The increased power transfer over the transmission network infuriates the power system to undergo low-frequency oscillations(LFO). The LFO typically lies between 0.2 to 3 Hz. Generally, these oscillations decay fast and the system remains stable if the system has positive damping at these oscillating frequencies. Depending upon certain operating conditions the LFO may grow/ sustain to reach the extent of causing synchronous machines to lose synchronism. The adverse effect of these oscillations intensifies the fluctuations The associate editor coordinating the review of this manuscript and approving it for publication was Flavia Grassi . in voltage, power flow, torque, and speed for a long duration. The damping of LFO is a function of system operating conditions [1], [2]. Detection of electromechanical oscillations is of utmost importance in the power system to mitigate the blackout scenario. Analysis of measurement data is very much essential in control rooms to have real-time awareness of such oscillations to initiate emergency control in view of maintaining the stability and to ensure secured power system operation.
Numerous approaches have been introduced for extracting details of low-frequency oscillations from the dynamic patterns of measurement data [3], [6]. The Prony method is a parameter fitting approach, that is used to identify dominant modes from measurement data and was initiated by Haur [7]. The extension of the Prony algorithm by Zhou et al. [8] assuming that the homogeneous responses of the linearized power system model is a linear combination of decaying complex exponential functions and is suited if generators swing in unison. Performance of the Prony method deteriorates in the presence of noises in the measured data. Fladrin et al. [9] introduced estimation algorithm empirical mode decomposition (EMD) which uses pre-processed data using a multi-bank filter to remove the noise and signal offsets. However, due to the frequency mixing phenomenon in multi-bank filter, EMD introduces artificial modes apart from the dominant mode of the system. Eigenvalue realization algorithm introduced by Kamwa et al. [10] suppresses the artificial modes using multi-band modal analysis. The use of the sliding window approach adapted by methods reported is computationally exhaustive and the memory requirement for the storage of data is extensively large. Kalman filter and extended Kalman filter-based approaches were introduced in [11] & [12] uses the dominant modes from preceding estimations and updates the estimation from new samples.
Further, Zhou et al. [13] introduce an autoregressive moving average(ARMA) algorithm which estimates dominant modes from the ambient responses of the power system. Identification of mode shapes from the ambient data is proposed by Trudnowski et al. [14] and Dosiek et al. [15] introduce the ARMAX to identify the eigenvalues that correspond to electromechanical oscillations and mode shapes from the multichannel measurement. However, the performance of ARMAX is reduced for the measured signal containing high damping oscillatory modes. The stochastic subspace identification Ni et al. [16] and recursive adaptive stochastic subspace identification Nezam et al. [17] were proposed to overcome the limitations of ARMAX and increase the computational speed of recursive methods during estimation of mode shapes and dominant mode.
The participation factors of the generators for dominant modes are computed using EA [1], [2], geometric measures, and modal participation ratio methods proposed in [18], [19] & [20], all these methods require the computation of the state matrix from the linearized power system model at a given operating point. Due to changes in operating points, the estimation of participation factors of the dynamic model of the power system in real-time is seldom guaranteed.
Senroy et al. [21] proposed the Hibert-Huang Transform to identify the coherent group of generators from the measured data for the inter-area mode of oscillations. Avdaković [22] also proposed the wavelet-based approach to extract the mode shapes among the generators associated with IAM. To identify the coherency of generators for multiple dominant modes the cluster-based data-driven approaches were proposed such as hierarchical clustering [23], principal component analysis [24], independent component analysis (ICA) [25], Dynamic mode decomposition(DMD) [26], Koopman mode (KMs) [27] and spectral clustering [28]. One of the prerequisites of cluster-based methods is that they necessitate a presumptive knowledge of many coherent groups in the system which is not available always. The comprehensive approach which provides information on all the modal parameters such as eigenvalues, mode shapes, participation factors, and generator coherency is proposed by Li et al. [30] in which the Eigenvalue Realization Approach (ERA) is used to extract the information and requires a large set of data to construct a block Hankel matrix. The validation of the ERA method for the system incorporated with the damping controller is not reported. The performance of ERA under noisy signal environment is not been investigated.
Xiong et al. [31] explains the advantages, limitations and applications of different approaches for system modelling and stability analysis.The extensive survey is carried out on small signal stability and large signal stability analysis for the power system dominated with VSC. The technical challenges and limitation of the small signal stability analysis using state-space modelling of such systems are addressed. Comprehensive study indicate the need of very effective modelorder reduction considering timescale features, the stability classification and the aggregation of VSCs with similar characteristics. In Impedance based approach the computation process is simplified compared to state space approach but it is prone to measurement noise. This demands very effective signal de-noising approach.
In [32], Li et al. proposed robust and effective unified controller for each DG's connected in multi bus microgrid which controls the power flow effectively during grid connected, islanded and operational transition. The robustness of the proposed algorithm is analyzed comparing the performance for the laboratory setup using Realtime digital simulator with experimental setup using MATLAB simulation are verified.
Li et al. in [33] introduces a discretized time delay model for power system. In this work, only retarded variable are considered for discretization to reduce the model order and introduces sparse techniques to reduce the CPU time required for the eigenvalue computation for time delay power system. However, the robustness of the algorithm in noisy measurement is not investigated in the paper.
In this paper, EA is performed on a linearized dynamic model of the power system and is considered as a base case for the proposed approaches. The data driven techniques IA and SS are utilized to estimate the eigenvalues, mode shapes, participation factors, and generator coherency for multiple electromechanical modes from the measurement data using holistic framework. The performance of IA & SS approaches are evaluated again with denoised signals. Data required to validate the performance of data driven techniques is obtained from MATLAB/SIMULINK model of 4 machine, 10 bus system.
In this paper, the data-driven techniques IA and SS methods are used for accurate estimation of the Inter-area mode and Swing modes as compared to the model-based EA approach. PEM algorithm with Levenberg-Marquardt (LM) search algorithm [38] is used to get the maximum fitness function for the measurement data in IA and Multivariable Output Error State Space based algorithm is used to compute the weighting function for singular value decomposition in SS method. Parallel computing with global optimization methods are applied in the present work to speed up the iterative algorithm and the performance of both approaches are compared. The performance of proposed data-driven methods in estimating the modal parameters and its accuracy is found to be consistent with the EA approach in extracting dynamic characteristics such as participation factor of the generator, mode shapes, and coherency of the generator for particular eigenmode. It is observed that, IA & SS approaches meticulously acquire dominant modes, participation factors, mode shapes, and generator coherency. The computation speed of the IA approach is increased by incorporating parallel computing using a global optimization search algorithm. Based on the results, it is shown that SS method is faster than IA. During noisy environment the data driven methods IA & SS are not very effective to estimate the state space due to convergence issues in IA and larger orders for estimation in SS. The performance of proposed approaches during noise environment is tested & necessity of noise detection and denoising the signal prior to the estimation is justified. To overcome the difficulty due to the effect of noise, the statistical based noise detection mechanism and wavelet based denoising techniques are incorporated in this paper. The wavelet denoising with threshold techniques are chosen and the efforts are made to select the efficient wavelet with best decomposition level using statistical analysis.
STATCOM supplementary modulation controller is designed to damp out the IAM with highest possible damping factors by placing STATCOM at the generator which has the highest participation factor for IAM. The IA & SS approaches are used to detect the eigenvalues of IAM to demonstrate the effectiveness of the STATCOM damping controller for the first time.
The work reported in this paper is mainly focused on 1. Implementation of variant of N4SID is carried out for subspace identification and estimation of eigenmodes describing the LFO from the time series data of Wide Area Monitoring System. 2. Detection of presence of noise in the PMU measurement is using Ljung-Box-Q-test & choice of appropriate wavelet by performing adequate statistical analysis for denoising the noisy signal & using box plot. 3. Implementation of IA & SS methods are used to identify the eigenvalues of Inter area mode (IAM) of the system with STATCOM damping controller incorporated for multimachine power system. The organization of the paper as follows. EA for a linearized power system model is reviewed in section II. The state-space identification using a data-driven techniques iterative approach IA and noniterative approach SS are used to extract the dominant modes, mode shapes, participation factors and coherency between the generators from the measurement data is developed in section III. The details of noise detection & wavelet based denoising techniques are discussed in section III.D, III.E & III.F. Section IV presents the results & performance of proposed approaches with model-based eigenvalue analysis as a base case for the simulated system of 4-machine,10 bus system in Section IV.A & Section IV.B. Section IV.C describes the usage of parallel computing techniques to speed up the iterative approaches. Based on the participation factors STATCOM parameters are tuned to damp the IAM in Section IV.D. In Section IV.E & Section IV.F, the performance of proposed approaches are verified during noisy environment & noise detection and wavelet denoising techniques are discussed. The conclusions are drawn in Section V.

II. EIGENVALUE ESTIMATION USING MODAL ANALYSIS
Stability analysis in a power system is performed using a set of differential-algebraic equations (DAE) described aṡ where x is a state vector x ∈ R m , input vector i ∈ R k and the output vector y ∈ R n . Linearization of equation (1) at an equilibrium point x 0 , i 0 results in the state space representation of system expressed as where A is the state matrix and the dimension is m × m, B is control matrix with dimension m×k,output matrix C has n × m dimension and the dimension of feedforward matrix D is n × k.
To implement eigenvalue analysis (EA) method on state matrix A, the eigenvectors and eigenvalues must satisfy following equation where λ j is the j th eigenvalue of A, U j is right eigenvector of A associated with λ j and dimension is m × 1, left eigenvector of A is represented as V j associated with λ j .
To represent the relationship between the system states and the modes, the participation factor is computed by where P r,j is the participation factor which measures the participation of j th state variable over r th mode; v rj is the r th entry of V j and M EA jr is the r th entry of U j . The small-signal oscillation in the power system causes system oscillation to increase due to the presence of lowfrequency oscillation modes. Based on the frequency, the low frequency oscillatory modes are categorized as, inter-area oscillatory modes and swing modes.
The mode shape matrix is used to estimate the coherency between the generators associated with low-frequency oscillatory modes [29], [34] i.e represented as where M EA is mode shapes matrix associated with multiple low frequency oscillatory modes.

III. DATA-DRIVEN APPROACHES FOR MEASUREMENT DATA FOR STATE-SPACE ESTIMATION
The procedures of extracting eigenmodes, mode shapes, participation factors, and generator coherencies for the dynamic model of the power system is discussed in Section II. Accuracy of estimation of such methods is related to the accuracy of the models and parameters of the system which need non-trivial effort in the real-time operation of the power system. Therefore the measurement based approaches are predominantly used in monitoring real-time oscillations considering the impact of the real-time operation of the power grid on small-signal dynamic stability. Estimation of the eigenmodes is identified accurately without depending on specific dynamic models. Due to the advent of phasor measurement systems, the measurement-based methods are the need of the hour. Hence this section describes the data-driven approaches used to estimate the state-space model of the system through the measured data.

A. AN ITERATIVE APPROACH USING A PREDICTION-ERROR MINIMIZATION ALGORITHM
In an Iterative approach (IA), the discrete data collected from various measurement units such as PMU's and SCADA are used to obtain the continuous-time state-space model of the system. Initialization of parameter estimates is performed either using an iterative rational function estimation approach. The parameter values are refined further to obtain minimum estimation error using the prediction-error minimization (PEM) algorithm to improve the closeness of fit.
Consider model structure M x which is the set of parametrized set models. And M x (θ) is the particular model in M x associated with the parameter vector θ. The prediction error e (t, θ) for any model M x (θ) is given by where y(t) is actual output which is p -dimensional column vector,ŷ(t|θ) is a predicted output which is p -dimensional column & e(t) also be a p-dimensional column vector. For any dataset Z N the errors can be computed for t=1, 2, 3 . . . . . . .N. Thus for a parameter estimation, prediction error is computed (7) based on Z t . Select the θ N at t=N so that prediction error become as small as possible A weighted norm of prediction error is defined as The prediction error sequence e (t, θ) which is represented as e (t) and is defined as the difference between the predicted output and measured output of the model, & it is prefiltered using linear filter h −1 (q) the error is defined as where h is the prefilter used for preprocessing the data with is a scalar and e (t) is a vector. The subscript N denotes that the cost function is a function of the number of data samples and accuracy of cost functions increases for larger values of N.
In [37], the iterative approach uses a general linear estimator with Langrage's method to minimize the cost function. In this work, to minimize the cost function, the PEM utilizes Levenberg-Marquardt (LM) search algorithm [38]. PEM is used to select the proper order of system which generates different Hankel singular values for integer values of order, the order with lower Hankel singular values are discarded. The state-space model of the measured data is now extracted after observing the close fitness between the predicted data and measured input by choosing the proper order of the system. [34], [35].
The estimated continuous-time state-space model is represented asẋ where A 1 is a state matrix with m×m dimension, B 1 is control matrix with dimension m × k, C 1 is output matrix n × m, D 1 is feedforward matrix n × k. y(t) and i(t) output and input respectively. Further, by using A 1 , B 1 and C 1 controllable matrix B c 1 and observable matrix C o 1 of (10) are computed as where U IA is the right eigenvector represented as M IA k is the mode shapes of all the measurement channels associated with eigenmode λ k are computed using observable matrix C o 1 , (12) where U IA k is the right eigenvector corresponds to the λ k & M IA k is the kth column of C o 1 . Hence using mode shapes M IA k , the participation factor p IA kq of k th state variable in the q th mode; is computed as To identify the coherency between the generators, the angle of cosine d iq between the measure channels i & q can be determined using (14) with the use of estimated mode shape matrix M IA which is defined by The proposed approach uses numerical algorithms for subspace (SS) state-space system identification [39]- [41]. The additional parametrization is required for initial conditions while estimating the state-space using iterative approach from the measured data with the non-zero initial state, wherein performance of subspace approach do not vary for zero and non-zero initial states. Since it is a non-iterative method, the convergence issues, and parametric sensitivities to initial state estimation do not affect the performance of the algorithm. Hence the non-iterative subspace methods are faster than iterative approaches as described in section V.
The state-space model for the measurement data from the multichannel can be expressed as discrete form as where A 2 ∈ R m×m is a discrete-time state matrix, B 2 ∈ R m×n is discrete-time input matrix, C 2 ∈ R l×m are an output matrix, the feedforward matrix D 2 ∈ R l×n &y(kT ) and i(kT ) output and input vectors respectively. T is the sampling interval.
To estimate A 2 , B 2 , C 2 and D 2 matrices the subspace identification algorithm is employed in this work to predict the exact state-space model from the measured data. The procedures for estimation is summarized as follows (1) Construct input block Hankel matrix H ik , output block Hankel matrix and error block Hankel matrix from the measurement data which can be represented as where i k , y k and e k are the measurement data collected at time k; m 0 is the initial order of the system. Construct a matrix Z P which includes past input I P and past output Y P [40] defines as From Z P and future input I f , Y ef i.e. the predicted future matrix of the future output matrix Y f is computed using the least square algorithm. Inspired by the linearity of the system, we combine the past (I P ) and the future inputs (I f ) linearly to predict the future outputs (Y f ). and linear combinations are denoted using linear operators L Z , and L i defines in (19) Consider Y ef is the estimated output matrix The least-square estimate of Using orthogonal projection of Y f row space in the row space of I f and Z P the estimated output Y ef computed as The optimal values of L Z and L i are obtained using Oblique projection [37].
where Z P,0|i−1 is the first row to the ith row of Z P Taking singular value decomposition of L z Hence the order 'n' of the system is non-zero singular values of S 1 .
The extended controllability matrix where the subscript i denotes the number of block row and state space metrics arê Therefore the system model is computed as follows The discrete-time state matrices A 2 , B 2 , C 2 &D 2 are calculated using (26)(27), (30) and the state-space model can be represented as (15) with an error vector ρ z ρ i is neglected.
The continuous-time state-space model is obtained from the discrete-time model (15). The continuous time transformation of the matrices A 2 , B 2 and C 2 are where T is the sampling interval. Eigenvalues λ ss and eigenvectors U ss are calculated using If λ j is one of the modes of low-frequency oscillation mode of the power system under consideration and its frequency f j and damping ratio ζ j are computed as The controllable and observable matrices are computed as similar to (9) The participation factor p ss kq of the machines for participation factor p IA kq of k th state variable in the q th mode and M ss , the mode shapes matrix of the eigenvalues are calculated using (34) and (35) respectively In this work, variant of N4SID is implemented for subspace identification.N4SID algorithm [37], [42] is used for initialization & MOESP based algorithm [37] is used to compute the weighting function for singular value decomposition.
The subspace method used in this work make use of orthogonal projection to remove the input effect and the oblique projection to remove the error. Therefore, the proposed method is suitable for estimation of output from the recorded data such as time series data for the system with both forced excitation and unforced excitation.

C. COMPUTATION USING PARALLEL COMPUTING
To speed up the estimation process of proposed approaches the parallel computing using the global optimization tool of MATLAB is used. The various parallelized numerical algorithms such as Non-linear least squares and sequential quadratic programming algorithms are used for parallel computing without CUDA or MPI programming due to special array types and high-level constructs parallel for-loops. The toolbox allows using parallel computing -enabled functions in MATLAB and other toolboxes [43]. In the proposed work the computation speed of the iterative method is verified with non-linear least square and sequential quadratic programming algorithms [44] as compared with the estimation algorithm. However, it is found that, the proposed subspace algorithm is faster compare to the iterative algorithm and it do not require parallel computation.

D. MEASUREMENT NOISE AND SNR RATIO
In real-time measurement units in power system experiences noise which is characterized as white noise [45]. Uncorrelated sample values are one of the feature of white noise [46], the first i th sample the value have no correlation with (i+1) th samples. The required while Gaussian noise in dB is added to the signal generated through the time domain simulations to have similarity with practical signals.
The required noise is added to the signal and Signal to Noise Ratio is computed using (36) and (37) SNR rms = 20 log rms value of signal rms value of noise dB (36) Root mean square value of a signal 'x'is represented as where L is the number of samples of measured signal x(k) Since the rms value of the signal is the square root of its power i.e X rms = X power [47] and equation (36) can be written as, SNR POWER = 10 log signal power noise power dB (38) Noise power in dB = Signal power in dB -SNR in dB (39) Noise power = 10 (Noise power in dB) 10 Power of White Gaussian Noise is its variance itself., MATLAB function used for generating noise is given by Noise = standard deviation of white noise * randn (1, N) (41) VOLUME 9, 2021 The command randn(1,N) creates a row vector of random distribution with zero mean and unity standard deviation. The length of the Noise array N is same as length of signal array x(k). The standard deviation of noise is computed using equation (39). Therefore the Noise signal is obtained using (42) The present work is continued, considering the data matrix consisting of noisy signal measured at generator 1 and noiseless signals measured at remaining generator buses for Two-Area, 4 machine system. The iterative approach fails to obtain the state space estimation due to convergence issues in the noisy environment and subspace method demands higher fitness order n > 50 & time. This drives the necessity of identification of noise & appropriate denoising techniques to process the data. In this paper, the added white Gaussian noise with SNR of 50 dB is chosen, since it is the minimum signal strength possessed by any ringdown data in wide area monitoring unit [47]. Noisy signals results in misidentification of state space model using proposed approaches. Therefore, it is necessary to identify the noisy signal and followed by denoising the signal prior to system identification process. White noise is most common measurement noise in power system [48].

E. LJUNG-BOX-Q-TEST FOR IDENTIFICATION OF WHITE GAUSSIAN NOISE IN MEASURED SIGNAL
Ljung-Box-Q-test is a statistical test based on null hypothesis to support the existence of autocorrelation [49], [50], [52] and defined by where r i is the estimated autocorrelation of the data with 'n' umber of samples and k is the number of lags being tested.
where γ 2 1−a,h is the chi-square distribution table with 'h' degree of freedom and significance level α.
Hypothesis 'h' and 'pvalue' are estimated for the measured signal. The signal lags are identified as heavily autocorrelated if Q value is greater the critical value of Chi-square and hence rejects the null hypothesis. If, Q< critical value, then supports null hypothesis with poor autocorrelation which indicates the existence of white noise. The autocorrelationr i of the signal is computed for time series using (44), where N is the length of the time series, Y is the measured signal, E is the expectation value [49].

F. WAVELET DENOISING USING THRESHOLDING TECHNIQUES
The one dimension noisy signal is represented as where x (t) is the noisy signal, s (t) is measured signal through PMU and e(t) is added White Gaussian Noise subject to zero mean value & unity standard deviation. When a wavelet transform applied on signal x (t), it extract the signal energy into large approximation coefficients and distribute the noise energy over detailed coefficients. Therefore the denoising is viewed as exploitation over detailed co-efficient to obtained smooth signal using non-linear filtering process [52]- [56]. Thresholding is the commonly used denoising techniques proposed by Donoho and Johnstone [57] which removes the detailed coefficients less than certain threshold. The denoising process consists of three major steps as represented in block diagram in Figure 1. Step1: Selection of mother wavelet & decomposition levels and computation of wavelet coefficients.
Step2: Selection of threshold function and extraction of estimated wavelet co-efficient Step3: Signal reconstruction using Inverse wavelet transform.
Universal thresholding is the most widely used thresholding techniques because of effectiveness and simple approach. The formula for computation of threshold using universal threshold is given by where Nj is length of detailed coefficient and σ j is the mean variance at jth scale. The calculation of σ j is by computed by median estimation given by where W 1,l is the wavelet coefficient in scale 1.
In this work, universal threshold method is used for denoising the signal. There are two major nonlinear filters are present in the thresholding of wavelet. First one is, denoising using Hard thresholding shown in Figure 2a is defines aš  Second one is denoising using Soft thresholding shown in Figure 2b is defines aš In hard thresholding the coefficients between the threshold range −λ to λ are set to zero if the absolute values of coefficient lower than threshold values prior to reconstruction. Soft thresholding is also termed as wavelet shrinkage method which shrinks all positive & negative coefficients which lies in the threshold range [57] prior to reconstruction.
To evaluate the performance of various mother wavelets while denoising the signal with thresholding techniques, the signal to noise ratio (SNR) using equation (37), Mean Absolute Error (MAE), Mean Square Error (MSE), Correlations and Fitness coefficients methods are employed in this work. The computation of all the parameters are defined using following equations.
Mean squared error is defined by where n=number of data points, X i = Measured noise free signal alues D i = Wavelet denoised signal values Mean Absolute Error (MAE) is defined by where n=number of data points, X i = Measured noise free signal, D i = Predicted values after wavelet denoising, Correlation between the measured and wavelet denoised signal is defined by where X i & D i are the Measured noise free signal and denoised signals respectively,X &D are the mean values of the sample means of X and D. Fit co-efficient is employed to check whether the important information are lost during denoising, the fit co-efficient is given by (53) where X i & D i are the Measured noise free signal and denoised signals respectively.
The experiment analysis of the statistical features described equation 50-51 are discussed in the result section IV.A.

IV. RESULTS
The performance and application of the proposed method are verified with simulated data from 4 machine-10 bus system shown in Figure 3 [1], [2]. The measurement of the angle at generator buses is collected with Generator 1( G1) as a reference at a step size of 0.01 sec for 25 seconds. It provides 2500 frames in 25 seconds with sampling rate of 100 frames per seconds. In real-time systems the PMU's used in Wide area Monitoring System (WAMS )has output frequency varying between 10 to 60 frames per second and in advanced PMU [58], the reporting rate is up to 200 frames per seconds. In this work, the reporting rate for advanced PMU is considered as 100 frames /sec. The model-based approach is used as a baseline with a small signal disturbance condition when the mechanical input to Generators is increased by 10% at 0.5 secs. Figure 4 illustrates the variation in bus angles of generators G1, G2, G3, and G4 for the smallsignal disturbance. The data between 4s to 20s window is used as an input for the proposed approaches to estimate the eigenvalues.

A. EIGENVALUE ANALYSIS
The model-based approach is applied on the 4-machine, 10 bus systems considering 1.1 model of synchronous machine [1], [2]. The modal parameters of the machines such as eigenvalues, eigenvectors, mode shapes, and participation factors are computed using (3)(4) and (5) when a small disturbance of 10% variation in mechanical input is activated. The computed values of modal parameters are compared with the proposed approaches in the sections to follow.

B. AN ITERATIVE APPROACH AND SUBSPACE METHOD
Applying an iterative approach (IA) to identify the state-space of a system using bus angle data collected at generator buses at area1 and area 2 of the study system. The parameter estimation minimization algorithm (8) & (9) is applied to find the maximum approximation between the measured input and processed output. The Order of the system is identified by computing Mean squared error of fitness. MSE is the average squared difference between the estimated value & measured value of a signal and is defined by (50). In general, MSE represent the difference between the actual observations and observation values predicted by the model. In this work, the difference between the measured values and predicted values of a signal is calculated using MATLAB function immse() and lesser the MSE more the closeness of fitting. It is observed from Figure 5 that, order n=21 has minimum mean square error.   from the chosen data and the matrices A 1 , B 1 , C 1 and D 1 are computed as per equation (10) (11).
The algorithm is modified to extract the eigenvalues lying between 0.2 to 3Hz as the frequencies of our interest is to identify the low-frequency oscillation modes.
The eigenvalues λ n and the right eigenvector  (12). The participation factors of the machines are computed from the measured data using (13).
In applying Subspace (SS) approach to the simulated data of the multimachine system bus angles between 4 sec to 20 s window. Initially, the Hankel Matrix associated with   the measures input is computed using (16) and the order of the system n=21 is validated against closeness of fitness by solving Singular value decomposition S1(25) which was found to be containing all non-zero singular values. The percentage of the closeness of fitness between measured data and predicted data is shown in Figure 6.2 which is very similar to the IA method. Further, the discrete state-space matrices A 2 , B 2 , C 2 and D 2 are extracted using (30) and computation of continuous-time state-space matrices A t 2 , B t 2 , C t 2 and D t 2 are achieved using (31). Eigenvalues of A t 2 and right eigenvectors U ss are estimated. Controllability and observability matrices B tc 2 and C to 2 are estimated using the right eigenvector of A t 2 via (33). The participation factor of generators is computed using (34). The verification of estimation accuracy of proposed data driven approaches IA & SS for n=21 are evaluated, and compared with critical LFO modes of EA method as the basis for correctness estimation.
The dominant modes of the system are computed using the proposed data-driven approaches and EA method are compared as shown in Table 1. It is observed from Table 1 that, the eigenvalues estimated with IA & SS methods closely matches with EA method. Further, the damping ratio (ζ ) and frequency of oscillatory modes are calculated using (32). The %error of damping ratio (ζ ) and frequency of oscillatory VOLUME 9, 2021   modes(f) is found to be significantly small using both the IA and SS method of state-space identification. The %error of ζ and f with SS is slightly high than the IA approach for mode 2 shown in Table 2. However, SS method is faster than the IA approach is discussed in section IV.C.
The mode shapes of different dominant modes for the IA approach is shown in Figure 7. For Inter-area mode in Figure 7a, G1, and G2 form coherent group to swing against coherent group formed by G3 and G4, and Generator 1 in area 1 show more dominance for this eigenmode. For mode M1 in Figure 7b, the coherent group constituted by G1 & G3 swing against G2 & G4 and the generator G2 shows more dominance over M1. For mode M2 in Figure.7c, the coherent group constituted by G1 & G4 swing against generators G3 & G2 and G3 shows more dominance over the mode under consideration.
The mode shapes of different dominant modes for the SS approach is shown in Figure 8. For Inter-area mode shown in Figure 8a,, G1, and G2 form coherent group to swing against coherent group formed by G3 and G4, and Generator 1 in area 1 show more dominance for this eigenmode. For mode M1 in Figure.8b, the coherent group constituted by G1 & G3 swing against G2 & G4 and the generator G2 shows more dominance over M1. For mode M2 in Figure 8c, the coherent group constituted by G1 & G4 swing    in Figure 9. This validation ascertains that the proposed data-driven techniques accurately estimates the mode shapes from the measurement data. Participation factors of the VOLUME 9, 2021  generators for all the oscillatory modes are computed using (13) and plotted in Figure 10 for EA, IA & SS methods.
It is interesting to observe from the Figure 10 that, participation factor of Generator1 is maximum for IAM irrespective of EA, IA and SS method. Similarly for swing mode 1 (M1), the G2 contributes the highest participation factor and for swing mode 2 (M2), the G3 Contributes the highest participation factor. It is to be noted that, although the magnitudes of participation factors computed using data-driven approaches are different compared to EA method, the prediction of the generator with the highest participation factor for the particular mode is accurate and is in agreement with EA. It is not surprising to see that, magnitude of participation factor obtained by IA & SS methods are found to be higher than EA method. This is mainly because, EA approach is derived from the mathematical model using state variables wherein the data driven approaches performed using measured data. Hence, for the chosen order the proposed approaches do not replicate all the eigenvalues of the dynamical system. However, it is evident from the results that, the low frequency oscillation components of the system are in agreement with EA approach.

C. PARALLEL COMPUTING USING GLOBAL OPTIMIZATION FOR ITERATIVE APPROACH
Though the estimation accuracy using the Iterative approach is more accurate as seen in section 3 but the speed of computation is less as compared to the SS approach. In this work, the parallel computing tools are used to verify the increase in the computation speed of IA estimation.
The speed of computation using the Nonlinear Least square method is comparatively faster than sequential linear programming. However, the estimation using the subspace method is faster than the IA approach. The computation time required by various methods for the completion of the estimation process is shown in Table 3 using a 16GB RAM and 1.99GHz processor in MATLAB environment.

D. DAMPING OF IAM USING STATCOM
As per the results obtained in section [IV.B], it is observed that the bus near the G1 is the best place to locate the STATCOM Supplementary modulation controller to damp out the IAM. Figure 11 depicts the block diagram representation of supplementary modulation controller (SMC) with lead/lag network for STATCOM.
The control signal considered is the Thevenin voltage and is synthesized by making use of the STATCOM current and magnitude of the voltage at the STATCOM bus. X 1 is the tunable reactance and K r is the reactive current modulation controller gain, I r is the reactive current injected by the STATCOM into the bus. V l is the magnitude of load bus voltage where the controller is connected. Dynamics of the STATCOM is characterized by the first order plant transfer function with time constant T p (0.02sec).T C is time constant of derivative circuit s 1+sT c and it is considered as 10msec. T 1r & T 2r are the time constants corresponds to the compensator and m r denote the number of stages associated with compensator block of STATCOM SMC [59].
Controller Transfer function is given by K r is the reactive current modulation controller gain m r is the number of stages in compensator circuit and T 1r &T 2r are the time constants for compensator circuit and the compensator design is not considered in the present work, therefore the value of m r = 0 and equation (54) is written as The design of STATCOM SMC for the 4 machines, 10 bus power system is discussed in APPENDIX-I.
The tunable parameters for STATCOM SMC are X 1 which is the value of the tunable entity to synthesize the Thevenin's voltage and K r is the reactive current modulation controller gain.Using Sequential Linear Programming as per [59], [60], the value of the parameters are calculated with respect to the required decrement factor σ = −.25 & σ = −.5 are K r = 46.29, X 1 = 0.01029 and K r = 46.29, X 1 = 0.014 respectively to damp the IAM. Figure 12a & 12b shows the response of the rotor angles of all generators and the STATCOM reactive current for σ = −. 25. Figure 13a & 13b shows the response of the rotor angles of all generators and the STATCOM reactive current for σ = −.5. It is evident from the response that, the increasing value decrement factor by properly tuning the parameters of STATCOM results in fast damping of Oscillations. The eigenvalues with improved damping ratio correspond to IAM are computed for the variation in STATCOM tuning parameters using EA. The accuracy of IA and SS are compared with EA by estimating the eigenvalues for IAM when the STATCOM controller with tunable parameters for σ = −.25 & σ = −.5 is located at the bus near to the G1.
It is observed from the results shown in Table 4, both of the data-driven techniques IA & SS approximates the eigenmode computed using the EA approach. Time elapsed by IA for the estimation is 6.84sec with parallel computing and SS is 3.345sec. The estimation accuracy is computed as shown in Table 5, which shows that, both  IA & SS estimates the damping ratios and frequency of modes effectively and found to be consistent with values estimated using EA.

E. SIGNAL WITH MEASUREMENT NOISE
In this section, the signal measured at generator bus 1 i.e 'δ 1 ' is treated with White Gaussian noise at SNR of 50dB using (42). The Figure 14 and Figure 15 shows the noiseless and contaminated signal. The existence of autocorrelation is computed for the noiseless signal measured at Generator 1 with its time lagged version as in Figure 16. Wherein for noisy signal with its time lagged version there is no autocorrelation as shown in Figure 17. In this work, significance level α= 0.1 is kept as a reference & the autocorrelations are computed for 20 lags for the samples under test. The chi-square distribution table is constructed by computing Ljung-Box-Q-test for both noiseless and noisy signals with confidence interval level of 90%.The MATLAB command in (48) is used to perform the test.
Here, h indicate the null hypothesis, pV is the p-value, Q is the Ljung-Box parameter computed using (43) and cV is Chi-squared critical value [h pV Q cV]=lbqtest(x, 'Lags', number of lags, 'alpha', significance value). Table 6 shows the observations of lbqtest conducted for noise free signal. It is observed that, Q value is greater than c-value for all sample lags & rejects the null hypothesis. Hence the alternate hypothesis of existence of strong correlation hold good. The results of lbqtest conducted for noisy signal of SNR 50dB is shown in Table 7. It is observed that, the Q value is lesser than c-value for all the sample lags. Hence it support the null hypothesis. Alternative hypothesis fails, due to poor autocorrelation between the measured noisy signal and the time lags. Therefore, the lbqtest proves the presence of white noise in the signal. The mother wavelet is selected based on the effect of denoising by various wavelets such as Haar, Daubechies, Symlet, Coiflets, Biorthogonal and Reversebior. For denoising studies most commonly used mother wavelets are Symlet 8 (sym8), Daubechies 10 (db10) and Coiflets 5 (coif5) due to their effectiveness in denoising the signal mentioned in [51], [53], [54]. The wavelet thresholding techniques is considered as a suitable tool to process the signal with noise. In this work, the efficacy of various wavelet thresholding methods to suppress the White Gaussian Noise are analyzed. Performance of sym8, db10 & coif5 for contaminated signal 'δ 1 ' is carried out under soft and hard universal thresholding environment with default threshold factor λ = 0.05 at different decomposition level as shown in Figure 18a and Figure 18b.
The box plot constructed for denoised signals of various mother wavelets & compared with noisy and noiseless signals. It is observed from the Figure 18a, the box representing noisy signal(box1) has red colored '+'symbol which indicates the outlier in the signal. Among the remaining boxes, the db10-L8s, db10-L8h, db10-L9s, db10-L9h, sym8-L8s, symL8h, sym8-L9s, sym8-L9h and coif5-L7s, coif5-L8s, Further, the efforts are made to select the potential mother wavelet & best decomposition level based on the different statistical measures such as, SNR, mean square error (MSE), percentage fitness between denoised & measured signal, Mean absolute error (MAE) and correlation between denoised and measured noiseless signal. The statistical features for denoised signal using sym8, db10 and coif 5 with different level of decompositions are computed using universal thresholding for soft & hard threshold levels as shown in Table 8 using MATLAB. In comparison with hard thresholding, the statical features obtained through soft thresholding for all the wavelet at majority of the decomposition levels yields better results. Hence, in this paper the soft thresholding method is chosen for the further analysis. It is observed that, db10-level8 soft threshold, symlet8-level8 soft threshold and coif5-level 9 soft threshold have better statistical features compared to other wavelets. Figure 19 a depicts that, db10-level8 soft thresholding has better SNR in dB compared to its counter parts. As shown in Figure 19 b & 19 c, %fit, correlation and MAE of db10-level 8 is better than coif5-level-9, But MAE & MSE of sym8-level8 soft threshold is better than the db10-level8 soft threshold to a small extent which can be neglected. However, since the computed SNR of the denoised signal using db10-level8 soft threshold is approximately 1.65 times of sym8-level8 soft threshold & 1.218 times of coif5-level-9 with considerably better SNR, MSE, MAE & correlation. Therefore, db-10 is chosen as mother wavelet & level 8 is the best decomposition level because of maximal SNR, minimal MAE, minimal MSE maximum %fit and maximum correlation. The denoised signal constructed by db10-level 8 soft threshold is chosen for the further system identification process.
Data vector is constructed again using denoised signal and analysis of state space identification in continued for the system under consideration using proposed approaches. Due to convergence issues, IA is not preferred method during noisy environment. SS method gives maximum fitness at   n=30 using (50) for denoised data vector. Comparison of the damping ratios and oscillation frequency of the dominant eigenvalues obtained from proposed methods with EA is shown in Table 9.
It is observed that, the %error of damping and frequency of IAM is negligibly small as compared to EA. However, the swing mode 1 and swing mode 2 has significantly small %error both in frequency and damping. Participation factors for LFO's and mode shapes of the denoised signal is approximately similar to the noiseless environment as shown in Figure 20 and Figure 21. The value of participation factors are equaling with eigenvalues analysis because in realtime scenario the operating points are not fixed whereas EA analysis is carried out for fixed operating point. However the highly participating machines for respective eigenmodes obtained by proposed approach is identical to that EA approach. Hence, it is found that, with presence of noise SS with wavelet denoising is potentially suitable method to identify the dominant oscillatory modes of the system.

V. CONCLUSION
In this paper, state-space model of the IEEE-4 machine-10 bus system is estimated using data-driven approaches from the measurement data. The performance of i) Iterative approach and ii) Structured subspace is validated by comparing the data driven estimated eigenvalues with the eigenvalues computed using traditional state space approach using MATLAB/Simulink. Data-driven techniques are applied to estimate dynamic characteristics of the studied power system from the measurement data. The proposed work is also focused on noise identification based on statistical analysis such as lbqtest & boxplots for the measured signal. Lbqtest and boxplots collectively enhances the reliability of noise detection by avoiding the chances of false prediction. On detection of noise, the wavelet based thresholding techniques are applied to denoise the noisy measurement signal. The performance of both IA and SS are evaluated further for denoised signals and results obtained are found satisfactory and consistent with eigenvalue analysis method. The performance of the system for proposed approaches are demonstrated by the synthetic data collected from the simulation of 4 machine −10 Bus system.
The performance and novelty of the work is summarized as follows: PEM algorithm with LM optimization is used to get the maximum fitness function for the measurement data in IA and MOESP based algorithm is used to compute the weighting function for singular value decomposition in SS method. The IA and SS methods can accurately estimate the Inter-area mode and Swing modes at par with EA approach. IA and SS data-driven approaches are capable of capturing comprehensive electromechanical oscillation from the dynamic patterns of the power system and reveal the inherent features of the Power system. The performance of data-driven methods estimate the dynamic characteristics of system which is in consistent with eigenvalue approach. The performance challenges of data driven methods under noisy signal environment is successfully addressed. Measured signals are successfully tested for the presence of white gaussian noise using Ljung-Box-Qtest and box plot. The wavelet basis db10 with decomposition level −8 is highly effective compared to other wavelet basis in denoising the signal. SS data driven approach with wavelet denoising is suggested as a potential method for identifying the dominant oscillatory modes of the power system under dynamic conditions. Data driven techniques are also helpful in identifying the location of STATCOM damping controller based on the estimated participation factor. IA & SS methods are capable of detecting the eigenvalues of IAM of the system with STATCOM damping controller and is incorporated for the first time. The effectiveness of damping controllers can be enhanced by tuning the controller parameters with the help of data driven techniques.

Derivation of Control Law for STATCOM SMC:
The power flow through ith transmission line is given by Linearizing the equation, we get (57) where ∅ i denotes the phase difference across ith branch, V ji , V pi denote the voltage at the two ends of ith branch. The instantaneous power consumptions at all the branches in the power system are given by The reactive current injected at l th node is given by i Rl = −V l b sh l + n pl kl=1 V l x kl − V kl cos ∅ kl x kl (59) where b sh l is the shunt susceptance at the bus l, ji, pi are the two noes of branch i, n pl indicates the number of branches connected to the bus l. V kl is the voltage at the node of branch kl and ∅ kl phase angle across the branch kl.
The non-linear reactive current given in (4) (60) Multiplying V l on both sides of equation (5) It can be interpreted from(12) that the average loss causes by the term n l=1 ∂i Rl ∂V kl V l . Therefore the average electromechanical loss is given by For injection of reactive current at bus l Based on equation (20) the equivalent circuit is drawn as shown in Figure 22. Due to the presence of the Thevenin impedance X th 1 jw m the term the power dissipation in the circuit is limited. To prevail over this, the modification is incorporated in the control law such that, Overall block diagram of the STATCOM supplementary modulation controller is as shown in figure 11.