Risk of Phase Incoherence in Noisy Power Networks With Delayed Feedback Control

7th IFAC Workshop on Distributed Estimation and 7th IFAC Workshop on Distributed Estimation and 7th on Distributed and Control Networked 7th IFAC IFACinWorkshop Workshop onSystems Distributed Estimation Estimation and 7th IFAC Workshop on Distributed Estimation Availableand online at www.sciencedirect.com Control in Networked Systems Control in Networked Systems 7th IFAC Workshop on Distributed Estimation and Groningen, NL, August 27-28, 2018 Control in Networked Systems Control in Networked Systems Groningen, NL, August 27-28, 2018 Groningen, NL, August 27-28, 2018 Control in Networked Systems Groningen, NL, NL, August August 27-28, 27-28, 2018 2018 Groningen, Groningen, NL, August 27-28, 2018

ScienceDirect

IFAC PapersOnLine 51-23 (2018) 142–147

Risk of Phase Incoherence in Noisy Power Risk of Phase Incoherence in Noisy Power Risk of Phase Incoherence in Noisy Power Risk of Phase Incoherence in Noisy Power Networks With Delayed Feedback Control Risk of Phase Incoherence in Noisy Power Networks With Delayed Feedback Control Networks With Delayed Feedback Control Networks With Delayed Feedback Control Networks With Delayed Feedback Control ∗ C. Somarakis ∗∗ Y. Ghaedsharaf ∗∗ F. D¨ orfler ∗∗ ∗∗ N. Motee ∗

