Estimating Participation Factors and Mode Shapes for Electromechanical Oscillations in Ambient Conditions

1 Estimating P articipation F actors and Mode Shapes for Electromechanical Oscillations in Ambient Conditions Xiaozhe W ang, Member , IEEE, Ilias Zenelis Abstract —In this paper , a new technique is applied to con- duct mode identification using ambient measurement data. The proposed hybrid measurement- and model-based method can accurately estimate the system state matrix in ambient conditions, the eigen values and eigen vectors of which r eadily provide all the modal knowledge including frequencies, damping ratios, mode shapes, and more importantly , participation factors. Numerical simulations sho w that the proposed technique is able to pro- vide accurate estimation of modal knowledge for all modes. In addition, the discrepancy between the participation factor and the mode shape is shown through a numerical example, demonstrating that using the mode shape may not effectiv ely pinpoint the best location for damping control. Theref ore, the proposed technique capable of estimating participation factors may greatly facilitate designing damping controls. Index T erms —mode identification, electr omechanical oscilla- tions, participation factors, small signal stability , phasor mea- surement units I. I NTRODUCTION Accurate estimation of low-frequenc y electromechanical oscillations is of paramount significance to ensure the se- cure operation of any power system. There are typically two approaches to conduct such estimation: linearizing the complete power system dynamic model around its steady- state operating point or utilizing measurement data to extract modal information (e.g., [1], [2]). Since an accurate model along with parameter v alues is usually una vailable in practice, much interests hav e been giv en to de velop measurement-based methods utilizing high-sampled and synchronized data from phasor measurement units (PMUs). PMU data can be further categorized into ring-down data and ambient data. Particularly , in this paper, the proposed technique utilizes ambient PMU data that can be easily obtained in normal operating condition without any large disturbance. In previous study , a large number of mode shape estimation methods using ambient data have been developed, including the spectral method [3], the frequency domain decomposition (FDD) method [4], the channel matching method [5], the transfer function (TF) method [6], and the subspace meth- ods [7], [8]. Particularly , it has been shown in [9] that the subspace methods may achieve more accurate estimations in high-damping situation and in the presence of oscillatory This work is supported by Natural Sciences and Engineering Research Council (NSERC) Discovery Grant, NSERC RGPIN-2016-04570. Xiaozhe W ang and Ilias Zenelis are with the Department of Electrical and Computer Engineering, McGill Uni versity , Montr ´ eal, QC H3A 0G4, Canada. email: xiaozhe.wang2@mcgill.ca, ilias.zenelis@mail.mcgill.ca disturbances compared with the other methods, and hence hav e receiv ed wide interests and applications. Howe v er , there are some limitations of the subspace methods. First, the estimated state matrix by subspace methods is not the true system state matrix; therefore, the left eigen vector associated with the mode of interest can not be estimated, leading to unattainability of the participation f actor [10]. Second, the subspace methods are computationally expensiv e because they require singular value decomposition (SVD) of large-dimensional matrices. In this paper, we will apply a novel method proposed in [11] [12] to carry out mode identifications for electromechanical oscillations, which may potentially overcome the challenges faced by the subspace methods. The proposed method is able to accurately estimate the true system state matrix and thus all the modal knowledge including frequencies, damping ratios, mode shapes, and more importantly , participation factors that carry crucial information for control. T o the knowledge of the authors, there have been no measurement-based methods on estimating participation factors without knowing the detailed network model. In addition, the proposed mode identification method is computationally efficient and can be implemented in near real-time. Simulation results demonstrate that the proposed mode iden- tification method can accurately estimate the modal knowledge (e.g., frequencies, damping ratios, mode shapes, participation factors) for all modes . Despite of the fact that the method is dev eloped based on the classical models, it sho ws good applicability to high-order generator models with detailed control. Furthermore, the proposed method is robust against measurement noise and thus is applicable in real-world ap- plications. The last but not the least, we show a potential discrepancy between the mode shape and the participation factor regarding the components with largest magnitudes, which implies that relying on the mode shape may not lead to the best locations for damping control. Therefore, the proposed technique capable of estimating participation factors may greatly facilitate damping control design by pinpointing the most problematic components, i.e., the best locations. The rest of the paper is organized as follo ws. Section II introduces the po wer system stochastic dynamic model and elaborates the proposed mode identification method. Section III presents the results of numerical studies. II. THE METHODOLOGY FOR MODE IDENTIF ICA TION In this paper, we are interested in the case when the power system operates in ambient condition and its dynamics can be 2 described by the classical generator model: ˙ δ = ω (1) M ˙ ω = P m − P e − D ω (2) with P e i = n  j =1 E i E j ( G ij cos( δ i − δ j )+ B ij sin( δ i − δ j )) (3) where δ , ω are generator rotor angles and relative rotor speeds w .r .t. the synchronous speed; M = diag ( M 1 , ...M n ) , D = diag ( D 1 , ...D n ) are inertia and damping constants, re- spectiv ely; P m and P e are generators’ mechanical po wer and electrical po wer; E = diag ([ E 1 , ..., E n ]) are the magnitudes of electromotiv e force behind the transient reactance; G ij + jB ij is the ( i, j ) th entry of the reduced admittance matrix that includes generators’ impedances. Similar to the approach in [13], we assume that the sys- tem loads modeled as constant impedances are experiencing Gaussian variations which can be translated to the diagonal elements of the reduced admittance matrix as below: Y ( i, i )= Y ii (1 + σ i dW i ) ∠ φ ii (4) where W is a standard W iener process and σ 2 i is the variance of load variations. Therefore, the system model takes the following form [11] [13]: ˙ δ = ω (5) M ˙ ω = P m − P e − D ω − E 2 G Σ ξ (6) where ξ = ˙ W is a vector of independent standard Gaussian random variables, and Σ= diag ([ σ 1 , ..., σ n ]) denotes the variation intensity . Note that the term E 2 G Σ ξ translates the stochastic v ariation at load side to the stochastic variation of electrical power at generator side [13]. Representing (5)-(6) in a compact form as a set of stochastic differential equations: ˙ x = A x + B ξ (7) where x =[ δ , ω ] T (8) A =  0 I n − M − 1 ∂ P e ∂ δ − M − 1 D  (9) B =[ 0 , − M − 1 E 2 G Σ] T (10) Specifically , x is a multiv ariate Ornstein-Uhlenbeck process [14]. For a generic Ornstein-Uhlenbeck process, the station- ary covariance matrix C xx =  C δδ C δω C ωδ C ωω  satisfies the following L yapunov equation if the system state matrix A is stable (typically satisfied in ambient conditions) [14] [15]: AC xx + C xx A T = − BB T (11) This relationship combines the statistical properties of states that can be estimated from PMU data and the physical model knowledge, providing an ingenious way to estimate model information from measurements. Substituting the detailed e x- pressions of A and B to (11), we obtain: ( ∂ P e ∂ δ )= MC ωω C − 1 δδ − DC ωδ C − 1 δδ (12) This deriv ed relation indicates that the dynamic state jacobian matrix ∂ P e ∂ δ can be estimated from the statistics of state variables extracted from the PMU measurements given that the damping factors D and inertias M of generators are known. Furthermore, we can readily construct the system state matrix A by (9) from which the electromechanical modes, mode shapes and participation factors can be obtained. A. The estimation of the covariance matrices The stationary covariance matrix is defined as: C δδ = ⎡ ⎢ ⎢ ⎢ ⎣ C δ 1 δ 1 C δ 1 δ 2 ... C δ 1 δ n C δ 2 δ 1 C δ 2 δ 2 ... C δ 2 δ n . . . . . . . . . . . . C δ n δ 1 C δ n δ 2 ... C δ n δ n ⎤ ⎥ ⎥ ⎥ ⎦ (13) where C δ i δ j =E [ ( δ i − μ i )( δ j − μ j )] , and μ i is the mean of δ i . H ow eve r, C δδ is typically unknown in practice due to limited data. Therefore, sample covariance matrix Q δδ is used to estimate C δδ , each entry of which is calculated as below [14]: Q δ i δ j = 1 N − 1 N  k =1 ( δ ki − ¯ δ i )( δ kj − ¯ δ j ) (14) where ¯ δ i denotes the sample mean of δ i , and N is the sample size. Likewise, C ωω , C ωδ can be estimated by Q ωω , Q ωδ . A window size around 500 s is used in the examples of this paper , which shows good accuracy . Like other measurement- based methods, there is always a tradeoff between estimation accuracy and sample size. As a result, the dynamic Jacobian matrix can be obtained from:  ∂ P e ∂ δ  = MQ ωω Q − 1 δδ − DQ ωδ Q − 1 δδ (15) if generator inertias M and damping constants D are known. B. The mode identification algorithm W e assume that PMUs are deployed at all the substations that generators are connected to. This assumption may be too optimistic currently , while is fairly reasonable in the near future considering the increasing rate of deployment of PMUs. Rotor angle δ and rotor speed ω in steady state can be calculated from PMU measurements. Discussion how exactly it is done is beyond the scope of this paper and the readers are directed to, e.g., [16]. If the inertias M and damping D of generators are kno wn, the following algorithm identifies and estimates electromechanical oscillations via ambient PMU measurements: Step 1. Compute δ and ω from the PMUs, and estimate the sample cov ariance matrices Q δδ , Q ωω and Q ωδ . Step 2. Estimate the Jacobian matrix ∂ P e ∂ δ by ( ∂ P e ∂ δ )= MQ ωω Q − 1 δδ − DQ ωδ Q − 1 δδ , and construct the system state matrix A =  0 I n − M − 1 ∂ P e ∂ δ − M − 1 D  . 3 Step 3. Compute the eigen v alues Λ , the right eigen v ectors Φ=[ φ 1 ,,,,. φ 2 n ] and the left eigenv ectors Ψ= [ ψ T 1 , ...ψ T 2 n ] T of A . The mode frequencies and the damping ratios can be directly obtained from Λ ; mode shapes correspond to Φ ; the participation factor of the mode λ i can be calculated by P i = [ P 1 i , ...P 2 n,i ] T =[ φ 1 i ψ i 1 , ..., φ 2 n,i ψ i, 2 n ] T . III. RESUL TS AND AN AL YSIS T wo cases in the IEEE 39-bus 10-generator system are to be presented. The first case is to v alidate the proposed mode identification method using the second-order generator models. The second case is to illustrate that the proposed method works well for high-order generator models with detailed control strategies, even though it is de veloped based on the classical generator models. In addition, it is shown that the method is robust against measurement noise, indicating good feasibility in practical applications. The parameters of the two cases can be found in: https://github .com/zenili/Mode-Participation-Estimation-2017. All simulations were done in PSA T -2.1.9 [17]. A. Case I: V alidation on Classical Gener ator Model In the 39-bus system, 10 generators are modelled as the second-order models as in (1)-(2); the angle of Generator 1 serves as the reference. W e set σ i =5 for all generators describing the stochastic load variations in the system. The sampling rate is 20 samples comparable with typical PMU sampling rate. The emulated trajectories of the rotor angle and the rotor speed of Generator 2 obtained from the PMU data are presented in Fig. 1, from which we see that system is stable and the state variables are fluctuating around their nominal operating states.                        Fig. 1: T rajectories of δ 2 and ω 2 in the 39-bus system using the classical models. Applying the proposed mode identification algorithm using the 450s sample data shown in Fig. 1, we obtain the estimation results of mode frequencies and damping ratios as presented in T able I-II. It is observed that the proposed method is able to correctly distinguish dif ferent modes and accurately estimate their frequencies and damping ratios, ev en though the frequencies of some modes are close to each other (e.g., Mode 2-4) and some damping ratios are high (e.g., Mode 1). T o further illustrate the accuracy of the proposed technique, we present a comparison between the estimated eigenv alues and the true eigen values of the system state matrix A in Fig. 2 (only the positive frequencies are shown), from which we Mode T rue Estimated Estimation frequency (Hz) frequency (Hz) error (%) 1 0.551 0.553 0.30 2 2.574 2.573 0.05 3 2.505 2.494 0.45 4 2.449 2.425 0.97 5 1.337 1.336 0.06 6 1.489 1.482 0.5 7 1.671 1.663 0.54 8 1.891 1.867 1.27 9 1.828 1.829 0.06 T ABLE I: Case I: the estimation results of frequencies for all modes. Mode T rue Estimated Estimation damping ratio(%) damping ratio (%) error (%) 1 6.12 6.03 1.40 2 0.20 0.19 4.55 3 1.14 1.15 0.89 4 2.44 2.48 1.37 5 2.07 2.16 4.39 6 2.34 2.36 0.96 7 1.09 1.03 5.56 8 2.58 2.62 1.77 9 2.99 2.97 0.49 T ABLE II: Case I: the estimation results of damping ratios for all modes.               Fig. 2: Case I: a comparison between the estimated eigen values and the true eigenv alues. see that all eigen values can be well estimated with relativ ely good accuracy . Regarding mode shapes and participation factors, Fig. 3- 4 present the results of the estimated mode shape and the participation factor for Mode 7 which is an inter-area oscil- lation mode. It can be seen that the estimated mode shape accurately identifies the two groups of oscillation—Generator 4 and 5 are oscillating against Generator 6 and 7. Moreover , the estimated participation factor correctly captures the most influential participants—Generator 5 and 7, implying that the correctiv e control measure should be adopted at these two generators to depress this particular mode. B. Case II: V alidation on High-Order Model with Detailed Contr ol and with Measurement Noise Since the method is dev eloped based on the classical generator model, it is necessary to justify its applicability to higher-order model with detailed control components in real- world power systems. In addition, the PMU data is ine vitably subject to measurement noise in real life, the impact of which needs to be inv estigated. W e still consider the IEEE 39-bus system while all gener - ators are modelled as the third-order models equipped with 4 0.2 0.4 0.6 0.8 30 210 60 240 90 270 120 300 150 330 180 0 (a) The actual mode shape 0.2 0.4 0.6 0.8 30 210 60 240 90 270 120 300 150 330 180 0 Gen 5 Gen 7 Gen 6 Gen 4 Gen 2 Gen 3 Gen 8 Gen 9 Gen 10 (b) The estimated mode shape Fig. 3: Case I: the actual and estimated mode shapes for Mode 7.                     Fig. 4: Case I: the actual and estimated participation factors for Mode 7. exciters and power system stabilizers (PSSs). In addition, Gaussian measurement noises with the standard de viation 10 − 3 for angles and 10 − 6 for speeds are added to the emulated PMU data according to the IEEE standard [18]. Like wise, 450 s noisy PMU data is employed by the proposed mode identification algorithm, and the estimation results for the mode frequencies and the damping ratios are shown in T able III-IV. It can be seen that the proposed technique still provides good estimation results for frequencies and damping ratios of all modes despite of the fact that higher-order models are used and measurement noises are incorporated. This observation is further confirmed by Fig. 5 that presents a comparison between the estimated eigenv alues and the true ones (only the positiv e frequencies are shown). Mode T rue Estimated Estimation frequency (Hz) frequency (Hz) error (%) 1 0.497 0.498 0.29 2 1.267 1.270 0.22 3 1.425 1.413 0.81 4 1.500 1.501 0.08 5 1.650 1.642 0.51 6 1.817 1.797 1.06 7 2.182 2.157 1.15 8 2.052 2.050 0.097 9 2.110 2.082 1.33 T ABLE III: Case II: the estimation results of frequencies for all modes. Particularly , we consider Mode 2—an inter-area oscillation mode with frequency 1 . 267 Hz and damping ratio 2 . 00% . In this mode, Generator 4-7 oscillate against the others. A comparison between the actual and the estimated mode shapes is presented in Fig. 6, showing that the estimated mode shape by the proposed technique well matches with the actual mode Mode T rue Estimated Estimation damping ratio(%) damping ratio (%) error (%) 1 6.84 6.90 0.92 2 2.00 2.06 2.86 3 2.47 2.49 0.94 4 1.39 1.29 7.22 5 3.27 3.29 0.83 6 2.75 2.72 1.07 7 2.68 2.77 3.21 8 0.21 0.27 28.57 9 1.38 1.35 2.17 T ABLE IV: Case II: the estimation results of damping ratios for all modes.               Fig. 5: Case II: a comparison between the estimated eigen v al- ues and the true eigen v alues. shape. In addition, Fig. 7 presents a comparison between the actual participation factor and the estimated one from which we see that the proposed technique accurately estimates the participation factor , indicating that Generator 5 and 9 take the most responsibility in each group for the e xcitation of this mode. 0.2 0.4 0.6 30 210 60 240 90 270 120 300 150 330 180 0 (a) The actual mode shape 0.2 0.4 0.6 30 210 60 240 90 270 120 300 150 330 180 0 Gen 5 Gen 9 Gen 4 Gen 7 Gen 6 Gen 2 Gen 3 Gen 8 Gen 10 (b) The estimated mode shape Fig. 6: Case II: the actual and estimated mode shapes for Mode 2. Note that mode shape and participation factor are different in the sense that mode shape gives the relativ e activity of each state in a particular mode while participation factor presents the influence of each state on a specific mode. Therefore, par- ticipation factor may be more effecti ve on designing damping control measures. In most cases and as seen from the previous results, the generators with the largest magnitudes in the mode shape correspond to the generators with the largest magnitudes in the participation factor . That is the reason for which people assume that mode shape can aid in the mitigation of a modal oscillation by pinpointing the right location where a particular mode is most energetic [9]. Nevertheless, we find that the generator with the largest magnitude in the participation factor may not align with the one with the largest magnitude in the mode shape. This discrepancy indicates that using the mode 5 shape may not ef fectiv ely find the best locations for damping control measures.                                     Fig. 7: Case II: the actual and estimated participation factors for Mode 2. For instance, Mode 6 in the considered 39-bus system is an inter-area mode in which Generator 8 and 10 oscillate against Generator 2 and 9. Fig. 8-9 present the comparisons between its actual mode shape/participation factor and the estimated ones by the proposed technique. It is observed from Fig. 8a and Fig. 9a that Generator 8 has the largest magnitude in the mode shape, yet its component in the participation factor is only the third largest (Generator 2 and 10 are top largest). This difference entails that Generator 2 and 10 hav e more influence than Generator 8 to Mode 6, ev en though state variables of Generator 8 are the mostly excited ones in this mode. Therefore, damping control for Mode 6 may be more effecti vely conducted at Generator 2 and 10 than at Generator 8. Using the mode shape rather than the participation factor may not effecti vely pinpoint the best location to conduct control measures. From Fig. 9b, it is observed that the proposed algorithm accurately captures the most influential participants—Generator 2 and Generator 10 in the mode, although the influence of Generator 10 is slightly ov erestimated. 0.2 0.4 0.6 0.8 30 210 60 240 90 270 120 300 150 330 180 0 (a) The actual mode shape 0.2 0.4 0.6 0.8 30 210 60 240 90 270 120 300 150 330 180 0 Gen 8 Gen 2 Gen 10 Gen 9 Gen 3 Gen 4 Gen 5 Gen 6 Gen 7 (b) The estimated mode shape Fig. 8: Case II: the actual and estimated mode shapes for Mode 6.           (a) The true participation factor            (b) The estimated participation factor Fig. 9: Case II: the actual and estimated participation factors for Mode 6. IV . CONCLUSIONS AND PERSPECTIVES In this paper, we have proposed a novel mode identification method which can accurately estimate the modal knowledge (e.g., frequencies, damping ratios, mode shapes, participation factors) of electromechanical oscillations using ambient PMU data. Particularly , participation factors can be estimated with reasonably good accuracy without the knowledge of detailed network model and parameters. The potential discrepancy be- tween the participation factor and the mode shape is discussed, showing that the proposed technique capable of estimating par- ticipation factor may greatly facilitate ef fectiv e damping con- trol design. Further efforts are expected to dev elop damping control measures utilizing the estimated participation factors. R EFERENCES [1] S. A. N. Sarmadi and V . V enkatasubramanian, Electr omechanical mode estimation using recursive adaptive stochastic subspace identification, IEEE Trans. Power Syst., vol. 29, no. 1, pp. 349-358, 2014. [2] G. Liu, and V . V enkatasubramanian, Oscillation monitoring fr om ambient PMU measur ements by frequency domain decomposition. in IEEE International Symposium on Circuits and Systems, 2008. [3] D. J. Trudno wski, Estimating electr omechanical mode shape from syn- chr ophasor measur ements . IEEE Transactions on Power Systems, vol. 23, no. 3, pp. 1188-1195, 2008. [4] L. Guoping, V . V enkatasubramaniar , Oscillation monitoring from ambient PMU measurements by fr equency domain decomposition . Proceedings of IEEE International Symposium on Circuits and Systems, pp. 2821-2824, May 2008. [5] L. Dosiek, J. W . Pierre, D. J. Trudno wski, and N. Zhou, A channel matching approac h for estimating electromec hanical mode shape and coher ence . In IEEE PES General Meeting, 2009. [6] N. Zhou, L. Dosiek, D. Trudnowski, and J. W . Pierre, Electr omechanical mode shape estimation based on transfer function identification using PMU measurements . In IEEE PES General Meeting, 2009. [7] S. N. Sarmadi, V . V enkatasubramanian, Electromec hanical mode estima- tion using r ecursive adaptive stochastic subspace identification. IEEE T ransactions on Power Systems, vol. 29, no. 1, pp. 349-358, 2014. [8] H. Ghasemi, C. Canizares, and A. Moshref, Oscillatory stability limit pr ediction using stochastic subspace identification. IEEE Transactions on Power Systems, vol. 21, no. 2, pp. 736-745, May 2006. [9] L. Dosiek, N. Zhou, J. W . Pierre, Z. Huang, and D. J. Trudno wski, Mode shape estimation algorithms under ambient conditions: A comparative r eview . IEEE T ransactions on Power Systems, vol. 28, no. 2, pp. 779- 787, 2013. [10] S. Mendoza-Armenta, and I. Dobson, Applying a formula for g enerator r edispatch to damp interar ea oscillations using synchrophasor s. IEEE T ransactions on Power Systems, vol. 31, no. 4, pp. 3119-3128, 2016. [11] X. W ang, J. Bialek, and K. Turitsyn, PMU-based estimation of dynamic state J acobian matrix and dynamic system state matrix in ambient conditions . IEEE Transactions on Power Systems, June 2017. [12] X. W ang, K. T uritsyn, PMU-based estimation of dynamic state J acobian matrix. IEEE International Conference on Circuits and Systems (ISCAS) 2017. [13] O. A, Theresa, M. L. Crow , An analysis of power system transient stability using stochastic energy functions. International Transactions on Electrical Energy Systems, vol. 23, no. 2, pp. 151-165, 2013. [14] C. Gardiner , Stochastic Methods: A Handbook for the Natural and Social Sciences. Springer Series in Synergetics. Springer, Berlin, Germany , 2009. [15] X. W ang, Estimating Dynamic Load P arameters fr om Ambient PMU Measur ements. In IEEE PES General Meeting, 2017. [16] N. Zhou, S. Lu, R. Singh, and M. A. Elizondo, Calibration of reduced dynamic models of power systems using phasor measurement unit (PMU) data. North American Power Symposium (NAPS), 2011. [17] F . Milano, An open source power system analysis toolbox. IEEE T ransactions on Power Systems, vol. 20, no. 3, pp. 1199-1206. [18] IEEE Standar d for Synchrophasor Measurements for P ower Systems- Amendment 1: Modification of Selected P erformance Requir ements . IEEE Std C37.118.1a-2014 (Amendment to IEEE Std C37.118.1-2011), pp. 1- 25, April 2014.