Forced Oscillation Source Location via Multivariate Time Series Classification

Forced Oscillation Source Location via Multivariate Time Series Classification Yao Meng 1,2 , Zhe Yu 1 , Di Shi 1 , Desong Bian 1 , Zhiwei Wang 1 1 GEIRI North America, San Jose, CA 95134 2 North Carolina State Un iversity, Raleig h, NC Emails: ymeng3@nc su.edu ; { zhe. yu , d i.shi , desong.bia n, zhiwei.wang } @geirina.net Abstract — Precisely locating low-freq uency oscillation sources is the prerequisite of suppressing sustain ed oscilla tion, which is an essential guarantee for the secure and stable operation of power grids. Using synchrophasor measur ements, a machine learning method is p roposed to locate the source of forced oscilla tion in power system s. Rotor angle and activ e po wer of each power plant are u tilized to construct m ultivariate time series (MTS). Applying Mahalanobis distance metric and dynam ic time w arping, the distance between M TS with different phases or lengths can be appropriately measured . The obtained distance metric, representing characteristics during the tra nsient phase of f orced oscillation under d ifferent disturbance sources, is used for offline classifier trai ning an d onli ne match ing to l ocate t he disturbance source. Simulation results using the four-machine two-area system and IEEE 39-bus system indicate that th e proposed locati on method ca n identify the pow er system force d oscillation source on line w ith high accuracy. Index Terms — Oscillation location, multivariate ti me series, Mahalanobis distance, dynamic time warping, machine learning I. I NTRODUCTI ON Low-frequen cy oscillation reduces po wer trans mission limit and may lead to damage of system eq uipment. It severe ly threatens the security and stabilit y of large-sca le interconnected power systems in real-time operations . Insufficient d amping of a syste m results in a majorit y of low- frequency oscillation, which can be suppressed by tun ing parameters o f po wer system stab ilizer or intertie line co ntrols . However, forced oscillation s caused by reso nance have been discovered recently [1 ] in p ower grid s . T he po sterior analysis shows the system was well damped when the oscillation occurred. Moreover, disturbance with a frequency approximating the intrinsic system freq uency w as i njected into the power system so mewhere. The resonance excited and ev en the s mall dist urbance could amplify and spread rapidly in the whole power system. T he tra ditional remediatio n actions, for example p utting p ower system stabilizer s into operation, etc. , are not applicable for suppr essing such osc illation s. T he most effective way to quench forced oscillations is to re move t he disturbance rapidl y and accurately. Locating the oscillation sou rce is the pr erequisite to eliminating the disturbance . W ith massive Phas or Measurement Units (P MUs) in stalled, m onitoring dynamic behaviors of the po wer system has beco me possible. T aking advantage of PMU m easure ments, the research community has proposed several forced oscillatio n so urce loca tion methods. Integrating PMU m easureme nts at different places and the wave speed map, a traveling wave based method is prop osed in [2]. B ased on the transient energy function, the energy flow direction can be calculated to locate th e f orced oscillation source, which has good location performance [ 3]-[5]. Authors of [6] locate the oscillation source by estimating t he mode shape, which represe nts the relative magnitude and p hasing of the oscillatio n throu ghout the system . Using tie -line power signal d uring o scillations, the real-time ap proximate entropy value i n a continuous time interval, which is correspo nding to the locatio n of disturbance in p ower systems , can be calculated [7]. I n [8], a machine le arning app roach is proposed. Measurement signals duri ng fo rced oscillations ar e mapped to a mu l ti-dimensio nal CE LL, and decision tree is utilized to identify the charac teristic par ameters of the CELL corresponding to d ifferent oscillation so urces. The discoveries and investigations of t hese methods provide a better underst anding of applying P MU measurements to locate the disturbance source. Nonethele ss, some appro aches suffer fro m a lack of t heoretical support. Moreover, m ost literature assumes the low-freque ncy oscillation ca n b e immediately detected after it occurs and utilizes the i nfor mation at t he beginning pha se of oscillation s , which are not pr actical in general. A machine learning based m ethod to locate the forced oscillation sources is proposed in this p aper. Rotor angle and active po wer of eac h ge neration m easured by PMUs are utilized to co nstruct multivariate time series (MT S). Mahalanobis distance and dynamic time warpin g (DT W) are applied to represent the distance b etween out - of -sync MTS. With t he obtained d istance metric, offline trai ning a nd online matching can b e conducted . In this appr oach, MT S with different len gth can be appropriately compared which relaxes the assumption of accurate detection of the beginning of oscillations. Case studies o n the four-machine two-area system and IEE E 39 -bus s ystem indicate the effectiveness of t he proposed ap proach . In addition, the e ffect of time dela y in oscillation detection o n locati ng accurac y is discussed . This work is funded by SGCC Science and Technology Program under contract No. 54 55HJ160007. II. M ULTIVA RIATE T IME S ERIES C LASSI FICATION A machine learning approach that utilizes P MU measurem ents of generator rotor angle and active power in power grids to locat e the source of forced oscilla tion geographica lly down to the substation level is proposed. It is assum ed that at least one PMU measu remen t is av ailable f o r each power p lant, and the model of the p ower system is also available . Fig.1 illu strates th e proposed ap proach proce dure. When forced oscillati on o ccurs, the r esponse of gene rators tends t o be distinct for diff er ent oscillat ion so urce locations . Based on this intuiti on, dy nam ics of generators under vari ous s cenarios of forced o scillati o n source locations are generated thr ough tim e domain sim ulation . S imu lated PMU m easurements of rotor angles and activ e power o f generators are utilize d to construc t signature multivariat e tim e series. Maha lanobis dist ance is trained offline through metric learnin g. When a forced oscillati o n is detected, the same coefficients are measured and compared with th e trained classif ier to determin e the sou rce location of the oscil lation. Dynam ic time warping is applied to handle the out- of -sy nc between testing data sets an d tr aining data sets due to the oscil lation detec tion del ay. Onli ne Matching Detect F orc ed oscil lation Obtain PMU Measurements Establis h T esting Datasets Locat e the oscil lation source Appl y K-Neares t Neig hbour Proce ed Dynamic T im e W a rping Offli ne T rai nin g Conduct T im e Dom ain S imul a tion Establis h Trai ning Datasets Obtain Mahalanobis Matrix Fig.1 Flow chart of forced oscill ation source iden tification based on multivariate time se ries classification . A. Ma halanobis Distance Mahalanobis Dis tance is a stan dard measure of the distan ce, which is characteri zed by a symm etr y Po sitive Sem i-Definite (PSD) matrix M . For two vectors x and y , where x , y ∈ R d , the squared Mahalanobis distance betw een x and y is defined as follow s: ( , ) ( ) ( ) T M d x y x y M x y    (1) If M = I , Mahalan obis distance degenerates to the s tandar d Euclidean d istance . Since Mahalan o bis distanc e considers corr elation among different va riables, an accurate relati o nship betw een variables and l abels of M TS can be established. It has two essential functions. T he fi rst one is to rem ove correlati on s among different vari ates and t o map the original space in to a new coordinate system . T he second one is to assign w eights to new variates [9]. With these two functions, Mahalan obis dis tanc e can measu re the distance b etween ve ctors effective ly. Next, the Mahalanobis m etric will b e extended to measure the distance between two multivariate time series. Given two MTS X and Y , 12 12 12 ( 1 ) ( 1 ) ( 1 ) ( 2) ( 2) ( 2 ) ( ) ( ) ( ) p p p x x x x x x X x h x h x h         (2) and 12 12 12 ( 1 ) ( 1 ) ( 1 ) ( 2) ( 2) ( 2 ) ( ) ( ) ( ) p p p y y y y y y Y y h y h y h         (3) where p is the num b er of featu res , h is the num ber of samplin g points . The lo cal dist ance measure is ex pressed as ( , ) ( ) ( ) i j i j i j T M d X Y X Y M X Y    (4) where X i represents the i th row in X , and Y j represents the j th row in Y . Then the distance betw een multivariat e time series X and Y is d efined as 1 ( , ) ( , ) h jj MM j D X Y d X Y    (5) B. Metric Lea rning The Mahalan obis metric represen ts the relev ance of tw o time series. T he goal of metric learning is utilizing the labeled training data to find an appropriate M such that the Mahalanob is distance emphasizes the relevant features while decreases the effect of irrelevan t featu res [ 10 ]. For any given triplet label { X , Y , Z } , where X and Y are in the sam e class while Z is in a different o ne , the Mahalanobis distance between instan ce X and instance Y should be smaller than the distance betw een X and Z . The metric lea rning probl em can be formu lated as an optimizati o n probl em as f o llow s: min . . ( , ) ( , ) , { } PSD M MM 1 s t D X Y D X Z X,Y, Z M       (6) Here, ρ >0 denot es the desired margin. The obje ctive of (6) is to find a PSD matrix to satisfy all triplet constraints { X, Y, Z }. T he number of trip let label [ 11 ] cons traints is the cubic of the number of the training samp les. In [12], th e authors have proposed a st r ategy to choo se th e m ost useful tripl ets f or training . To solve the optimizati on proposed in (6), an iterativ e process is proposed, wh ich can be w ritten as: Algor i thm 1 Train Mah alanobis Mat rix th rough M etric Le arning 1: I nitializ e Ma halanob is ma trix M 0 2: For i = 1:M AX 1 3: For j = 1:M AX 2 4: Choo se the most u seful triple t s co n str aint { X j , Y j , Z j } 5: If co nstraint ( 6) is v i ola t ed, calculate loss f unction l ( M j ) 6: Update Mahala nobi s ma t rix M j 7: End fo r 8: Calcul ate to tal lo ss func tion L i = ∑ j l(M j ) 9: Chec k whether L i is le ss t han th res hold , if y es, brea k 10: End fo r MAX1 and MA X2 are set t o lim it the number of ite rations. If th e trip let constraint is viola ted , the loss function is define d as : ( ) ( , ) ( , ) M j j M j j l M D X Y D X Z     (7) How to update Mahalanobis matrix is the key problem in Algorithm 1 . The Mahalanobis matrix M j should be update d to reduce th e valu e of l oss fu nction. What’s more, to avoid unstable learning process, a regularizati on term which restric ts the matrix diverg ence between two differen t iterations should be added to the m etric learnin g objectiv e function. So the Mahalanobis matrix updating equati o n can be expresse d as: 1 0 arg m in ( , ) ( ) jj M M div M M l M    (8) where λ is the regulari zation param eter balancing the loss function l ( M i ) and the regulariz ation function div ( M j ,M j+1 ). The regulari zation fun ction div ( M j ,M j+1 ) measu res the matri x divergenc e. It can be exp ressed as: 11 1 1 1 ( , ) ( ) lo g( de t( )) j j j j j j di v M M t r M M M M n        (9) where tr () denotes the trace of a matrix, n is th e dimension of M [1 3 ]. We solv e (8) in an iterative way. T o ensure the obtained Mahalanobis matrix M j +1 is a PSD matrix in each iterat ion, we require 1 ( ) 0 0 TT jj j j j j j P P Q Q M           (12) Several tools can s olve this st andard linear matrix inequalit ies (LMIs). If the o btained resul t is 𝜆 𝑗  , λ j ∈ [0, 𝜆 𝑗  ] ensures tha t updated M j +1 is a PSD matrix. T herefore , w e so lve LMIs f irst in each iter ation and selec t λ j in t he feasib le rang e. Given λ j , equation (8) reaches its minimum when its gradient is zero. By setting the gradient of (8) to be zer o, we get: 11 1 ( ( )) TT jj j j j j j M M P P Q Q       ( 10 ) where P j = X j -Y j , Q j = X j -Z j . We can solve ( 10 ) by applyin g Woodbury matrix i dentity, th e ite r ative ex pression o f M j is: 1 1 1 () () TT j j j j j j j TT j j j j j j j j j j j j j j j M M P I P M P P M M Q I Q Q Q                      ( 11 ) where γ j =( M j -1 + λ j P j P j T ) -1 . Through this process , Mahalanobis matrix can be u pdated if triplet constraint is violated. When calculate d total loss function is smaller than a predefined thresh o ld or MAX 1 is reached, the algorithm termin ates. Metri c lea rning is the most compute- intensive part of th e proposed a pproach, but it can be done offline, which is not a t ime-critical t ask. C. Dynamic Time Warping When apply ing multivariate tim e ser ies classifi cation , it is unreasona b le to assum e the begin ning and ending time points of interes t c an be correctly identifie d , especial ly d uring the later deploym ent. In th e o scillati on locatin g context , det ection of the beginning of the f orced oscillati on event is not guarantee d . Thus, w e need to deal with time series analy sis with differe nt phases and leng ths. Dy namic Time War ping (DTW) is an algorithm w hich con d ucts nonlinea r mapping of one tim e ser ies to another b y minimizin g the d istance [1 4] . T hrou gh calculating the optim al warp p ath, tw o tim e ser ies w ill be extende d and placed in one- to -one cor responden ce, which m akes their similarity can be measur ed eas ily. Given two time s eries, Q ( i ) , i = 1,2,…,m and C ( k ) , k = 1,2,…,n , define a n optim al warp path W as () ( ) , 1 , 2, ... , () wj Q W j j s wj C      ( 13 ) where w Q ( j ) ∈ [1 , m ] denotes the index in se quence Q , w C ( j ) ∈ [1, n ] denotes the index in sequence C, and s is the leng th of the warp path . ( w Q ( j ) , w C ( j )) T means the w Q ( j ) th element of Q and the w C ( j ) th elemen t of C co rrespon d to each other. To reduce the number of paths during the search , a vali d warping path sh ould s atisfy several well- known conditi ons . Boundary condition ensures all in dices of each time series are used in the w arping path. Continuity c ondition requires the warping path to be made of o nly adjacent cells. Moreover, monotonic ity conditi o n restrict s the feasibl e warping pa th only increase monotonically . Th ese thr ee cond itions can b e expressed as ( 1 ) ( 1 , 1 ) ' ( ) ( , ) ' ( , ) ( ( ), ( )) min( ( 1 , 1 ), ( 1 , ), ( , 1 )) W W s m n D i k d Q i C k D i k D i k D i k               (1 4) where d ( i,k ) is the distance f ound in the current cell, D ( i,k ) represents the min imum sub war p path distance , and the l ength of the w arping path s ∈[max( m, n ) ,m + n ] . The optimal warping path ca n be found usin g dynam ic programm in g . Af ter all o f the elements in the distance matrix D ( i,k ) are ca lcula ted , the correspon d ing w arp ing path is the optimal warping path W . Fig.2 illustrat es an example of optimized warping path betw een two given time series. The shadow ce lls denote the correspon d ence r elations hip betw een these two time se ries. Traditional dynamic time warping al gorithm is only applicable to univariate ti me series. T o utili ze d ynamic time warping in multivariate ti me series and apply Mahalanob is distance, the local d istan ce d ( i , k ) is defined as () , ( , ) i j i j M d X Y d X Y  (15) where X i repr esents the i th row in (2), and Y j represents the j th row in (3). Fig.2 . An op timal warping pa th. Now the distance bet ween ou t- of -sync MTS is measured correct ly . K -nearest n eighbors algorithm ( k - NN) is u sed for classification. The input of k -NN is the multivariate time series constructed by PMU measure ments of generator rotor angle and active po wer; the output i s a class membership. An ob ject is class ified by a majority vote of its neighbors, with the object being assigned to the class most co mmon a mong its k nearest neighbors. III. N UMERI CAL E XPERIMENTS In this section, num erical result s of the proposed method are presented . T he PMU measurem ent d ata sets u sed in case stud ies are generate d by tim e-domain simu lations using DSATools TM . The sam p ling rate is 25 Hz. A sinusoidal signal serve d as the forced oscillation disturbance is in jected int o excitati o n sy stems . Rotor angle an d active pow er o f each generat or serve as features . For trainin g data sets , the t ime series begins at 0s and en ds at 15s , w hile for testing data sets, there is a delay d correspon d ing to the tim e needed for o scillati on d etection. Therefore, the testing time series starts d se conds l ater than the b eginning of the o scill ation and the w indo w s ize of time is 5 seconds. T he cons tructe d multivari ate time series is the same as (2), where p equals to two tim es the number of generators in th e power gri d , h eq uals to 25× 15=375 for training sequences, an d 25× 5=125 for the testing s equence. A. F our-mach ine Two-area Model The detail ed model param eters can be found in [15]. All generators are in a second- o rder m odel, while the load is in a constant impedance m od el . W ith the SSA T softw are, a detailed modal analysis on the system is con duct ed . The analysis results show th at a ll natural modes h ave reasonably good damping and the sys tem has a n atural mode at f 0 =0.62 08Hz . A sinusoidal signal Δ ref = k× si n ( 2πft ) emulating the oscillati o n distu rbance is a dded to the refe rence sign al of excitati on sy stems, w here k= 0.03.The oscillation begins at t =0s. For each instance, o nly one of the generato rs is inject ed with the oscil lation d isturban c e. To gene rate enough n umber of valid sam ples, the sy stem load varies randomly betw ee n 90 percent and 110 percent of the original load. PSAT is applied to solve the power flow equations and the infeasible load conditions are remov ed. In reality, system param eters a re tim e- varyin g and random , which leads to the discrepancy between the model and the real system . What’s more, the frequency of oscillati o n source can cause resonance fluctua tes with in a specific range. To em ulate practical situ ation, th e damping factor of each gener ator w hich influences the low-frequenc y oscillati o n most is chosen randomly in [0,4] . T he fre quency of oscillati o n source f fluctu ates in the range of 90 percent to 110 percent of f 0 . Since only one of the generators acting as the oscillation source, there are four scenari o s in total. Now the oscillati on source locati ng problem converts to a multiclass classific ation problem . Follow ing rules mentione d above, 800 samples (200 samples for each scenario) are generated. T he ratio of tra ining data sets and testing data sets is 1:1. To make the simulati o n more practical , G au ssian White Noise with the sign al- to -noi se ratio as 13dB is su perimposed to the sim ulated PMU measurem ents. As sume the low-frequency oscill ation can be detected in 3 seconds , so the te sting tim e series begins at d =3s . The d imension of train ing and testing mu ltivariate time series is 375 - by -8, 125 - by - 8 respectively . In the metric learni ng process, dynamic triplet co nstr aint building strategy is applied to select triple t constr aints. With the obtain ed Mahalan obis matrix, the distance between training samples and testing samples can be calculated . Af ter that , k -neares t neigh bor classifi cation is applied t o sele ct the label of th e n earest train ed classifi er as the catego r y of the test s ample with k =1 . Table I presents the accu racy of o scillati on source loc ation . According to Ta b le I, the a ccuracy of a ll four situa tions is 100 % ， which in dicates th e effectivene ss of th e propos ed approach. T ABLE I P ERFORMANCE T EST OF O SCILLATION S OURCE I DENTIFICATION FO R F OUR -M ACHINE T WO - AREA M ODEL Ge nerator wit h disturb ance Corre ct Error Accur acy/% G 1 100 0 100 G 2 100 0 100 G 3 100 0 100 G 4 100 0 100 B. IE EE 39-bu s System The detail ed model param eter s can b e found in [16]. Generato rs are in a fou r th-orde r mo del, w hile the const ant impedance load m odel is adopt ed. Based on the de tailed modal analysis of the sy stem, a natural mode with the frequ ency at f 0 =1.32 17Hz exists . Similar to the four-m achine tw o-area c ase, a sin usoidal signal Δ ref = k × si n ( 2πft ) w ith k =0.6 is adde d to the referen ce signal of excitati o n system s. Si nce each time only one generat o r acting as the f orced o scillati on source, there are ten scenarios in total. The oscillati on d istur bance adds to the system at t =0s . Considerin g randomn ess in sy stem load, the frequency of forced o scill ation sourc e and damping factor of each generator, 4000 samples are gene rated. The rati o of training data sets and testing d ata sets is 1 :1. Gaussi an W hite Noise with the signa l- to -noise ratio equalin g to 13dB is superimposed to the simu lated PMU measurem ents. T he testin g time seri es begins at d =4s. Table I I pres ents the oscillati o n sou rce location perform ance of p ropose d approach. Another machine learning approach , CELL&Decisi on tree ap proach from [8], is employ ed as comparison . From table II, th e overall accuracy of proposal is 97.8%, w hile the ac curacy o f s everal scena rios reaches 100%, sa tisfy ing the accurate positioning request of engineerin g practice. In each scenario, the proposal outperfo rms the CEL L&Decis ion metho d. T ABLE II P ERFORMANCE T EST OF O SCILLATION S OURCE I DENTIFICAT ION FOR I EEE 39 - BUS S YSTEM Ge nerator wit h distu rba nce Accur acy of propo sal /% Accur acy of CELL &Decisio n/% G 1 100 95.7 G 2 100 92.8 G 3 100 98.7 G 4 100 99.3 G 5 92 91.3 G 6 100 100 G 7 100 97.3 G 8 92 90.1 G 9 94 93.3 G 10 100 95.7 Ave rage 97.8 95.4 C. Influence of Oscillatio n Detection Dela y Here the influence o f time del ay in oscillat ion d etection is analyzed. When forced oscillation occurs, software o r sys tem operators need som e time to detect it . In our work, test ing time series b egins at t=d second t o emulate the time requir ed for oscillation detection. For I EEE 39 -bus s ystem, we ha ve conducted simulation s with different value s of d . Fig. 3 illustrates the relationship between the lo cation acc uracy and the dela y d . As shown in Fig. 3, as delay d increases, the location accurac y declines. T his is because the pr oposed approach localizes oscillation sources mainly based on the dynamic character istics of features i n the transient proces s. When the oscillatio n tends to be stable, the characteristic of features b ecomes indistinct. So the sooner the oscillation is detected, the higher is the acc uracy . Fig.3 The re lationship betwe en location accuracy and the time de lay IV. C ONCLUSION In this work, a multivariate time series classifica tion method to loca te the forced oscillation source s in power systems has been proposed . Mahalanobis distance of MTS constructed b y roto r angle a nd ac tive po wer of e ach generat or is o btained through metric lea rning. D ynamic time warping is applied to find the optimized warping path for time series with different p hases o r lengths. After that, t he r eal-time measurements ca n be utilized to determine the forc ed oscillation source by compar ing t he Ma halanobis distance between the te sting data and t he training d ata . B ecause MT S integrates all information d uring the transient pr ocess, i t can represent the force d oscillation characteristic ca used by different oscillation disturbance ver y well. Experi mental results show the p roposed ap proach has hig h location accuracy. Considerin g practical ap plications, the relationshi p between oscillation d etection delay and detection accuracy is investigated . The p roposed approach works well for forced oscillation . In practice, poorly damped oscillation is another principal mechanism of o scillations. Ho w to utilize the prop osed approach to deal with poorl y d amped oscillation source locati on problem remains to b e further investigated . R EFERENCES [1]. D. Song, X. Yang, Q. D ing, et al., “A su rvey on low-frequency oscillation in large -scale inter-conne ct ed powe r grid and its control measures,” Power System Tech nology , vol. 35, no. 10, pp. 22- 27 , 2011. [2]. J.S. Thor p , C.E. Se yler, and A.G. Phadke , “ Electromechanical wave propagation i n large electric pow er systems,” IEEE Trans. Cir cuits Syst. I, Fundam. The ory Appl. , vol.45, no. 6, pp.614 – 622, Jun. 1998. [3]. L. Chen , Y. Min , and W. Hu, “ An energy-based method for location of power system oscillation so urce,” I EEE Trans. Power S yst., vol.28, n o.2, pp. 828 – 836, May 2013. [4]. C. Jing, J. D. McCalley, and M. Kommare ddy, “ An energy approach to analysis of inter area oscill ations in pow er systems,” IEEE Tra ns. Power Syst ., vol .11, no.2, pp.734-740, M ay 1996. [5]. R. B. M ye rs, and D. J. Trudnowski, “ Effects of forced oscillations on spectral-based mode-shape e stimation,” I EEE P ES Ge n eral Meeting , Vancouver , BC, Canada, pp. 21 – 25 , 2013. [6]. L. Dosiek, N. Zhou, J.W . Pierre, et al., “ Mode shape estimation algorithms under ambient c onditions: a comparative review,” IEEE Trans Power Syst. , vol.28, no.2, p p. 779 – 787, Sep. 2012. [7]. X. Yang, D. Liu, Q. Liao, et al. , “ Application of dynamic approxim at e entropy in re al- time detection of forced pow er oscillation,” Powe r System Technolo gy , vol.36, no.11, pp.129-133, 2 012. [8]. C. Jiang, J. Liu, Y. L i u, et al., “Onl ine forced oscillation d etection and identification based on wide are a me asurement system and CELL theory ,” Electric Power Automation Equipment , vol . 35 , no.2, pp.125- 132 , Feb. 2015. [9]. J. M ei, M. Liu, Y. Wang, and H. Gao, “Learning a M ahalanobis distance - based dy namic time warping measure for multivariate time se ries classification,” IEEE Tr ans. Cybern. , vol.46, no.6, pp.1363-137 4, May 2015. [10]. M. Liu, and B. C. Vemuri, “A robust and efficient doubly regular ized metric learning appro ach,” in Proc. of the 11th European conferen ce on Computer Visi on. Springer- Ve rlag , 2 012, pp. 646-659. [11]. J. B i, D. Wu, L. Lu, et al. , “ Adaboost on low-rank psd matrices for metric learning,” in IEEE Conference on Computer Vision and Patter n Recognition . I EEE, 2011, pp. 261 7 – 2624. [12]. J. Mei, M. Liu, H. R. Karimi, and H. Gao, “Logdet d ivergence based metric learning with triplet constraints and its applicatio ns,” IEEE Trans. on Image Process. , vo l. 23, no. 11, pp . 4920 – 4931,20 14. [13]. B. K u lis, M. Sustik, and I. Dhillon, “L earnin g low - rank kernel matrices,” in Proceedings of the 23rd international conference on Machine learning. ACM, 2 006, pp. 505 – 512. [14]. M. Muller , “Dynamic time warping,” Informati on Retrieval fo r Music and Motion , pp. 69-84, 2007. [15]. P. Kundur, Power system stability and control . The EPRI power system engineering se ries. McGr aw-Hill, 1994. [16]. M. A . Pai . Energy function an alysis for po wer system stability. T h e Kluwer internation al series in enginee ring and computer science. Pow er electronics and power systems. Boston: Kl uwer Academic Publisher s, 1989.