∗ Y. Ghaedsharaf ∗ ∗ F. D¨ ∗∗ N. Motee ∗ ∗ Somarakis o rfler ∗∗ ∗ ∗ ∗ ∗∗ N. Motee ∗ ∗ Somarakis o Somarakis ∗∗ Y. Y. Ghaedsharaf Ghaedsharaf ∗∗ F. F. D¨ D¨ orfler rfler ∗∗ ∗∗ N. Motee ∗ Somarakis Y. Ghaedsharaf F. D¨ o rfler N. Motee ∗ ∗ Department of Mechanical Engineering and Mechanics, of Mechanical Engineering and Mechanics, ∗ ∗ Department Department of Engineering and ∗ Department of Mechanical Mechanical Engineering and Mechanics, Mechanics, Lehigh University, Bethlehem, PA, 18015, USA ∗ Lehigh University, Bethlehem, PA, 18015, USA Department of Mechanical Engineering and Mechanics, Lehigh University, Bethlehem, PA, USA Lehigh{csomarak,ghaedsharaf,motee}@ University, Bethlehem, PA, 18015, 18015, USA (e-mail: lehigh.edu). (e-mail: {csomarak,ghaedsharaf,motee}@ lehigh.edu). ∗∗ Lehigh University, Bethlehem, PA, 18015, USA (e-mail: {csomarak,ghaedsharaf,motee}@ lehigh.edu). Ctrl Lab., ETH, Z¨ u rich, SUI (e-mail: [email protected]) ∗∗ Automatic (e-mail: {csomarak,ghaedsharaf,motee}@ lehigh.edu). ∗∗ Ctrl Lab., ETH, Z¨ u rich, SUI (e-mail: [email protected]) ∗∗ Automatic (e-mail: {csomarak,ghaedsharaf,motee}@ lehigh.edu). Automatic Ctrl Lab., ETH, Z¨ u rich, SUI (e-mail: [email protected]) ∗∗ urich, SUI (e-mail: [email protected]) ∗∗ Automatic Ctrl Lab., ETH, Z¨ Automatic Ctrl Lab., ETH, Z¨ urich, SUI (e-mail: [email protected]) Abstract: A systemic risk framework is proposed for robustness analysis in linear post-fault Abstract: A systemic risk framework is proposed for robustness analysis in linear post-fault Abstract: risk framework is proposed for analysis in post-fault dynamics ofA power grid. We investigate scenarios where perturbations, induced by exogenous Abstract: Aaa systemic systemic riskWe framework is scenarios proposed where for robustness robustness analysis in linear linear post-fault dynamics of power grid. investigate perturbations, induced by exogenous Abstract: A systemic risk framework is proposed for robustness analysis in linear dynamics of a power grid. We investigate scenarios where perturbations, induced by exogenous dynamics of a power grid. We investigate scenarios where perturbations, induced by exogenous sources of of disturbances, disturbances, steer steer the the power power network network away away from from operation operation equilibria. equilibria. In In post-fault addition, sources addition, dynamics ofdisturbances, a power We investigate scenariosaway where perturbations, induced byIn exogenous sources the network from operation equilibria. addition, sources ofcontrol disturbances, steer the power power network away from operation equilibria. Inas addition, feedbackof control laws grid. thatsteer are prone prone to communication communication time-delays are implemented implemented as means feedback laws that are to time-delays are aaa means sources ofcontrol disturbances, steer theperturbations. power networkSystemic awaytime-delays from equilibria. Inas addition, feedback laws that are prone to communication are implemented means to stabilize the network against riskoperation methods are then employed to feedback control laws that are prone to communication time-delays are implemented as a means to stabilize the network against perturbations. Systemic risk methods are then employed to feedback laws that arenetwork prone toto communication time-delays are implemented as a due means to stabilize the network against perturbations. risk methods are employed to assess thecontrol possibility of the deviateSystemic into undesirable modes of then operation to to stabilize the network against perturbations. Systemic risk methods are then employed assess the possibility of the network to deviate into modes of operation due to stabilize the network against perturbations. risktime-delay methods are employed assess the possibility the network to deviate into undesirable modes of operation due to assess theadditive possibility ofand the network tonegative deviateSystemic into undesirable undesirable modesimplemented of then operation due to to both the additive noiseof and the possibly negative effect of the the implemented feedback both the noise the possibly effect of time-delay feedback assess the possibility of the network to deviate into undesirable modes of operation due to both the additive noise and the possibly negative effect of the time-delay implemented feedback both the additive noise and the possibly negative effect of the time-delay implemented feedback control. Numerical explorations reveal possible limitations regarding the efficiency of the closedcontrol. explorations reveal possible limitations regarding the efficiency of the closedboth theNumerical additive and the possibly negative effect ofgains, the time-delay implemented feedback control. Numerical explorations reveal possible limitations regarding the of closedloop system, as a noise result of network connectivity, control time-delay, and noise. control. Numerical explorations reveal possible limitations regarding the efficiency efficiency of the the closedloop system, as a result of network connectivity, control gains, time-delay, and noise. control. Numerical explorations reveal possible limitations regarding the efficiency of the loop system, as a result of network connectivity, control gains, time-delay, and noise. loop system, as a result of network connectivity, control gains, time-delay, and noise. closed© 2018, IFAC as (International of Automatic Control) by Elsevier Ltd. rights reserved. loop system, a result ofFederation network connectivity, control Hosting gains, time-delay, andAll noise. 1. INTRODUCTION INTRODUCTION induce spurious dynamics that invalidate the model itself. 1. induce spurious dynamics that invalidate the model itself. 1. INTRODUCTION induce spurious dynamics that invalidate the model itself. This is because fluctuations of dynamics may steer the 1. INTRODUCTION induce spurious dynamics that invalidate the model itself. This is because fluctuations of dynamics may steer the 1. INTRODUCTION induce spurious dynamics that invalidate the model itself. This is because fluctuations of dynamics may steer the solution outside the basin of attraction of the nonlinear is because fluctuations of dynamics may steer the Power grids are constantly supervised by system operators This solution outside the basin of attraction of the nonlinear Power grids are constantly supervised by system operators This is because fluctuations of dynamics may steer the solution outside the basin of attraction of the nonlinear Power grids are constantly supervised by system operators solution outside the basin of attraction of the nonlinear system. Secondly, Secondly, mainstream mainstream network network synthesis synthesis methods methods to prevent or constantly confront potential faults or unexpected Power grids are supervised by system operators system. to prevent or confront potential faults or unexpected solution outside the the basin of attraction of theadmittance nonlinear system. Secondly, mainstream network synthesis methods Power grids are constantly supervised by system operators to prevent or confront potential faults or unexpected system. Secondly, mainstream network synthesis methods regard adjusting power voltage of the grid to prevent or confront potential faults or unexpected disturbances. Mainstream emergency indicators regard regard adjusting the power voltage of the grid admittance disturbances. emergency indicators regard system. Secondly, mainstream network synthesis methods adjusting power voltage of to preventnadir or Mainstream confront potential faults or of unexpected disturbances. Mainstream emergency indicators regard parameters. Thisthe is in reality a costly, ifgrid not admittance unfeasible, disturbances. Mainstream emergency indicators regard regard adjusting the power voltage of the the grid admittance frequency phenomena or maximum maximum rate change of regard parameters. This is in reality aa costly, if not unfeasible, frequency nadirMainstream phenomena or rate of change of regard adjusting the power voltage of the grid admittance parameters. This is in reality costly, if not unfeasible, disturbances. emergency indicators regard frequency nadir phenomena or maximum rate of change of solution as it would verbal of of parameters. This is in require reality aa costly, ifreadjustment not unfeasible, frequency. nadir Whenphenomena these indicators exceed pre-defined safety frequency or maximum rate of change solution as it would require verbal readjustment of frequency. When these indicators exceed pre-defined safety parameters. This is inlines. reality aa costly, not unfeasible, solution as would require verbal of frequency nadir phenomena or maximum rate of change of physical frequency. When these indicators exceed pre-defined safety solution as it it would require verbal ifreadjustment readjustment of transmission margins, they typically trigger corrective measures such frequency. When these indicators exceed pre-defined safety physical transmission lines. margins, they typically trigger corrective measures such solution as it would require a verbal readjustment of physical transmission lines. frequency. When these indicators exceed pre-defined safety margins, they typically trigger corrective measures such margins, they typically trigger corrective measures such physical transmission lines. as intervention control policies or disconnection of power Contribution. Motivated by our recent results on systemic as intervention control policies or disconnection of power physical transmission lines. margins, they trigger or corrective measures such Contribution. as intervention control policies disconnection of Motivated by recent results systemic as intervention control policies or disconnection of power power generators fromtypically the grid. grid. Contribution. Motivated by our our recentsystems results on on systemic Contribution. Motivated by our recent results on systemic generators from the measures for networked control (Somarakis as intervention control policies or disconnection of power risk risk measures for networked control systems (Somarakis generators from the grid. Contribution. Motivated by our recent results on systemic generators from the grid. risk measures for networked control systems (Somarakis et al. (2018a,b, 2017)), we tackle the aforementioned risk measures for networked control systems (Somarakis Shortcomings of the the grid. aforementioned methods have moti- et al. (2018a,b, 2017)), we tackle the aforementioned issues generators from issues Shortcomings of the aforementioned methods have motirisk measures for networked control systems (Somarakis et al. (2018a,b, 2017)), we tackle the aforementioned issues Shortcomings aforementioned motion aa two step approach of the noisy swing equations. Shortcomings ofof the aforementioned methods have moti- et al. (2018a,b, 2017)), we tackle thelinear aforementioned issues vated aa number numberof ofthe researchers (such as asmethods Poolla et ethave al. (2017); on two step approach of the noisy linear swing equations. vated researchers (such Poolla al. (2017); et al. (2018a,b, 2017)), we tackle the aforementioned issues on a two step approach of the noisy linear swing equations. Shortcomings of the aforementioned methods have motivated a number of researchers (such as Poolla et al. (2017); on a two step approach of the noisy linear swing equations. first we propose a state feedback law that will attempt Tegling et al. (2015); Paganini(such and Mallada (2017)) to in- At vated a number of researchers as Poolla et al. (2017); At first propose aa state feedback law that will attempt Tegling et al. (2015); Paganini and Mallada (2017)) to inon astabilize twowe step approach ofinthe noisy linear swing equations. At first we propose state feedback law that will attempt vated a number of researchers (such as Poolla et al. (2017); Tegling et al. (2015); Paganini and Mallada (2017)) to inAt first we propose a state feedback law that will attempt to the system, the face of disturbances. The vestigateetnovel performance measures for power networks Tegling al. (2015); Paganini and Mallada (2017)) to in- to stabilize the system, in the face of disturbances. The vestigate novel performance measures for power networks At first we propose a state feedback law that will attempt to stabilize the system, in the face of disturbances. The Tegling et al. (2015); Paganini and Mallada (2017)) to investigate novel performance measures for power networks feedback gains will duplicate the grid topology into a convestigate novel performance measures for power networks to stabilize the system, in the face of disturbances. The leveraging recent advances in distributed multi-agent co- feedback gains will duplicate the grid topology into a conleveraging recent advances in distributed multi-agent coto stabilize the system, in the face of disturbances. The feedback gains will duplicate the grid topology into a convestigate novel performance measures for power networks leveraging recent advances in distributed multi-agent cotrol scheme dictated by a presumable higher-level commuleveraging recent advances in distributed multi-agent cofeedback gains will duplicate the grid topology into a conoperative algorithms (see Fardad et al. (2014); Summers trol scheme dictated by aa presumable higher-level commuoperative algorithms algorithms (see Fardad Fardad et al. al. (2014); (2014); Summers feedback gains will duplicate the grid topology into a control scheme dictated by presumable higher-level commuleveraging recent advances in distributed multi-agent cooperative (see et Summers trol scheme dictated by a presumable higher-level commulayer. Flow of information on this layer may, howet al. (2015); Siami and Fardad Motee et (2016); Bamieh et al. nication operative algorithms al. (2014); Summers nication layer. Flow of on this may, howet al. (2015); Siami (see and Motee (2016); Bamieh et al. trol scheme dictated byinformation a presumable higher-level communication layer. Flow of information on this layer layer may, howoperative (see Fardad et al. (2014); Summers et al. (2015); and Motee (2016); Bamieh et al. nication layer. Flow of information on this layer may, however, occur with time-delay, taken constant and uniform (2012); D¨ oalgorithms rfler Siami and Bullo (2012) and references therein). et al. (2015); Siami and Motee (2016); Bamieh et al. ever, occur with time-delay, taken constant and uniform (2012); D¨ o rfler and Bullo (2012) and references therein). nication layer. Flow of information on this layer may, however, occur with time-delay, taken constant and uniform et al. (2015); and (2012) Moteeand (2016); Bamieh et al. ever, (2012); D¨ oorfler and references therein). over all interconnected generators. We will investigate the (2012); D¨ rfler Siami and Bullo Bullo (2012) and references therein). occur with time-delay, taken constant and uniform all interconnected generators. We will investigate the ever, occur with time-delay, taken constant and uniform In their D¨ majority, these works study the transienttherein). stability over over all interconnected generators. We will investigate the (2012); o rfler and Bullo (2012) and references closed-loop network dynamics and we will evaluate the risk over all interconnected generators. We will investigate the In their majority, these works study the transient stability closed-loop network dynamics and we will evaluate the risk In their majority, these works study the transient stability all interconnected generators. We will investigate the problem on powerthese networks, focusing ontransient post-fault swing over closed-loop network dynamics and we will evaluate the risk In their majority, works study the stability closed-loop network dynamics and we will evaluate the risk of rotor angle incoherence between synchronous generaproblem on power networks, focusing on post-fault swing of rotor angle incoherence between synchronous generaIn their majority, these works study theon transient stability problem on power networks, focusing post-fault swing problem on(see power networks, focusing onor post-fault swing network dynamics and we will evaluate the risk dynamics (see Sauer and Pai (2006) or Varaiya et al. closed-loop of angle incoherence between synchronous generaof rotor rotor angle incoherence between synchronous generators. These events are identified as systemic events. At the dynamics Sauer and Pai (2006) Varaiya et al. tors. These events are identified as systemic events. At the problem on(see power focusing onor swing dynamics and (2006) Varaiya et al. dynamics (see Sauer and Pai Pai (2006) orpost-fault Varaiya et al. same of rotor angle incoherence between synchronous genera(1985); Poolla Poolla etSauer al.networks, (2017); Vu and Turitsyn (2017)). Retors. These events are identified as events. At time, we assess the effectiveness of the proposed (1985); et al. (2017); Vu and Turitsyn (2017)). Retors. These events are on identified as systemic systemic events. At the the same time, we assess on the effectiveness of the proposed dynamics (see Sauer and Pai (2006) or Varaiya et al. (1985); Poolla et al. (2017); Vu and Turitsyn (2017)). Retors. These events are on identified as is systemic events. At the cent approaches to the problem regard multi-agent frametime, we assess the effectiveness of the proposed (1985); Poolla etto al.the (2017); Vu and Turitsyn (2017)). Re- same time-delayed feedback. This work authors’ preliminary same time, we assess on the effectiveness of the proposed cent approaches problem regard multi-agent frameis authors’ preliminary (1985); etto al.the (2017); and Turitsyn (2017)). Re- time-delayed cent problem regard multi-agent framesame time, wefeedback. assess onThis the work effectiveness of the proposed worksapproaches onPoolla the linearized swingVu equations, evaluated around time-delayed feedback. This work is preliminary cent approaches to the problem regard multi-agent frametime-delayed feedback. This work is authors’ authors’ preliminary attempt towards developing a risk-oriented risk-oriented performance works on the linearized swing equations, evaluated around towards developing a performance cent approaches to the Unexpected problem regard multi-agent frame- attempt works on the linearized swing equations, evaluated around works on the linearized swing equations, evaluated around time-delayed feedback. This work is authors’ preliminary a desired equilibrium. disturbances, typically attempt towards developing risk-oriented performance attempt towardsframework developingforaa efficient risk-oriented performance and robustness network synthesis a desired equilibrium. Unexpected disturbances, typically robustness framework for network synthesis works onas the linearized swing equations, around and a desired equilibrium. Unexpected disturbances, typically attempt towards developing a efficient risk-oriented performance modeled additive white noise, excite theevaluated linear dynamics dynamics and robustness framework for efficient network synthesis a desired equilibrium. Unexpected disturbances, typically techniques in power grids. excite the linear and robustness framework for efficient network synthesis modeled as additive white noise, techniques in power grids. a desired equilibrium. Unexpected disturbances, typically modeled as additive white noise, excite the linear dynamics excite the linear dynamics and robustness framework for efficient network synthesis and divert their trajectories into fluctuation around the techniques in power grids. modeled as additive white noise, techniques in power grids. and divert their trajectories into fluctuation around the modeled as additive whiteauthors noise, excite the aspects linear dynamics and their into fluctuation around the and divert their trajectories trajectories into fluctuation around the techniques 2. in power grids. stabledivert equilibrium. The measure of these MATHEMATICAL NOTATION stable equilibrium. The authors measure aspects of 2. MATHEMATICAL NOTATION and divert their into fluctuation around the stable equilibrium. The authors measure aspects of these 2. stable equilibrium. The authors measure aspects of these these fluctuations (suchtrajectories as aggregate variability) utilizing con2. MATHEMATICAL MATHEMATICAL NOTATION NOTATION fluctuations (such as aggregate variability) utilizing constable equilibrium. The authors measure aspects of these fluctuations (such as aggregate variability) utilizing con2. MATHEMATICAL NOTATION cepts such as(such H2 or H norms. Moreover, H /H based n fluctuations as aggregate variability) utilizing con∞ norms. Moreover, H2 /H∞ based By R we understand the n-dimensional Euclidean space or H n cepts such as H 2 ∞ 2 ∞ 2 asHaggregate ∞ norms. variability) 2 ∞ fluctuations concepts as H Moreover, H /H based n By R understand the n-dimensional space 2 ∞ n we network synthesis are investigated in order to T n Euclidean ∞ 2 ∞ By R we understand the n-dimensional Euclidean space cepts such such as(such H222 or oralgorithms H∞ Moreover, Hutilizing based n ∞ norms. 2 /H ∞ with elements y = (y , . . . , y ) ∈ R , represented as network synthesis algorithms are investigated in order to T n By R we understand the n-dimensional Euclidean space 1 n cepts such as H or H norms. Moreover, H /H based n with elements y = (y , . . . , y ) ∈ R , represented as network synthesis algorithms are investigated in order to T n 2 ∞ 2 ∞ 1 n network synthesis algorithms are investigated in order to T n minimize the effect of noise on the linear equations. n×n By R we understand the n-dimensional Euclidean space with elements y = (y , . . . , y ) ∈ R , represented as T n 1 n with elements y = (y , . . . , y ) ∈ R , represented as column vectors. A square M ∈ R is represented by minimize the effect of noise on the linear equations. 1 n n×n 1 n column vectors. A square M ∈ R is represented by network in order to column T n×n n minimize the of the linear minimizesynthesis the effect effectalgorithms of noise noise on onare theinvestigated linear equations. equations. n×n n with elements y = (y , . . . , y ) ∈ R , represented as vectors. A square M R is represented by n×n 1 = n∈ column vectors. A square M ∈ R is represented by n vectors m ∈ R as M [m | m | . . . | m ] = [m n i ∈ Rn as M = [m1 | n×n 2 | . . . | mn ] = [m ij ]] There are two main with H minimize the effect ofshortcomings noise on the linear equations. 2 /H ∞ control.  by n vectors m m i ∈ 1 R 2 |is n ] = [m ij ] n /H control. There are two main shortcomings with H column vectors. A square M ∈ represented n vectors m R as M = [m | m . . . | m 2 ∞ n  th 1 By, 2 ij control. There are two shortcomings with H 2 ij ∞ n vectors miii ii∈th =m [m . . .diag | mnnn{m ] =i } [m  we with m element /Hmodel There are two main main shortcomings withlinear H222 /H Firstly, stochastic perturbation of the may 1 | mM 2 |= ij ] th R ij the j ..1 ∞ control. n as M of ∞ with m of = we th Firstly, stochastic perturbation of the may ij the j .1 By, n vectors mthe Relement as Mmatrix =m [m | with mM . . .diag | mn{m ]elements =ii } [m th i i∈ 2 |diagonal ij ] m the element of m By, M = diag {m } we thdiagonal There two main shortcomings withlinear H2 /Hmodel Firstly, stochastic perturbation of linear model may ij jj .M ∞ control.  understand Firstly,are stochastic perturbation of the the linear model may with ij i with m the i element of m By, M = diag {m } we ij i thdiagonal matrixj M with diagonal elements understand the T T {m with m the i element of m . By, M = diag } we understand the diagonal matrix M with diagonal elements Firstly, stochastic perturbation of the linear model may ij j i  m Also, y (y .. .. ,, yywith and M = [m ]] T = T elements understand diagonal matrix diagonal research is supported by NSF CAREER ECCS-1454022 and 1 ,, .. .. .. ,, m n ..the 1 ,, .. M n)  This m Also, y (y M T = T = [mji This research is supported by NSF CAREER ECCS-1454022 and 1, . . . , m n)  understand diagonal diagonal m Also, y = (y .. M .. .. ,, the yywith ) and and M = [m ]  This research is 1 n ji m ..,m mnnnn ..the Also, yTTT vector =matrix (y1111 ,,and M TTT elements =matrix, [mji  will the row transpose ONR N00014-16-1-2645. ThisYIP research is supported supported by by NSF NSF CAREER CAREER ECCS-1454022 ECCS-1454022 and and 1 , . denote n ) and ji ] This research is supported by NSF CAREER ECCS-1454022 and 1 n ji will denote the row vector and the transpose matrix, ONR YIP N00014-16-1-2645.  m , . . . , m . Also, y = (y , . . . , y ) and M = [m ] will denote the row vector and the transpose matrix, ONR YIP N00014-16-1-2645. This research is supported by NSF CAREER ECCS-1454022 and 1 n 1 n ji will denote the row vector and the transpose matrix, ONR YIP N00014-16-1-2645. ONR YIP N00014-16-1-2645. will denote the row vector and the transpose matrix, ONR YIP N00014-16-1-2645. 2405-8963 © 2018, IFAC (International Federation of Automatic Control) Hosting by Elsevier Ltd. All rights reserved.

C. C. C. C.

Copyright © 2018 IFAC 142 Copyright 2018 142 Peer review© under responsibility of International Federation of Automatic Control. Copyright © 2018 IFAC IFAC 142 Copyright © 142 Copyright © 2018 2018 IFAC IFAC 142 10.1016/j.ifacol.2018.12.025 Copyright © 2018 IFAC 142

IFAC NecSys 2018 Groningen, NL, August 27-28, 2018

C. Somarakis et al. / IFAC PapersOnLine 51-23 (2018) 142–147

respectively. The n × n zero and identity matrices are denote by On and In , respectively. The n × 1 vector 1n is the vector of all ones.   An undirected, weighted graph is the triple G = V, E, w , with V the set of n vertices, E the set of edges and w = {wij } the set of edge-weights with: wij = 0 whenever i, j ∈ V are not connected by an edge, and wij = wji > 0 otherwise. We assume that wii ≡ 0. The graph Laplacian LG = LTG is the square  matrix with elements lij = −wij if i = j and lii = j wij . We will be considering only connected graph topologies. In such case we can write LG = QΛQT with Λ = {λ1 , . . . , λn }, 0 = λ1 < λ2 ≤ · · · ≤ λn and Q = [q1 | . . . | qn ] is an orthonormal matrix of the eigenvectors of LG . In particular, q1 = √1n 1n . The complete incidence matrix of n(n − 1)/2 × n dimensions is denoted by Bn . Finally, by y ∼ N (µ, σ 2 ) we understand a random variable y that is normally distributed with mean value µ and variance σ 2 . 3. POWER NETWORK DYNAMICS Let a network of n ≥ 1 synchronous generators (vertices) over m transmission lines (edges). The i’th generator attains the positive triple of parameters (Ji , βi , Ei ) consisting of the rotational inertia, the damping constant, and the voltage magnitude, respectively. The dynamic phase of voltage Ei is denoted by θi = θi (t). The phase (rotor angle) d θi =: ωi characterize together with the rotor frequency dt the dynamic state of generator i. The n generators drive (through a transient reactance structure) the electric network that consists of transmission lines and loads modeled as constant impedances. The power network topology is represented via the reduced complex admittance matrix Y = [Yij ∠φij ] where the angles φij ∈ [0, π2 ] are the phase shift representing the energy loss due to transfer conductance. In this work we assume that the transfer conductances are all zero, i.e. φij = π2 for i = j. Due to space limitations, proofs of technical results are either omitted or briefly sketched.

3.1 The Classical Model The dynamic state of the i’th generator is dictated by the nonlinear differential equation Ji θ¨i = −βi θ˙i +

n  j=1

Ei Ej Yij sin(θj − θi ) + pi

(1)

where pi = Pim − Ei2 {Yii } is the effective power input, i.e. defined as the difference between the mechanical input power and the electrical internal voltage power. Model (1) has been extensively used for wide-area control of power networks (see for example Eq. (2.3) of Varaiya et al. (1985) or Poolla et al. (2017); Vu and Turitsyn (2017) and references therein) The state of the network is then characterized through the vector (θ, ω)T ∈ R2n . The main objective of transient stability is to investigate the conditions under which (1) will converge to a desirable   operating equilibrium θ (s) , ω (s) ∈ S with

143

S=





 θ, ω ∈ R2n : ω = 0 pi =

n  j=1

143

and

 Ei Ej Yij sin(θi − θj ), i = 1, . . . , n

As Varaiya et al. (1985) point out, model (1) makes unwanted simplifications in its formulation of the transmission lines and loads on the one hand, and of the generators on the other. Despite the drawbacks, (1) is a rich enough benchmark, for the study of post-fault power stability. 3.2 Linear Approximations Let (θ ∗ , 0)T ∈ S be a desirable equilibrium. The dynamics of (1) operating around (θ ∗ , 0)T can be determined by its first order approximation. It can be readily determined that the steady-state power flow equations relating the effective power input p = (p1 , . . . , pn )T relates to θ ∗ via p = LG θ ∗ where LG = [lij ] the graph laplacian with off-diagonal elements lij = J1i Ei Ej Yij . We note that in view of the symmetry of Y , LG is symmetric by construction. The linearized (1) reads n  βi lij (θj − θi ) + pi (2) θ¨i = − θ˙i + Ji j=1 The next condition is a widely used assumption in the literature (see Tegling et al. (2015); Siami and Motee (2016); Poolla et al. (2017)), that significantly simplifies the study of (2). Assumption 1. All synchronous generators have uniform damping and inertia parameters, βi ≡ β > 0, Ji ≡ J > 0. In view of the latter condition we define the damping over inertia ratio βi /Ji ≡ d > 0.

Furthemore, we assume that the network effective power injection p is prone to disturbances. These are envisioned as n independent stochastic sources of white noise dξ t = (1) (n) (dξt , . . . , dξt )T infused in the power network as a disturbance of additive type. The resulting dynamics in vector formbecome      On dθ t On In θt dξ t = dt + −1 (3) dω t −LG −D ω t J In

where D = diag({d}). While (3) is a higly stylized approximation of (1), it is considered as a reliable approximation and a fruitful source for meaningful analysis (see for example Eliasson and Hill (1992)). 3.3 Distributed Delayed Feedback Control

Disturbances typically steer the system dynamics away from the nominal operating conditions and into uncertainty. Rotor angles can diverge to the point of challenging the validity of the linear surrogate (3) or even threatening the stability of the overall power network. We propose to stabilize the network via a state-feedback law, assumed to act on the communication layer, i.e. through higher level controllers that attempt to coordinate the generators without any intervention on the

IFAC NecSys 2018 144 Groningen, NL, August 27-28, 2018

C. Somarakis et al. / IFAC PapersOnLine 51-23 (2018) 142–147

power network admittance. The control process existentially bypasses the electric grid, but it comes at the cost of communication time-delay that is assumed here to be constant and uniform. More specifically, (3) with the additive control term, ut , reads dθ t = ω t dt   dω t = − LG θ t − Dω t + ut dt + dξ t . The distributed state-feedback law takes the form ut = −Lθ θ t−τ − Lω ω t−τ (4) for Lθ = κθ LG , and Lω = κω LG , where κθ , κω > 0 are control parameters, and θ t−τ and ω t−τ are lagged state observations, with some time-delay τ > 0. The feedback law (4) applied to (5) yields the system under investigation in this paper. This is posed below as an initial value problem      dθ t On In θt dt− = −LG −D ω t dω t      (5) O On On θ t−τ − dt + −1n dξ t , κθ LG κω LG ω t−τ J In   for t ≥ 0, under initialization φ(t) = φθ (t), φω (t) t ∈ [−τ, 0], that is assumed statistically independent from the white noise. 1 4. PRELIMINARIES The assessment of risk of systemic events in dynamical systems requires two preliminary stages. The first stage regards investigation of the stability region of the unperturbed closed loop network, i.e. (5). The second stage requires exploration of the statistical behavior of an output observable vector y ∈ Rn(n−1) . In this work, we take as output observables the differences of the rotor angles y = Bn θ. Any element of y is of the form θ(i) − θ(j) . Our objective will be to introduce systemic sets for the elements of y, as red-zone ranges of values. Then we will assess the risk of the rotor angles to approach or land on such undesirable ranges. We will highlight the interplay between spectral properties of the network, the feedback parameters and the delay. 4.1 Internal Stability The uncontrolled deterministic network (3) is exponentially stable wrt the origin 2 , when no disturbance is assumed. Implementing a feedback law with delays may jeopardize the stability the network. The next result is considered to be the first result of this work. For its exposition we define ∆d,κω (λ) := (κω λ)2 + 2λ − d2 , and    1 2 2 2 2 ∆d,κω (λ) ± ∆d,κω (λ) − 4λ (1 − κθ ) , (6) γ± (λ) := 2 whenever the right-hand side delivers a non-negative number. Next, ϑ± = ϑ± (λ) ∈ [0, 2π) determined through 2 2 κω λdγ± + κθ λ(λ − γ± ) cos(ϑ± ) = − 2 2 2 2 λ (κω γ± + κθ ) 2 ) κθ λdγ± − κω λγ± (λ − γ± . sin(ϑ± ) = 2 2 2 2 λ (κω γ± + κθ ) 1

For the well-posedness of such initial value problems we refer to Mohammed (1984). 2 That is wrt (θ ∗ , 0)T after appropriate axes translation.

144

For parameter values in (6) that deliver two distinct solutions it can be verified that γ+ > γ− . Then for the values 2π ϑ− 2π ϑ+ τk,+ := +k and τk,− := +k , k = 0, 1, 2, . . . γ+ γ+ γ− γ− it can be easily deduced that τk+1,+ − τk,+ < τk+1,− − τk,− and that there exists at least one index k ∗ > 1 such that τk∗ −1,− < τk∗ ,+ < τk∗ +1,+ < τk∗ ,− . We are now ready to state the first result of this work. Theorem 2. Assume the dynamics of the unperturbed verT  sion of (5). The solution θ t , ω t , t ≥ 0 is exponentially stable if and only if for all non-zero eigenvalues of LG , λi , i > 1, one of the following conditions is true: (1) κθ ≥ 1 and τ ∈ [0, τ0,+ ). (2) 0 ≤ κθ < 1, ∆d,κω (λi ) > 0, ∆2d,κω (λi ) > 4λ2i (1 −   κ2θ ), so that τ0,+ < τ0,− and τ ∈ 0, τ0,+ or τ ∈  k∗ −1  τk,− , τk+1,+ . k=1 (3) 0 ≤ κθ < 1 and ∆d,κω (λi ) < 0 or ∆2d,κω (λi ) < 4λ2i (1 − κ2θ ). Proof. [Sketch] The symmetry condition on the grid, implies that LG is symmetric as well, can be written as LG = QΛQT . Invoking the transformation z := QT θ, we obtain the decoupled dynamics that can be written as the following system of scalar second order delayed differential equations (i)

(i)

(i)

(i)

(i)

z¨t + λi zt + d z˙t + (κθ λi ) zt−τ + (κω λi ) z˙t−τ = 0 (7) for i = 1, . . . , n. For i = 1, we have z¨(1) = −dz˙ (1) , (1) equivalent to z˙t = e−dt z˙0 (1) . For i > 1 the delayed terms contribute in the dynamics, as well. The proof investigates the existence of solutions of the form eht , h ∈ C for (7). The result establishes necessary and sufficient conditions under which {h} < 0. The derivation of all three stability conditions involves lengthy algebraic manipulations that are omitted due to space limitations. We remark that Condition (3) of Theorem 2 is a delayed independent stability result. For instance, if damping over inertia ratio d > 0 is large enough, the implemented feedback law poses no threat to the stability of the solution. It may however affect the performance of the network in the face of disturbances. In any case, if Φ(t) is the fundamental solution of (5), then Theorem 2 provides necessary and sufficient conditions for limt→+∞ Φ(t) = O2n .

We remark that the results of Theorem 2 hold for any type of distributedfeedback law that attains the form  (i)  T Lθ = Qdiag {κθ } QT and Lω = Qdiag {κ(i) ω } Q , where Q is the eigenvector matrix of LG . For the sake of (i) simplicity, we restrict our study to the case κθ := κθ λ(i) (i) and κω := κω λ(i) , i ≥ 1. 4.2 State and Output Vector Statistics The main implication of Theorem 2 is that the full scale dynamics (5) generate a stochastic process (θ t , ω t )T , t ≥ 0, that can be represented with the use of Itˆo calculus      t On θt Φ(t − s) −1 (8) = Tt + dξ s . ωt J In 0

IFAC NecSys 2018 Groningen, NL, August 27-28, 2018

C. Somarakis et al. / IFAC PapersOnLine 51-23 (2018) 142–147

where qlk is the lth element of the eigenvector qk that corresponds to the k th eigenvalue of LG . The latter observation characterizes the variability of the phase difference. The distributional properties constitute the basis for the risk evaluation of systemic events that will be introduced and discussed in the next section. 5. SYSTEMIC SETS & EVALUATION OF RISK

Table 1. The vector α(λ) ˆ of Theorem 3. i 0 1 2 3 4

(P )

α ˆ i (λ)/τ 4−i λ2 (κ2θ + 1) 0 d2 − 2λ + (κω λ)2 0 1

(C)

(λ)/τ 4−i 2κθ λ2j 0  2(d κω − κθ λ 0 0 α ˆi

145

(S)

(λ)/τ 4−i 0 2λ(κω λ − d κθ ) 0 −2κω λ 0 α ˆi

0 The transient term Tt = −τ Φ(t − s − τ )φ(s) dµ(s) is represented in a compact form as a Stieltjes integral for some appropriate measure µ of bounded variation (see Mohammed (1984)). For any t < ∞, (θ t , ω t )T is normally distributed   as      t O O θ Φ(s) n n ΦT (s) ds . ∼ N Tt , On I n ω 0 The calculation of the covariance matrix is not feasible in the time-domain. The next results exploits the frequency domain methodology to provide some insight when t → ∞ for the vector of phase incoherence yt = Bn θ t where Bn is the complete incidence matrix, that it is of interest in this work. For the exposition of the result we define the multivariate function f (α) = f (α(P ) , α(C) , α(S) ) : W ⊂ R15 → R+ with  dr f (α) = . 4  (p) i (c) i (s) i R i=0 αi r + αi r cos(r) + αi r sin(r) The domain of definition W includes all points in R15 that make the integral in f to converge. Conditions of Theorem 2 assure that W is non-empty and open. Theorem 3. Assume either condition of Theorem 2 to hold. For the solution (θ t , ω t )T of (5), the output vector yt = Bn θ t converges as t → ∞ in distribution to     T T τ3 y ∼ N 0, B Q diag f ( α(λ ˆ )) Q B n k n 2πJ 2 where α(λ) ˆ is the vector, the elements of which are listed in Table 1, and λk are the eigenvalues of the graph Laplacian matrix LG . Proof. [Sketch] Recall the transformation z(1) = QT θ and z(2) = QT ω, that fully diagonalizes (5). In view of the Assumptions considered so far, z1 can be shown to  ∞ be a Gaussian with covariance matrix diag 0 h2i (t) dt , for hi (t) representing the transient response from the ith element of z (2) to the ith element of z(1) . Then Parseval’s theorem is utilized to calculate the ith diagonal element through its frequency domain form. Indeed  ∞  1 d 2 hi (t) dt = 2 2π χ (j ) i R 0 2 −sτ where χi (s) = s + d s + λi κω se + λi (1 + κθ e−sτ ), is the characteristic function of the i’th decomposed state. The result follows in view of y = Bn Qz(1) and application of basic properties of normal distributions (see Tong (1990)). A direct implication of Theorem 4 is that every element of y, being in the form θ(i) − θ(j) for some i, j ∈ V, follows a normal distribution with parameters   n   τ3  (i) (j) 2 θ −θ ∼ N 0, (qik − qjk ) f α(λ ˆ k) . 2πJ 2 k=2

145

Network (3) is an efficient approximation of (1) in the vicinity of (θ ∗ , 0) ∈ S. On the other hand, when variation of θ(i) − θ(j) exceeds certain margins, nonlinear terms (e.g. higher order terms that are present in (1)) may prevail on the linear dynamics steering the system away from the desirable equilibrium. 5.1 Systemic Sets We proceed with the setting of desirable, unsafe and undesirable ranges of values for the elements of y, y i   that corresponds to a difference θ(i ) − θ(i ) between two synchronous generators. Clearly, values in the vicinity of the origin are desirable. Let η > 0 be a cut-off beyond which |y (i) | should be not be found. We then pick c ≥ 1 we define the intermediate cut-off η/c. Any value of |y (i) | in (0, η/c) will be considered safe and acceptable, any value of |y (i) | in (η, ∞) will be marked as unsafe, and finally, any value of |y (i) | in (η/c, η) should trigger a warning that the system state may exceed the safety cut-off η. An example of such safety profile is motivated by the form of (1) when compared to the linear approximation (5). It is clear that the dominating linear term is highly challenged for any directly pair i, j with difference values |θ(i) − θ(j) | to exceed π/6, whereas (5) becomes obsolete for |θ(i) −θ(j) | beyond π/2. For fixed safety parameters η > 0, c ≥ 1 we consider the δ > 0-parametrized collection of sets defined as   1+δ Uδ = η , +∞ . c+δ The two extremes of collection {Uδ }δ , i.e. U0 = (η/c, +∞) and U∞ = (η, ∞) signify the warning set and the (redzone) systemic set of our configuration, respectively. The parameter values η and c are also known in the literature as security or stability constraints. In the next section, we will leverage quantile-based measures of risk, and apply Theorem 3 to evaluate the possibility of the aforementioned systemic events. 5.2 Risk Evaluation Our objective is to provide closed form expressions of the risk of an element of y to develop a systemic event (i.e. to reach close to or on the systemic set U∞ ). The value-at-risk metric we will use in this work is   (9) Rε (|y i |) = inf δ > 0 | P (|y (i) | ∈ Uδ ) < ε for some confidence level ε ∈ (0, 1). The choice of ε determines the rigidity of the risk index. The closer ε in 1, the less reliable, the index becomes. On the other hand for ε close to zero, it may become unnecessarily strict. The next result provides a closed form expression of the risk of systemic events generated by extreme deviatiosn of

C. Somarakis et al. / IFAC PapersOnLine 51-23 (2018) 142–147

k=2

The essence of Theorem 4 is readily captured should one η 1 view Rij ε as a function of σij . When σij lies below c κε , the risk of critical deviation is zero. In this range of values, the event of observing the absolute difference of generators i and j phases, above the first cut-off, is smaller than ε. For small level of confidence ε, zero value of risk implies that with very high probability the absolute difference lies within the range of acceptable fluctuations. A positive (but finite) risk value implies that the absolute difference will fluctuate dangerously close to the systemic set with an unacceptable high probability. For large of value of σij we obtain an infinite risk value that signifies the occurrence of systemic events with dangerously high probability. Theorem 4 can be successively applied for any pair of generators of interest or for all of them, {Rij ε }i
λ= λ= λ= λ= λ= λ= λ= λ= λ= λ=

1.2

1

0.8

σ

the rotor phase between two fixed generators. Its proof is totally omitted due to space limitations. Theorem 4. The risk of systemic events for the dynamics of (5) with output observable |θ(i) − θ(j) | is in the equilibrium statistics evaluated as  η 1  0, if σij ≤   c κε    1 − σij ηc κε η 1 1 ij if < σij < η (10) Rε = 1 c κε κε   σij η κε − 1   1  +∞ if σij ≥ η κε √ κ 2 where κε is the solution of −κε ε e−t /2 dt = 2π(1 − ε), and marginal variance n   τ3  2 = (qik − qjk )2 f α(λ ˆ k) . σij 2 2πJ

0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0

0.6

0.4

0.2 0

2

4

6

8

10

12

τ

7 λ= λ= λ= λ= λ=

6

5

R0.1

IFAC NecSys 2018 146 Groningen, NL, August 27-28, 2018

1.0 1.5 2.0 2.5 3.0

4

3

2

1

0 2

4

6

8

10

12

14

τ

Fig. 1. Evolution of marginal variance σ and the associated risk with acceptance level ε = 0.1, as a function of τ . Plots are for various λ. We remark that R0.1 of λ = 0.5 is infinite, and R0.1 of any λ > 3.0 is zero.

enough, and it disappears as they attain higher values. We begin our exploration with phase feedback control and we will proceed with frequency feedback, to conclude with a case of homogeneous control (i.e. κω = κθ = κ).

6.1 Stated Feedback, Network, and Time-Delay Interplays. We separate our study into schemes of phase feedback control (i.e. κω = 0) and schemes of frequency feedback control (i.e. κθ = 0). Delay-Independence. In the absence of frequency control, Theorem 2 distinguishes between κθ ∈ (0, 1) and κθ > 1. We investigate this case in Figure 1, where for κω = 0 we set κθ = 0.5 and tested the behavior of σ and R0.1 as functions of time-delay τ and different λ values that vary from 0.5 to 5.0. Marginal variance increases up to an asymptotic value, and risk behaves similarly. Delay-Dependence. Critical delays can appear in either of the two schemes we investigate. In the phase control, maximal acceptable delays emerge when κθ > 1, while in the frequency control such boundaries exist for arbitrary range in the gain κω > 0. We illustrate relevant results in Figure 2. In the first plot we have phase control with κθ = 1.5. Time-delay is assumed constant to τ = 0.1 and we varied λ from 1.2 to 19.5. The curves of σ and risk are plotted together for convenience. While both types of control provide qualitatively similar behavior, we remark that the phase control delivers zero risk for a wider range of network parameter λ.

3

In this paper, we focus on scalar events. Extensions to multiple events can be done by calculating bounds for the risk measure, which is beyond the scope of this paper. For a detailed discussion of vectors of risk metrics and distributional properties of multiple joint events we refer to Somarakis et al. (2018b).

146

6.2 Combined Feedback and Delay-Induced Limitations Now we will numerically explore the effect of feedback law, of the form κθ = κω = κ > 0 that is prone to time-

IFAC NecSys 2018 Groningen, NL, August 27-28, 2018

C. Somarakis et al. / IFAC PapersOnLine 51-23 (2018) 142–147

147

deeper investigation. Our results hold promise for the development of optimal control design techniques on the basis of risk metrics. REFERENCES

1.5

σ R0.005

σ, R0.005

1

0.5

0 0

2

4

6

8

10

12

14

16

18

20

λ

σ R0.005

0.2

σ, R0.005

0.15

0.1

0.05

0 3

4

5

6

7

8

9

10

11

12

13

λ

Fig. 2. Delay-Dependent Control. Up : κθ = 1.5, κω = 0. Down: κθ = 0 and κω = 1.5. σ R0.0001

0.5

σ, R0.0001

0.4

0.3

0.2

0.1

0 0

1

2

3

4

5

κ

Fig. 3. Combined control with identical gains κω = κθ = κ. Here λ = 3, τ = 0.1. delay value. We fix τ = 0.1, λ = 3 and plot (11) and the associated risk. The results are presented in Figure 3. The feedback law ensures a zero risk at the 0.0001 level of confidence within the range [1.88, 5.65]. Stronger control gains yield higher values of risk. 7. CONCLUSION Rotor phase deviation is a critical quantity to assess power grid robustness against disturbances. We developed a riskoriented framework to study the interplay between grid parameters, stochastic perturbations and control strategies. We showed that risk depends on the variability of the phase deviation that, in turn, encapsulates network connectivity, feedback gains and time-delay. Theoretical findings, as well as numerical explorations reveal a perplexed connection among these parameters. In this first round of results we indicate possible fundamental tradeoffs among network connectivity (that depends on feedback gains), time-delay and the risk phase incoherence. This implies that using larger feedback gains may result in risk of phase deviation. This is a point that requires 147

Bamieh, B., Jovanovi´c, M., Mitra, P., and Patterson, S. (2012). Coherence in large-scale networks: Dimensiondependent limitations of local feedback. Automatic Control, IEEE Transactions on, 57(9), 2235 –2249. D¨orfler, F. and Bullo, F. (2012). Synchronization and transient stability in power networks and nonuniform kuramoto oscillators. SIAM Journal on Control and Optimization, 50(3), 1616–1642. Eliasson, B.E. and Hill, D.J. (1992). Damping structure and sensitivity in the nordel power system. IEEE Transactions on Power Systems, 7(1), 97–105. Fardad, M., Lin, F., and Jovanovic, M.R. (2014). Design of optimal sparse interconnection graphs for synchronization of oscillator networks. IEEE Transactions on Automatic Control, 59(9), 2457–2462. Mohammed, S.E.A. (1984). Stochastic Functional Differential Equations. Pilman Advanced Publishing. Paganini, F. and Mallada, E. (2017). Global performance metrics for synchronization of heterogeneously rated power systems: The role of machine models and inertia. In 2017 55th Annual Allerton Conference, 324–331. Poolla, B.K., Bolognani, S., and D¨orfler, F. (2017). Optimal placement of virtual inertia in power grids. IEEE Transactions on Automatic Control, 62(12), 6209–6220. Sauer, P. and Pai, M. (2006). Power System Dynamics and Stability. Stipes Publishing L.L.C. Siami, M. and Motee, N. (2016). Fundamental limits and tradeoffs on disturbance propagation in linear dynamical networks. IEEE Transactions on Automatic Control, 61(12), 4055–5062. Somarakis, C., Ghaedsharaf, Y., and Motee, N. (2017). Aggregate fluctuations in time-delay linear consensus networks: A systemic risk perspective. In The 2017 American Control Conference. Somarakis, C., Ghaedsharaf, Y., and Motee, N. (2018a). Risk of collision in a vehicle platoon in presence of communication time delay and exogenous stochastic disturbance. In The 57th IEEE CDC. (to appear). Somarakis, C., Ghaedsharaf, Y., and Motee, N. (2018b). Time-delay origins of fundamental tradeoffs between risk of large fluctuations and network connectivity. CoRR, abs/1801.06856. URL http://arxiv.org/abs/1801.06856. Summers, T., Shames, I., Lygeros, J., and Drfler, F. (2015). Topology design for optimal network coherence. In 2015 European Control Conference (ECC), 575–580. Tegling, E., Bamieh, B., and Gayme, D.F. (2015). The price of synchrony: Evaluating the resistive losses in synchronizing power networks. IEEE Transactions on Control of Network Systems, 2(3), 254–266. Tong, Y.L. (1990). The Multivariate Normal Distribution. Springer, 1st edition. Varaiya, P., Wu, F.F., and Chen, R.L. (1985). Direct methods for transient stability analysis of power systems: Recent results. Proceedings of the IEEE, 73(12), 1703–1715. Vu, L.T. and Turitsyn, K. (2017). A framework for robust assessment of power grid stability and resiliency. IEEE Transactions on Automatic Control, 62(3), 1165–1177.