Probabilistic Robust Small-Signal Stability Framework using Gaussian Process Learning

1 Probabilistic Rob ust Small-Signal Stabilit y Frame work using Gaussian Process Learning Parikshit Pareek, and Hung D. Ng uyen ⋆ Abstract —While most power system small-signal stability as- sessments rely on the reduced Jac obian, which depends non- linearly on the states, uncertain operating points intro duce nontrivial hurdles in certifying the systems stability . In this paper , a nov el probabilistic robust small-signal stability (PRS) framewo rk is dev eloped f or a power system based on Gaussian process (GP) learning. The proposed PRS assessment prov ides a robust stability certificate for a state subspace, su ch as t h at specified by the error bounds of the state estimation, with a giv en probability . With such a PRS certificate, all i nner poin t s of the concerned subspace will be st ab l e with at least the corres pondi ng confidence level. T o this end, the b eha vior of the critical eigen value of the r educed Jacobian with state points in a state subspace is learned usin g GP . The proposed PRS certificate along with the Subspace-b ased Search and Confidence-based Search mechanisms constitute a holistic framework catering to all scenario s. The pr oposed fra mework is a p owerful approach to assess th e stability und er uncertainty because it d oes not require input uncertainty distributions and other state-specific input-to- output approximations. Further , th e critical eigen va lue behavior in a state subspace is analyzed using an upper bound of the eigen value variations and their inferences are discussed in detail. The re sults on three-machine n ine-bus WSCC system show that the p roposed certificate can fin d the robust stable state subsp ace with a giv en probability . Index T erms —Probabilistic Robust Small-Signal St ab i lity (PRS), Gaussian P rocess (GP) Learning, Stability und er Uncer- tainty I . I N T RO D U C T I O N The small-signal stab ility is an integral p art of power system stability assessment, refer ring to th e system’ s c a p acity to withstand small-distur bances and remain in synchro nism [1], [2]. Further, stab ility asses sment under un certainty is a crucial issue with the r e cent tren d of integration of uncertain renewable sources and mobile e lectric vehicle loads. The necessary and sufficient conditio n fo r small-signal stability is that all eigenv alues of reduced system Jacob ia n must have a negativ e real par t. This m ethod of stability assessment requir es ev aluation of all or a set of eigen values of redu ced Jacobian , which inv olves a matrix in verse. The red uced Jacobian based stability assessment work s for a specified operating state, but becomes inapplicable in th e presence of state unce r tainties, mostly due to the non-line arity induced from the power flow Jacobian matrix inv ersion. Some imp o rtant, r e lev ant works on D AE and stability a re presented in [ 3]– [10]. ⋆ Correspondi ng Author Authors are with School of E lectri cal and Electroni cs Engineerin g, Nanya ng T echnological Univ ersity , Singapore. e-mail: pare0001@ntu.edu.sg, hunghtd@n tu.edu.sg Parik shit Pareek a nd Hung Ng uyen are supported by NTU SUG, MOE, and NRF fundings. The pr obabilistic small-signa l stability ( PSSS) assessments have been developed to deal with the issue of stability u nder uncertainty [ 1 1]–[13]. Th e co re id ea is to der iv e the pr obability density fu nction ( pd f ) of a n uncertain system o utput, PSSS measure like critical eigenv alue or minim um da m ping ratio (MDR), based on a known uncertain ty of inpu t state [1 4]. The methods of PSSS assessment br oadly fall under two categories: i) num erical methods, and ii) an alytical methods. In nume r ical app roaches, the Monte Carlo simulatio ns (MCS) c a n be employed to determ ine the pd f of the various stability indices [15]–[17]. T o improve computationa l p e rfor- mance, metho ds such as quasi-Monte Carlo [18] have also been used to assess PSSS. Nevertheless, an essential requ ire- ment of a large n umber of simulations leads to com putationa l burden too high to be imp le m ented for detailed p robab ilistic studies. Analytical methods, free from the imposition of p a rametric output distributions, are also used for PSSS assessment. Th ese methods includ e probab ilistic collo cation metho d ( PCM) [1 9], point estimate metho ds ( PE M s) [20], analytical cumulan t- based app roaches [21]. A comp rehensive co mparative stud y is p resented in [ 14]. Th ese meth ods suffer fr om issues related to com p lex f ormulatio n and inaccu rate first and second-or der approx imations of sensitivities of eig en values. The detailed revie w o f these tech niques and approaches can be fo und in [22]–[25]. Regardless of the metho d employed , all th ese PSSS assessment metho ds tried to solve specific uncertain ty issues considerin g a typical pd f for u ncertainties in the wind , solar , or load. Th e m odeling er ror in inpu t u ncertainty d escriptions and effects of assumptions taken get p ropaga te d to output pd f . Further, these metho ds do not provide insight into critical eigenv alue behavior in state-space and are mostly limited to one or two-dimen sional subspace at a time. In this p aper, we pr esent a novel prob abilistic ro bust small- signal stability (PRS) framework for small-signal stability assessment in a state subspace. This state subspace can be defined based on measurem ent error s or th e level of external disturbanc e s that move the system within such a sub-domain . Further, we de fin e PRS a s fo llows. If all the in n er poin ts of a subset of the state- sp ace are small-sign al stable with a t least a g iven le vel of confidence, then such a subset is called probab ilistic robust stable subspace with re sp ect to such a confidenc e level. In other words, PRS is co ncerned th a t for a given uncer tain state- subspace X in state space, what the probab ility with that any state point x ∈ X is stable. T o the best of our knowledge, this is the first instance when the PRS cer tificate ha s been d ev eloped an d pro posed using GP learning for a p ower system. Different from some existing Monte Carlo- based approa c hes, our proposed PRS framework 2 does n ot ch aracterize how m a ny stable operating points in X . Specifically , o ur g oal is to quan tify the p robability with that any op e r ating po int x ∈ X is stable. An advantage of this PRS is th a t it provides a non -parametr ic, computatio n ally efficient, and less complex mod e ling alternative for pro babilistic small- signal stability assessment. More impor tantly , th e prop osed PRS fram ew ork does not require un certainty modelin g , such as a prior distribution of uncertain inpu ts, an d hence provid es a ge n eralized framework for the stability und er uncertain ty . Th e main focus of this PRS framework is provid ing a PRS certificate for a given state subspace being PRS. Along with the c e rtificate, we formulate Subspa c e-based Sear ch and Confi dence-b ased Sear ch prob- lems for cases wherein PRS cer tification can not be verified . The framework is b uilt upo n the Gaussian p r ocess upper confidenc e bound ( GP-UCB ) search algorith m [26]. The GP- UCB is used for sam pling the state points inside state subspace X to learn the b e h aviors of th e critical eig en values, wh ich are closest to the imaginar y axis for a small-signal stable system . The main con tributions o f the work are as follows: • Defining pr obabilistic robust small-signal stability (PRS) and dev elopin g a certificate fo r a subspace to be PRS with a given confidence le vel. • De velopment of a PRS f ramework with Su b space-b a sed Sear ch and Confiden ce-based Searc h fo r a state subspace that is no t satisfyin g the PRS c ertificate with a given probab ility . • De veloping a novel and fast GP learning scheme to learn and analyze the cr itical eig en value behavior in a multi- dimensiona l state- su bspace. The objective o f this work is to p resent the novel idea of PRS for the power system. W e first build up the backgro und in Section II by providing a br ief r evie w of small-sign al stability assessment, GP , and G P -UCB method s. From this, we build the PRS f r amework, which is presen ted in Section III. This section in cludes the main results on PRS ce rtificate with details on the S u bspace- b ased Sear ch and Confid ence-ba sed Sear ch mechanisms. The simulations an d discussions of resu lts are provided in Section IV while con clusions a r e dr awn in Sectio n V with future scop e . W e use stability to r e fer to small-signal stability for b revity . I I . B A C K G RO U N D The prop osed PRS framework in this work has three backgr o und building blocks. This section p resents these thr ee blocks, namely: Power System Mod eling and Stability Assess- ment, GP and GP-UCB. A. P ower S ystem Modeling an d Stability Assessment Consider the power system d ynamics that c an b e expressed as semi-explicit D AE as: ˙ x = f ( x , y ) 0 = g ( x , y ) (1) Here, x , y are dyn amic and algebraic variable vector s respec - ti vely and f ( · ) and g ( · ) are sets of d ifferential and algebraic equations respectively , expressing the power system’ s behav- ior . Th e D AE system can be linearized an d expressed a s: δ ˙ x = A δ x + B δ y 0 = C δ x + D δ y (2) Con ventionally , to assess the stability of the DAE system, we rely on the red u ced Jacob ian matrix ob tained by eliminat- ing the algebraic v ariables y [1]. With an in vertible matrix D , we ha ve the reduced Jacobian J r = A − B D − 1 C . The linearized D AE system is small-sign al stable at the b a se oper- ating po int x ⋆ if an d only if th e re duced Jacobian evaluated at x ⋆ is a Hur witz stable matrix , i.e . , max i ( ℜ ( λ i ( J r ))) < 0 , where λ i ( J r ) be the i th eigenv alue of J r [27]. W e use λ c ( x ) to denote the fu nction d escribing the real part of th e maximum eigenvalue of J r ( x ) as the state x varies. The refore, the beh avior of this fun ction λ c ( x ) d ictates the small- sig n al stability of the system . In the power system, som e attempts have been made to learn and estimate the movement of critical eigenv alue [28], [2 9]. The numerical results o f these work s indicate that eigenv alues of power system Jaco bian are continuo us in continu ous system state space. An imp ortant research direction related to these numerical app roaches is based on eige nvalue sen siti vities with respect to states [ 28], [29]. Th e downsides of th is directio n have been discussed in [27] an d oth ers. This paper , howev er, does no t rely on sensitivities but attempts to lea rn the function representin g the critical eigenv alues wh en the states chang e. The learned fu n ction can benefit other oper a ting p rocedu res in power systems concerning eigenv alues. In ord er to apply GP , we rely on the fact tha t the eigenv alues of a matrix are continuous f u nctions of its entries. This proper ty is sho wn in Theo rem 5.2 [30]. The details and discussion on co n tinuity of eigenv alues can also be obtained from Section 5.2.3 of [30]. An intuitive argument c an b e found in [31]. The roo ts o f a polynom ial equation are shown to depend continu ously on coefficients of the polynomial [32]. The eigenv alues are in fact r oots o f Ch aracteristic P olynomial of a matrix M and given as h ( λ ) = det { λI − M } , a mon ic polyno mial. Further, th e co efficients of h ( λ ) can be expressed in terms of the sum of prin cipal minors of M . Each of th ese principle minor s dep end on the coefficient m ij of matrix M establishing th e co ntinuou s depen dency of r o ots of h ( λ ) o n entries m ij of matr ix M . Theref ore, roots of Cha racteristic P olynomial , i.e., eige n values are a c o ntinuo u s fu nction of entries o f the matrix . As th e maximu m oper ator over a set of continu ous values is continuo u s, the critical eigenvalue of the matrix J r is continu ous with variations in matrix entries. It is essential to n ote the following. I n the p resent work, we are interested in explor in g the neighbo rhood su bspace of an operatin g point. This subspace can be taken as contin uous as the ab rupt state ch anges are not u nder co nsideration o f this work. In various power system ope rations and co nditions, the state-space canno t be co nsidered as contin uous. Th ese situations arr i ve main ly du r ing discrete changes such as co n- tingencies or instan c es of con tr oller saturation . Identification of continuo us subspace is cr itical in such cases bef ore applyin g the propo sed method. 3 The red uced Jacobian is a non-lin e a r functio n o f states J r ( x ) , and so is critical eige nvalue λ c ( x ) . Due to inv olvement of the matr ix inverse, obtainin g analytical expression for λ c ( x ) is n ot possible. Fur ther, in the situatio n where th e state vector x is uncertain and can vary in a state subsp a ce X , the λ c ( x ) will be unce r tain and mo re difficult to estimate. Ther e fore, we propo se this GP learnin g based m ethod to learn the be h avior of λ c ( x ) with x ∈ X . B. Ga ussian Pr ocess Regr ession In the Bayesian optim ization paradig m, wh ere the exact objective fun ction is not known, the Ga ussian process is used extensi vely as a modeling tool. The GP is dev eloped as an extension of multi-variate Gaussian and can b e consider ed as a distribution over random functions [33]. Th e Gau ssian p rocess is a no n-par a m etric metho d, hen ce suitable as a modeling method for pr obability distributions over functio ns. In power systems, the use o f GP h as been restricted to the fo r ecasting applications for power and load. Th e GP ha s been app lied fo r wind p ower fore c ast [34]–[36], solar power forecast [37], an d electricity de m and forecast [3 8]–[40]. Other then these forecasting works, the id ea of u sing GP to learn the dynamics and stability index behavior has not b een explored by the p ower system research community . First, we de fine a gene ral framework for GP regre ssion . Let, a tra ining d ata set D = { x ( i ) , ˆ λ c ( x ( i ) ) } m i =1 where ˆ λ c ( x ( i ) ) is the ob served function value for input vector x ( i ) ∈ R n at th e i th step. Then the GP regression mod el can be g iv en as [33]: ˆ λ c ( x ( i ) ) = λ c ( x ( i ) ) + ǫ ( i ) , i = 1 . . . m (3) Here, ǫ ( i ) are in depend ent and iden tically distributed noise variable with zero-mean, σ n standard de viation no rmal dis- tribution  N (0 , σ 2 n )  . Interested read er can look in to [33] for details of GP f undam e ntals. In this work, the un k nown critical eigenv alue fun c- tion is approximated b y GP regression. The covariance or kernel fu nction k ( x , x ′ ) , bring s our under standing and assumptions abou t the function into GP . A sample set ˆ λ c m = [ ˆ λ c ( x (1) ) , . . . , ˆ λ c ( x ( m ) )] T at o p erating points D m = { x (1) , . . . , x ( m ) } with Gaussian n oise ǫ , and th e analy tic formu la set can be obtain ed for posterio r distribution co rre- sponding to (3). µ m ( x ) = k m ( x ) T ( K m + σ 2 n I ) − 1 q m (4a) k m ( x , x ′ ) = k ( x , x ′ ) − k m ( x ) T ( K m + σ 2 n I ) − 1 k m ( x ′ ) (4b) σ 2 m ( x ) = k m ( x , x ) (4c) Here, µ m ( x ) is mean , k m ( x , x ′ ) is covariance a n d variance is indica ted by σ 2 m ( x ) . The x and x ′ are two sample o perating points from set D m . The k m ( x ) = [ k ( x (1) , x ) , . . . , x ( m ) , x )] and K m = [ k ( x , x ′ )] . In th is work, we infuse our prio r knowledge abo ut critical eigenv alue fun ction as the squ ared exponential cov ariance func tion with zer o mean and un it characteristic leng th . C. GP-UCB The Gau ssian pr ocess upp er confiden c e b ound (GP-UCB) algorithm is an intuitiv e Bayesian metho d [41] f or sampling. For a given δ ∈ (0 , 1) , in a state subspace X , ou r g oal is to learn the mean µ ( x ) for critical eigenv alue function λ c ( x ) with least standard d eviation σ ( x ) and confid e nce le vel 1 − δ . A combin ed strategy to strike th e balance between explor a tio n and exploitatio n can be used for sampling th e p o int x ( i ) . With β i taken independ ent of state vector the sampling strategy will be: x ( i ) = arg max x ( i ) ∈X n µ i ( x ) + p β i +1 σ i ( x ) o (5) Here, ( 5) sug gests that x i is selected where µ i ( x ) + β 1 / 2 i +1 σ i ( x ) maximizes. The µ i ( x ) contributes to enlarging the lev el set of critical eigen value while σ i ( x ) help s in minimizin g the u n certainty . Interested read er can obtain the details of this sampling strategy and GP-UCB fro m [41]. I I I . M A I N R E S U LT : P R S F R A M E W O R K In this sectio n, we will pr esent the main resu lt for devel- oping the PRS fram ew ork. Prior to th is result, we note tha t various work s ha ve been proposed under the generic name pr obabilistic r obust contr ol to perfo rm ro bust control under stochastic uncertainties [42]. The sy stem p erform ance is said to be robustly satisfied (with a fixed probab ility) if a guar antee is provid e d ag ainst almost all of th e po ssible uncertain ties [43]. A similar analog y in th e co ntext o f p ower system small-signal stability can be u sed here for th e PRS framework. Giv en δ ∈ (0 , 1) , if all the points in a su bspace are stable with the prob ability of at least 1 − δ , th en that subspa c e is called PRS subspace. For developing the basis of PRS certificate, we first present a result o n r egret boun d , i. e . , error bound , in Theorem 1 for GP- UCB represented in ( 5). Theorem 1. Let λ c ( x ) be the critical eigen value functio n in state x with the no ise ǫ bou n ded by σ n . Then the follo win g holds with the pr oba bility of at least 1 − δ with δ ∈ (0 , 1) | λ c ( x ) − µ m ( x ) | ≤ p β m +1 σ m ( x ) , ∀ x ∈ X . (6) Her e X is the sample space wherein the states x lie, and m is the number of samplin g po ints. β m = 2 k V k 2 k + 300 γ m ln 3 ( m/δ ) defined in [26]. Pr o of. The results follow dir e ctly from Theorem 6 in [44] and Theorem 1 in [26]. Most impo rtantly , f ollowing T heorem 1, we o btain a pro b - abilistic robust stability certificate as th e following. Theorem 2. ( PRS certificate) F or a given δ ∈ (0 , 1) , the linearized DAE system (2) is pr oba bilistic r o bust small-signa l stable (PRS) in a state subspa ce X with p r o bability 1 − δ if p m ( x , δ ) < 0 (7) wher e p m ( x , δ ) = max x ∈X  µ m ( x ) + p β m +1 σ m ( x )  with m sampling poin ts. Pr o of. The d efinition of PRS cer tificate p m ( x , δ ) < 0 can be obtained directly from T heorem 1. 4 The PRS cer tificate p m ( x , δ ) < 0 verifies that with the probab ility 1 − δ any state vector x in the given state subspac e X is small-signal stab le . This idea does not q uantify the number of points that are stable in a sub space but gives the probab ility with th at a ny inner point is stab le in the consid e r ed subspace. Based u pon Th e o rem 2, we continue th e PRS framework, which covers v arious con ditions as: 1) If the PRS certificate ho lds or p m ( x , δ ) < 0 : The state sub sp ace X is pro babilistic ro bust stable with the con fidence le vel of 1 − δ . 2) Otherwise : a) Subs pace-ba sed Search: Find X ′ ( X su c h that p m ( x , δ ) < 0 , ∀ x ∈ X ′ . b) Confidence-based Search: Find δ ′ > δ such that p m ( x , δ ′ ) < 0 , ∀ x ∈ X . The ratio n ale behin d the seco nd step of th e PRS fram ew ork is that in case the PRS certification can not be established f o r a gi ven set X , one can proceed with two alternati ve searches. Subspa c e-based Sea r ch results into a subspace X ′ , of state space, which is PRS with probability 1 − δ . This X ′ is also a subspace of the original state subspace X . Th e other a lter native Confiden ce-based S ear ch results in a lower confiden ce level 1 − δ ′ , by wh ich the or iginal subspace X is PRS. Therefore , with the PRS certificate and two alternative searches, the PRS framework is complete as it cater s to all po ssible scena r ios. Remark. The e xistence of X ′ in Subspace-ba sed Sea r ch is guaranteed if ∃ x ⋆ ∈ X such that ma x i ( ℜ ( λ i ( J r ( x ⋆ )))) < 0 . In many ca ses, x ⋆ is selecte d as the base stable operating point ar ound that the subspa ce X is co nstructed. Ther efor e, we can at least find X ′ containin g th e stable base poin t. For Con fidenc e - based Searc h , by in creasing δ to δ ′ , we obtain a lower confiden ce lev el 1 − δ ′ . Howe ver , we do not guaran tee th at the new c o nfidence level 1 − δ ′ can be large enoug h for any practical me aning. W ithin our PRS fram ew ork, it is importan t to determine when the GP-UCB search can be terminated or the search is completed . As the GP-UCB attempts to minimize th e lear ning uncertainty related to σ m ( x ) in (5), th e ch ange in values o f σ m ( x ) indicates the level of learn ing err ors with in creasing number of training samp les m . Fur ther, this lear ning level will also be reflec te d in variation o f p m ( x , δ ) while incre a sing m . Therefor e, either of the two indicato rs, i.e., σ m ( x ) and p m ( x , δ ) , c an be used to decide whether the GP-UCB search for learning λ c ( x ) comp letes. Mo reover , we prop ose to use the real par t of critical eige nvalue as a stab ility index instead of minimum dampin g ratio (MDR) in this work. Although, the pro posed framework can b e extended to accomm odate any stability index such as MDR. I V . S I M U L AT I O N S A N D D I S C U S S I O N S In this work, we have used IEEE three- machine nin e-bus system [ 4 5] for tes ting and validation of the pro posed PRS framework. The system h as three P V buses with conventional generato r s at bus num b er 1 , 2 , and 3 . The gene rator 1 is considered as slack bus and has the highest inertia co nstant. 0.96 1 1.04 1.08 |V 3 | -0.04 -0.02 0 0.02 0.04 c (|V 3 |) 95% Confidence region m (|V 3 |) Training Data Fig. 1. An illustrat ion on GP-UCB worki ng on 1-D state subspace of | V 3 | . The blue line denotes mean v alues of the unkno wn function λ c ( | V 3 | ) , and the gray shade describes the 95% confidence interv al. The three red triangles represent three dif ferent sampling points. 1.0 1.2 1.4 1.6 1.8 2.0 P g1 1.0 1.5 2 Q g1 PRS Subspace Base Point Fig. 2. PRS certifica te testi ng for X 1 = { 0 . 909 ≤ P g 1 ≤ 2 . 119; 0 . 852 ≤ Q g 1 ≤ 1 . 988 } indicat ed with black dashed rectangl e W e co n struct the reduced system Jacobian ( J r ) a n d test PRS certificate and fra m ew ork for various P V bus sub spaces. As the PRS framework, in Section I I I, is developed for a gen eral n-dimension al sp a ce vector, the extensio n to the load bus subspace is straig ht fo rward. W e use | V i | to indicate the node voltage magnitu de at i th bus, while P g i and Q g i to indicate real and reactive power o utput of g enerator con n ected to i th bus. All nu merical values are in per un it ( pu ). A general set describing a sub space is ind icated with X while a set for which the PRS certificate (T heorem 2) can b e satisfied is ind icated with X c . T he values of the small-sig n al stable base point are giv en in App endix T able II. First, we discuss the GP learnin g process in on e - d imen- sional voltage magnitude | V 3 | state subspace. Figure 1 shows the learnin g mechanism for λ c as a functio n of | V 3 | with three sampling po in ts shown as red triangles ( m = 3 ). The me a n µ m ( | V 3 | ) is sho wn in blue an d ± 2 σ m ( | V 3 | ) is indicated with gray ar ea covering 95% confide n ce interval. It is clear that for | V 3 | ≥ 1 . 04 , th e un certainty is hig her, leading to a larger gr ay area for achie ving the 9 5% co n fidence level. Therefo r e, it can be co ncluded that w ith a sufficiently large n umber o f sampling points, the unc ertainty in learning λ c ( x ) can be d ecreased to an acceptable le vel. Furth er , the objective of the PRS framew ork is to c e rtify a neigh borho od state subsp a ce of a stable base 5 2.6 3.0 3.4 P g2 0.95 1.15 1.35 Q g2 PRS Subspace Base Point Fig. 3. PRS certifica te testing for X 2 = { 2 . 338 ≤ P g 2 ≤ 3 . 503; 0 . 895 ≤ Q g 2 ≤ 1 . 341 } indicat ed with black dashed rectangle point. It is not task ed to cer tify or ev aluate the stability in the entire state space. This allo ws minimizing the uncertainty values e fficiently with less n umber of samp lin g points. For p ower sy stem stability , th e state space with real and reactive power variables is crucial as they reflect th e ge nerator dynamics impact directly . Further , d ifferent gen e rators have a d ifferent level o f inertia, thus lead ing to a d ifferent stable subspace around the ba se p o int. T o study PRS with real and reactive p ower variations, the PRS certificate verification is done for each g enerator f o r different subspace. Th e Fig u res 2, 3 and 4 shows results for PRS testing in p ower sub space correspo n ding to gen erators 1 , 2 and 3 respec tively . The b lack rectangle indica tes the re gion covered by different sub sp ace set X 1 , X 2 and X 3 while the region covered by blue line is subspace for which µ m ( x ) + 2 σ m ( x ) ≤ 0 . All these figures highligh t the difference in the PRS subspace o f different generato r s. The higher area covered by th e blu e line in Figure 2 is indicative of the fact th at generator 1 can handle higher variations in power befor e leading into in stability . As p m ( x , δ ) ≮ 0 for sets X 1 , X 2 and X 3 , we can perfor m Subspa c e-based S ear ch to find a PRS subspace. The Figure 5 shows two d ifferent subspace where PRS certificate validation is perform e d in P g 1 − | V 1 | subspace. Further, it shows a subspace in blue dash line for which p m ( x , δ ) < 0 . Th erefor e , every state point inside the b lue rectangle h as a 95% chance of being stable. This PRS subspac e is ob ta in ed u sing Subsp a ce-based Se a r ch . Here, we use a subspace description with rectan g ular g e ometry for conv eying the main idea, although any set descrip tion can be used with the propo sed PRS framework. Figure 5 also suggests that a maximum PRS subspace can also be f ound, which will hav e a lesser r ange in | V 1 | and high er in P g 1 state dimension. T h is result is also indicative of the fact th a t the PRS subspace is non-uniq ue, and there exist multiple such subspaces in state space varying in d escriptions. Therefo re, Subspac e-based Sear ch can be tasked to find the maxim um range in the dimension which is most uncer tain o r exhib it larger variations than other states. This paper is focused on presenting the PRS framework id ea; the r efore, such customization is no t in scope. Similar to the result shown in Figure 5, the T able I co ntains 1.1 1.3 1.5 1.7 1.9 P g3 0.40 0.60 0.80 Q g3 PRS Subspace Base Point Fig. 4. PRS certifica te testi ng for X 3 = { 1 . 143 ≤ P g 3 ≤ 1 . 914; 0 . 354 ≤ Q g 3 ≤ 0 . 824 } indicat ed with black dashed rectangl e 1 1.5 2 P g1 1.03 1.038 1.046 |V 1 | 0.80 P g1 2.20 P g1 =1.51, |V 1 |=1.040 1.00 P g1 1.94 1.030 |V 1 | 1.049 Fig. 5. P g 1 − | V 1 | subspace X such that p m ( x , δ ) ≮ 0 ∀ x ∈ X with red line and subspace X c ( X such that p m ( x , δ ) < 0 ∀ x ∈ X c with blue line the dimension s of different X c for different variables. W e observed th a t | V | is the limiting variable and has the least acceptable v ariations in comparison to the power su bspace. Here, the objective is again to show that the PRS framework has b een able to obtain meaning ful X c in d ifferent d imensions and variable subspa c e an d X c is not ind ic a ting largest possible subspace in T able I. The pro posed GP lea rning m ethod can be used to und erstand the behaviour of λ c ( x ) . The Figure 6 depicts variations in µ m ( x ) + p β m +1 σ m ( x ) in | V 2 | − | V 3 | subspace around the base point. From Theorem 1, it is clear that Figu re 6 repr esents upper bo und of critical eigenv alue with probability 1 − δ . This upper bo und, in th is s ubsp a ce, is smooth and an edge n ear | V 2 | = | V 3 | = 1 . 0 5 cro ss the stab ility bound ary . Furthe r, from Theorem 1 it can be concluded that if lea rning un certainty is very less i.e. max { σ m ( x ) } ≪ 1 ∀ x ∈ X , then e ven with higher value of β m +1 the d ifference between actual λ c ( x ) and its estimated mean µ m ( x ) will be very sma ll. The Figure 6 indicate one such situation with max { σ m ( x ) } ∼ 10 − 8 . Therefo re, in this situatio n, the pla n e sh own very clo sely represent λ c ( x ) . T ime Consum p tion Rep ort: All simulatio ns in this work are perfor m ed using GPML [4 6] with MA TL A B 201 8b on a mac hine with Intel X e o n E5- 1630v 4 h aving 3 . 70 GHz clock 6 T ABLE I P R S C E RT I FI E D S U B S PAC E X c T Y P E A N D D I M E N S I O N S I N pu Space T ype V aria bles Minimum Max imum ∆ | V | Space 3-Dimensiona l | V 1 | 1.032 1.048 0.017 | V 2 | 1.017 1.033 0.016 | V 3 | 1.017 1.033 0.016 P g Space 3-Dimensiona l P g 1 1.478 1.554 0.076 P g 2 2.849 2.995 0.146 P g 3 1.486 1.562 0.076 P g − Q g Space 6-Dimensiona l P g 1 1.478 1.554 0.076 P g 2 2.849 2.995 0.146 P g 3 1.486 1.562 0.076 Q g 1 1.385 1.457 0.071 Q g 2 1.092 1.148 0.056 Q g 3 0.576 0.605 0.030 P g − Q g − | V | Space 9-Dimensiona l P g 1 1.478 1.554 0.076 P g 2 2.849 2.995 0.146 P g 3 1.486 1.562 0.076 Q g 1 1.385 1.457 0.071 Q g 2 1.092 1.148 0.056 Q g 3 0.576 0.605 0.030 | V 1 | 1.035 1.045 0.010 | V 2 | 1.020 1.030 0.010 | V 3 | 1.020 1.030 0.010 Fig. 6. Cri tical eigen value ( λ c ) upper bound plane in two-dimensional | V 2 | − | V 3 | subspace with 95% confidence lev el, with max { σ m ( x ) } ∼ 10 − 8 speed and 16.0 GB of RAM. It takes ap p roxima te ly 2 . 3 se c to o btain results of Figure 6 while PRS certificate is verified in appro x imately 2 . 5 sec in Figure 5. Th e time consumption depend s on the initial sampling point and can be im proved using better initialization strategies, w h ich are not in the scope of this work. Furth e r, fo r higher dim ensions and larger data sets, spa r se Gaussian p rocesses [4 7] can b e used for less time consump tion. V . C O N C L U S I O N In this p aper, we pr esent a n ovel proba b ilistic robust small- signal stability (PRS) framework based on the n o n-par ametric Gaussian pro cess to learn the behavior of the cr itical eigen- value of the reduced J acob ia n . I n particular , we use GP to estimate the real- part of the critical eigenv alue by learning its upp er boun d an d thu s assess the system stability . The robustness of system stability is verified in the state subspace with a given probability by establishing a PRS certificate. I n case the PRS cer tificate is not satisfied, we pro pose Subspa ce- based S ear ch and Confid ence-ba sed S ear ch to complete PRS framework. The simulation results on the WSCC n etwork il- lustrates the perfo rmance of the p roposed method in certifying the PRS certificate an d learning the critical eigen values. The pr e sent work op ens u p a ne w d imension in th e ar e a of p robab ilistic stability assessment o f the power system. T he certificate presente d he re ca n be applied to larger test systems and can be furth er op tim ized, fo r example, to be less ti me- consumin g . The learning of the critical eigen value b ehavior can provide a refined unde rstanding of the stability under uncertainty . The r e also exists a potential for solving prob lems to find th e largest PRS su b space f or a g iv en confidence level. Based on this work , one can de sign a general PSSS assessment method that ne ith er req uire th e pd f o f input unc e rtainty nor rely on state-specific appro ximations. The se idea s will b e explored further in upc o ming works. A C K N O W L E D G E M E N T W e than k Prof. Ashu V erm a an d Dr . Zha i Chao for their suggestions. A P P E N D I X T ABLE II S M A L L - S I G N A L S TA B L E B A S E P O I N T Bus Number 1 2 3 Real Power ( P g ) 1.515 2.922 1.523 React iv e Power ( Q g ) 1.421 1.119 0.590 V oltage Magnitude ( | V | ) 1.040 1.025 1.025 R E F E R E N C E S [1] P . Kundur , N. J. Balu, and M. G. Lauby , P ower system stability and contr ol . McGraw- hill New Y ork, 1994, vol. 7. [2] J. Macho wski, J. Bialek, and J. Bumby , P ower system dynami cs: stability and contr ol . J ohn Wile y & Sons, 2011. [3] H. D. Nguyen, T . L. V u, J.-J. Slotine , and K. Turi tsyn, “Contra ction analysi s of nonlinear dae systems, ” arX iv preprint , 2017. [4] F .-E. W olter and B. Be rger , “Dif ferent ial geometric foundations for po wer flow computations, ” arXiv pr eprint arXiv:1903.11131 , 2019. [5] T . L. V u and K. T uritsyn, “L yapunov functions family approach to transien t stabili ty assessment, ” IEEE T ransac tions on P ower Systems , vol. 31, no. 2, pp. 1269–1277, March 2016. [6] Z . Zhang, H. D. Nguyen, K. Tu ritsyn, and L. Daniel, “Probabilistic po wer flow computa tion via lo w-rank and sparse tensor recov ery , ” arXiv pre print arXiv:1508.02489 , 2015. [7] D. Mehta, H. D. Nguyen, and K. T uritsyn, “Numerical polynomial homotop y continuation m ethod to locate all the power flow solution s, ” IET Gener ation, T ransmission & Distribution , vol. 10, no. 12, pp. 2972– 2980, 2016. [8] D. K. Molzahn, V . Dawa r , B. C. Lesieutre, and C. L. DeMarco, “Suf ficient conditions for powe r flow insolv ability considering reacti ve po wer limited generators with applica tions to volt age stability margins, ” 2013 IREP Symposium , pp. 1–11, 2013. [9] L . A olaritei, S. Bolognani , and F . D ¨ orfler , “ A distrib uted voltage stab ility margin for powe r distri bution net works, ” IF AC-P apersOnLine , vol. 50, no. 1, pp. 13 240–13 245, 2017. [10] M. Ali, A. Dymarsky , and K. Turit syn, “Transv ersality enforced newto n raphson algorithm for fast calcula tion of maximum loadabilit y , ” IE T Genera tion, T ransmi ssion & Distributi on , 2017. 7 [11] R. C. Burchett and G. Heydt, “Probabilisti c methods for power system dynamic stabil ity s tudies, ” IEEE Tr ansactions on P ower Apparat us and Systems , no. 3, pp. 695–702, 1978. [12] C. Pans, Z. Dong, P . Zhang, and X. Y in, “Probabilist ic analysis of po wer system small signal stability region, ” in 2005 International Confer ence on Contr ol and Automati on , vol. 1. IEEE, 2005, pp. 503–509. [13] R. Arrieta, M. A. R´ ıos, and A. T orres, “Contin gency ana lysis and risk assessment of sm all signal insta bility , ” in 2007 IEEE Lausanne P ower T ech . IEEE, 2007, pp. 1741–1746. [14] R. Preece, K. Huang, and J. V . Milano vi ´ c, “Probab ilistic small- disturban ce stabil ity assessment of uncertain power systems using ef ficient estimation m ethods, ” IEEE T ransact ions on P ower Systems , vol. 29, no. 5, pp. 2509–2517, 2014. [15] Z. Xu, Z . Dong, and P . Z hang, “Probabil istic small s ignal analysis using monte carlo simulation, ” in IEEE PESGM, 2005 . IEEE, 2005, pp. 1658–1664. [16] J. L. Rueda, D. G. Colom ´ e, and I. Erlich, “ Assessment and enhance ment of small signal stabilit y consid ering uncer tainties, ” IEEE Tr ansactions on power systems , vol. 24, no. 1, pp. 198–207, 2009. [17] R. Preece and J. V . Milano vi ´ c, “Risk-ba sed small-disturba nce securit y assessment of po wer systems, ” IEEE T ransactio ns on P ower Delivery , vol. 30, no. 2, pp. 590–598, 2014. [18] H. Huang, C. Chung, K. W . Chan, and H. Chen, “Quasi-monte carlo based probabil istic small signal stabi lity analysis for po wer systems with plug-in electric vehicl e and wind power inte gration, ” IEEE T ransac tions on P ower Systems , vol. 28, no. 3, pp. 3335–3343, 2013. [19] R. Preece, N. C. W oolle y , and J. V . Mi lanovi´ c, “The pro babilisti c colloc ation method for po wer-syste m damping and voltage c ollapse studies in the presenc e of uncertainti es, ” IEEE T ransactions on P ower Systems , vol. 28, no. 3, pp. 2253–2262, 2012. [20] J. Rueda and I. E rlich, “Probabil istic framew ork for risk analysis of po wer system s mall-signa l stabili ty , ” Proc eedings of the Instit ution of Mec hanical Engineer s, P art O: J ournal of Risk and Reliability , vol. 226, no. 1, pp. 118–133, 2012. [21] K. W ang, C. Chung, C. T se, and K. T sang, “Improv ed probabili stic method for power system dynamic stabili ty studies, ” IET GTD , vol. 147, no. 1, pp. 37–43, 2000. [22] J. Xu, P . K. Kanying i, K. W ang, G. Li, B. Han, and X . Jiang, “Probabi listic small s ignal stabil ity analysis with lar ge scale inte gration of wind power consideri ng dependence, ” R enew able and Sustainabl e Ener gy R eview s , vol. 69, pp. 1258–1270, 2017. [23] Z. W ang, C. Shen, and F . Liu, “Probabili stic analysis of small signal stabili ty for power systems with high penetra tion of wind generation , ” IEEE T ransacti ons on Sustainable Energy , vol. 7, no. 3, pp. 1182–1193, 2016. [24] R. Krishan, A. V erma, S. Mishra, and P . R. Bijwe, “ Analysis of hopf bifurca tion with forecast uncerta inties in load/ge neration, ” IET GTD , vol. 11, no. 6, pp. 1531–1538, 2017. [25] K. N. Hasan, R. Preece, and J. V . Milano vi ´ c, “Existing approaches and trends in uncert ainty modelling and probabilist ic stability analysis of po wer systems with rene wabl e generation, ” Renewa ble and Sustainable Ener gy R eview s , vol. 101, pp. 168–180, 2019. [26] C. Zhai and H. D. Nguyen, “Re gion of attra ction for power systems us- ing gaussian process and con verse lyapuno v function –part i: Theoretic al frame work and off -line study , ” arXiv prepri nt arXiv:1906.03590 , 2019. [27] P . Pareek, K. T uritsyn, K. Dvijotham, and H. D. Nguyen, “ A suf ficient conditi on for small-signa l stabi lity and con struction of robu st stabi lity regi on, ” in IEEE PESGM, 2019 . IEE E, 2019. [28] D. Y ang and V . Ajjarapu, “Critical eigen v alues tracing for po wer system analysi s via continuati on of in vari ant subspaces and projected arnoldi method, ” in 2007 IEEE PE SGM . IEEE, 2007, pp. 1–9. [29] C. Luo and V . Ajjara pu, “In va riant subspace based eigen v alue tracin g for po wer system small-signal stabilit y analysis, ” in 2009 IE E E PESGM . IEEE, 2009, pp. 1–9. [30] D. Serre, Matrices Theory and Applications- Second Edition . Springer Science and Media, 2010. [31] C. D. Meyer , Matrix analysis and applied linear algebra . Siam, 2000, vol. 71. [32] G. Harris and C. Martin, “Shorter notes: The roots of a polynomia l v ary continuo usly as a functio n of the coef ficients, ” P r oceeding s of the American Mathematic al Soci ety , pp. 390–392, 1987. [33] C. K. Will iams and C. E. Rasmussen, Gaussian pro cesses for machine learning . MIT press Cambridge, MA, 2006, vol. 2, no. 3. [34] D. Lee and R. Baldick, “Short-term wind power ensemble predicti on based on gaussia n proce sses and neural networks, ” IEEE T ransacti ons on Smart Grid , vol. 5, no. 1, pp. 501–510, 2013. [35] N. Che n, Z. Q ian, I. T . N abney , and X . Meng, “W ind power forecasts using gaussian processes and numerical weather predicti on, ” IE EE T ransac tions on P ower Systems , vol. 29, no. 2, pp. 656–665, 2013. [36] J. Y an, K. Li, E.-W . Bai, J. Deng, and A. M. Fol ey , “Hybrid probabi listic wind po wer forecast ing using temporally local gaussian process, ” IEEE T ransac tions on Sustainable Energy , vol . 7, no. 1, pp. 87–95, 2015. [37] H. Sheng, J. Xiao, Y . Cheng, Q. Ni, and S. W ang, “Short-t erm solar po wer forecasti ng based on weighted gaussian process re gression, ” IEE E T ransac tions on Industrial Electr onics , vol . 65, no. 1, pp. 300–308, 2017. [38] M. Blum and M. Riedmill er , “Electricity demand forecast ing using gaus- sian processes, ” in W orkshops at the T wenty -Seve nth A AAI Confer ence on Artificial Intellig ence , 2013. [39] Y . Y ang, S. Li, W . L i, and M. Qu, “Power load probability density forecast ing using gaussian process quantile regressio n, ” A pplied Energy , vol. 213, pp. 499–509, 2018. [40] D. W . v an der Meer , M. Shepero, A. Svensson, J. W id ´ en, and J. Munkha mmar, “Proba bilistic forecasting of electri city consumptio n, photo voltai c powe r generati on and net demand of an indi vidual building using gaussian processes, ” Applied energy , vol. 213, pp. 195–207, 2018. [41] N. Sriniv as, A. Krause, S. M. Kakade, and M. Seeger , “Gaussian process optimiza tion in the bandit s etting: No regre t and experi mental desi gn, ” arXiv prepri nt arXiv:0912.3995 , 2009. [42] G. C. Calafio re and M. C. Campi, “The scenario approach to robust control design, ” IEEE T ransacti ons on Automatic Cont r ol , vol. 51, no. 5, pp. 742–753, 2006. [43] G. C. Cala fiore and F . Dabbene , “Probabilisti c rob ust control, ” in 2007 American Contr ol Confer ence . IEEE, 2007, pp. 147–158. [44] N. Srini vas, A. Krause, S. M. Kakade, and M. W . Seege r , “Information- theoret ic re gret bounds for gaussian process optimiz ation in the bandit setting , ” IEEE T ransacti ons on Information Theory , vol. 58, no. 5, pp. 3250–3265, 2012. [45] M. Pai, Energy funct ion analysis for power system stability . Springer Science & Business Media, 2012. [46] C. E. Rasmussen and H. Nickisch, “Gaussian processes for machine learni ng (gpml) toolbox, ” Journal of machine learning resear ch , v ol. 11, no. Nov , pp. 3011–3015, 2010. [47] E. Snelson and Z. Ghahrama ni, “Sparse gaussian processes using pseudo-inp uts, ” in Advance s in neur al information proce ssing systems , 2006, pp. 1257–1264.