Close
About
FAQ
Home
Collections
Login
USC Login
Register
0
Selected
Invert selection
Deselect all
Deselect all
Click here to refresh results
Click here to refresh results
USC
/
Digital Library
/
University of Southern California Dissertations and Theses
/
Data-driven aggregated fractal dynamic load modeling for early warning of voltage collapse
(USC Thesis Other)
Data-driven aggregated fractal dynamic load modeling for early warning of voltage collapse
PDF
Download
Share
Open document
Flip pages
Contact Us
Contact Us
Copy asset link
Request this asset
Transcript (if available)
Content
DATA-DRIVEN AGGREGATED FRACTAL DYNAMIC LOAD MODELING FOR EARLY WARNING OF VOLTAGE COLLAPSE by Laith Shalalfeh A Dissertation Presented to the FACULTY OF THE USC GRADUATE SCHOOL UNIVERSITY OF SOUTHERN CALIFORNIA In Partial Fulfillment of the Requirements for the Degree DOCTOR OF PHILOSOPHY (ELECTRICAL ENGINEERING) May 2017 Copyright 2017 Laith Shalalfeh Dedication To my parents, my wife, and my children for their love and support. ii Acknowledgments First and foremost, all praises and gratitude to Allah for all strengths and blessings that he bestowed on me to complete my research work successfully. I would like to thank my thesis advisor, professor Edmond Jonckheere, for giv- ing me the opportunity to be a doctoral student under his supervision. During my graduate study, professor Jonckheere was a great teacher, mentor, and friend. He was always available to help and ready to answer any difficult questions faced in my research. Besides that, our regular group meetings were very important to succeed in research because of the excellent comments and feedbacks provided by my advisor and our group members. Truly, my graduate studies would be much harder without his guidance. I would like also to thank my co-advisor, professor Paul Bogdan, for his super- vision during my work on statistical analysis of phasor measurement unit (PMU) data. I have learned a lot from him during our regular meetings and discussions. I am grateful to professor Mohammed Beshir for the very helpful comments on the work related to voltage collapse and load modeling. His insight and practical experience consistently added a great value to my research work. iii I am also grateful to my parents, Abdel-majeed and Jehad, for encouraging me to continue my graduate study. They have always inspired me to be a better person. Thank you for your unconditional love and support. I would like to thank my wife, Safa’a, who stood beside me during the graduate study journey. I have to admit that having her made my graduate study journey much easier. My children, Miral and Yahia, were constantly giving me hope and courage when I needed it the most. iv Contents Dedication ii Acknowledgments iii List of Tables vii List of Figures viii Abstract xii Introduction xiv 1 Load Modeling 1 1.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.2 Berg [7] Model versus Hill Model [33] . . . . . . . . . . . . . . . . . 2 1.3 Static-Dynamic Gap Load Modeling . . . . . . . . . . . . . . . . . . 4 1.4 Power System Model . . . . . . . . . . . . . . . . . . . . . . . . . . 7 1.5 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8 2 Voltage Collapse in the Power System 9 2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9 2.2 Existence of a Voltage Collapse Solution . . . . . . . . . . . . . . . 10 2.2.1 Analyticity . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 2.2.2 Voltage Collapse Condition . . . . . . . . . . . . . . . . . . 12 2.3 The possibility of voltage collapse at different loads . . . . . . . . . 14 2.3.1 Synthetic loads with p v =q v ,p w =q w . . . . . . . . . . . . . 15 2.3.2 Berg load models . . . . . . . . . . . . . . . . . . . . . . . . 17 2.3.3 Correct initial conditions . . . . . . . . . . . . . . . . . . . . 21 2.4 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22 3 Load Aggregation Effect in Power Grid 24 3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24 3.2 Berg load modeling and fractional dynamics . . . . . . . . . . . . . 25 v 3.3 Passivity of Berg loads . . . . . . . . . . . . . . . . . . . . . . . . . 26 3.3.1 Synthetic loads . . . . . . . . . . . . . . . . . . . . . . . . . 26 3.3.2 Berg loads . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27 3.4 Basic motif . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29 3.5 Basic motif interconnection . . . . . . . . . . . . . . . . . . . . . . 29 3.6 Aggregation effect: one single Strongly connected component . . . . 32 3.7 The effect of contingencies on the power system graph . . . . . . . . 38 3.7.1 Single-element contingency . . . . . . . . . . . . . . . . . . . 38 3.7.2 Double-element contingency . . . . . . . . . . . . . . . . . . 39 3.8 Formal result . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41 3.9 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41 4 Evidence of Long-Range Memory in Power System 43 4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43 4.2 Prior work and novel contribution . . . . . . . . . . . . . . . . . . . 44 4.3 Statistical properties of the power system . . . . . . . . . . . . . . . 45 4.3.1 Description of the power system measurements . . . . . . . . 45 4.3.2 Statistical analysis of PMU measurements . . . . . . . . . . 47 4.4 Existence of long-range dependence . . . . . . . . . . . . . . . . . . 49 4.4.1 Detrended Fluctuation Analysis (DFA) . . . . . . . . . . . . 49 4.4.2 Analysis of the power system data using DFA . . . . . . . . 50 4.4.3 DFA analysis of the surrogate power system data . . . . . . 54 4.5 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55 5 Kendall’s Tau of Frequency Hurst Exponent as Blackout Proxim- ity Margin 56 5.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56 5.2 Related work and novel contribution . . . . . . . . . . . . . . . . . 57 5.3 Power system data . . . . . . . . . . . . . . . . . . . . . . . . . . . 58 5.3.1 Description of the PMU data . . . . . . . . . . . . . . . . . 58 5.3.2 PMU data fractality . . . . . . . . . . . . . . . . . . . . . . 60 5.4 Kendall’s tau as Blackout Proximity Margin . . . . . . . . . . . . . 61 5.4.1 Kendall’s tau of AR(1) . . . . . . . . . . . . . . . . . . . . . 62 5.4.2 Kendall’s tau of the Hurst exponent as dynamic proximity margin to blackout . . . . . . . . . . . . . . . . . . . . . . . 64 5.4.3 Kendall’s tau of AR(1) coefficient versus Kendall’s tau of the Hurst exponent under normal conditions . . . . . . . . . . . 65 5.5 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67 6 Conclusion 69 Reference List 70 vi List of Tables 1.1 Typical values of characteristic-load parameters . . . . . . . . . . . 5 1.2 Impedances of various loads . . . . . . . . . . . . . . . . . . . . . . 6 4.1 Scaling exponents of voltage magnitude, frequency, and phase angle for four data sets at three locations: Baylor University, Harris Sub- station, and McDonald Observatory . . . . . . . . . . . . . . . . . . 53 vii List of Figures 1.1 Load Models: (a) Hill Model, (b) Berg Model . . . . . . . . . . . . 3 1.2 Power system model: (a) Circuit diagram, (b) Feedback model . . . 7 1.3 Block Diagram of the harmonic oscillator . . . . . . . . . . . . . . . 8 2.1 Block diagram of the power system including the harmonic oscillator 10 2.2 Voltage collapse region . . . . . . . . . . . . . . . . . . . . . . . . . 16 2.3 Voltage collapse phenomnon : (a) Voltage magnitude (V), (b) Phase angle of the generator bus (δ), (c) Frequency (ω) . . . . . . . . . . . 17 2.4 Sigma (σ) and frequency (w) for different synthetic loads . . . . . . 17 2.5 Sigma (σ) and frequency (w) for induction motor loads (half/full) . 18 2.6 Sigma (σ) and frequency (w) for regulated aluminum plant . . . . . 19 2.7 The relationship between the transmission line coefficient (K Line ) and Sigma (σ) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20 2.8 Therelationshipbetweentheactivepowercoefficient(K p )andSigma (σ) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21 2.9 The relationship between the reactive power coefficient (K q ) and Sigma (σ) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21 viii 3.1 Realpartoftheadmittance(Y)asfunctionofsigma(σ)fordifferent Berg loads: (a) Filament lamp (b) Fluorescent lamp (c) Induction motor-halfload(d)Inductionmotor-fullload(e)Reductionfurnace (f) Aluminum Plant . . . . . . . . . . . . . . . . . . . . . . . . . . . 28 3.2 Realpartoftheadmittance(Y)asfunctionofvoltage(V)andsigma (σ) for different Berg loads: (a) Filament lamp (b) Fluorescent lamp (c) Induction motor- half load (d) Induction motor- full load (e) Reduction furnace (f) Aluminum Plant . . . . . . . . . . . . . . . . 28 3.3 Simple 2-bus power system: (a) Circuit diagram (b) Feedback model 29 3.4 Hypothetical decomposition of a grid digraph into strongly con- nected components. Each strongly connected component (dotted circle)isitselfamultivariableversionofthebottommotifofFig.3.3(b). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30 3.5 3-bus power system: 1 generator, 2 loads . . . . . . . . . . . . . . . 33 3.6 3-motif case (1 generator, 2 loads): (a) Feedback model, (b) Graph model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33 3.7 6-bus power system . . . . . . . . . . . . . . . . . . . . . . . . . . . 34 3.8 6-motif case: (a) Feedback model, (b) Graph model . . . . . . . . . 36 3.9 IEEE 14 bus power system . . . . . . . . . . . . . . . . . . . . . . . 37 3.10 Graph model in Small World configuration . . . . . . . . . . . . . . 37 3.11 Graph model of 6-bus power system after disconnecting line 5-6 . . 39 3.12 Graph model of 6-bus power system after a 3-phase fault at bus 2 . 39 3.13 Graph model of 6-bus power system after disconnecting lines 2-3 and 5-6 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40 3.14 Graph model of 6-bus power system after faults at buses 2 and 5 . 40 ix 4.2 (a) The voltage magnitude at Baylor University on 05/27/15 (12:00- 1:00PM)(b)ThefrequencyatBaylorUniversityon05/27/15(12:00- 1:00PM)(c)ThevoltagephaseangleatBaylorUniversityon05/27/15 (12:00-1:00 PM) (d) The unwrapped phase angle (e) The empirical CDFs of the voltage magnitude at four intervals (900 seconds each) (f) The empirical CDFs of the frequency at four intervals (900 sec- onds each) (g) The empirical CDFs of the phase angle at four inter- vals (900 seconds each) (h) The empirical CDFs of the unwrapped phase angle at four intervals (900 seconds each) . . . . . . . . . . . 47 4.3 Log-log plots of the rms fluctuation function F (n) versus the box sizen ofV,f, andθ of different data sets: (a) Data set 1 at Baylor University (b) Data set 2 at Baylor University (c) Data set 3 at Baylor University (d) Data set 4 at Baylor University (e) Data set 1 at Harris Substation (f) Data set 2 at Harris Substation (g) Data set 3 at Harris Substation (h) Data set 4 at Harris Substation . . . 52 4.4 Log-log plots of the rms fluctuation function F (n) versus the box size n of original and surrogate data at: (a) Baylor University (b) Harris Substation (c) McDonald Observatory . . . . . . . . . . . . . 54 5.1 Analysis of PMU#4 data recorded at EPFL on 01/15/2015 (8:00- 8:33 AM) during normal operation. First row: raw data recorded; second row: corresponding log-log plot of rms fluctuation function versustheboxsizen. First(red)column: voltagemagnitude; second (blue) column: frequency; third (green) column: phase angle . . . . 59 5.2 Hurst exponent histograms of EPFL PMU#4 data during all 4 sea- sons (January, April, June, and October): (a) Voltage magnitude (b) Frequency (c) Phase angle . . . . . . . . . . . . . . . . . . . . . 61 x 5.3 (a) Raw (blue) and Gaussian kernel smoothed (red) frequency time series before 2012 Indian blackout (b) Autocorrelation coefficient at lag 1 of residual differences between raw and smoothed frequency time series (c) Hurst exponent of the raw frequency time series . . . 62 5.4 (a) Box-plots of Kendall’s tau of frequency AR(1) coefficient for data collected in EPFL during January, April, June, and October (b)HistogramofKendall’stauoffrequencyAR(1)coefficientduring these four months (c) Box-plots of the Kendall’s tau of frequency Hurst exponent for data collected in EPFL during January, April, June, and October (d) Histogram of Kendall’s tau of frequency Hurst exponent during these four months. Due to blackout data scarcity, a red histogram and a test for random draws from different distributions couldn’t be obtained. . . . . . . . . . . . . . . . . . . 66 5.5 Decision tree acting on the multi-layer decomposition of Fig. 5 of [27] as response to increase of Hurst exponent . . . . . . . . . . . 68 xi Abstract Voltage collapse, one of the critical phenomena that threatens the power infras- tructure, manifests itself by a sudden and fast collapse of the system voltage. The simplest bus system—one generator, one line, and one load—reveals a feedback structure, which, depending on the nature of the load, could go in a voltage col- lapse resembling an “anti van der Pol" behavior. A static-dynamic compromise for the load model is provided by an “impedance" describing function that has the advantage of being derived from experimental data collected in a real grid environ- ment. This reveals, among other things, a noninteger exponent of the frequency requiring an extension for the model to be usable for a collapsing solution away from the imaginary axis. We argue that an aggregation effect of all of the loads in a grid endows a single load with grid-characteristics properties in addition to the usual load-specific properties. Topologically, thehiddenfeedbackstructureofabusmodelrevealsthat the resulting digraph is strongly connected, meaning that all loads are intertwined in a single system that cannot be decomposed into islands. To build more accurate dynamic load models in the power system, we examine the mathematical characteristics of the power system variables. Toward this end, we demonstrate that the voltage magnitude, frequency, and voltage phase angle of xii Phasor Measurement Unit data exhibits long-range dependence. Our findings call for the development of a fractal modeling approach of the smart grid. Last but not least, we introduce a novel method to assess the distance to black- out or other instability of the smart grid. Based on the existence of long-range cor- relationinthePMUdata, weexhibitanincreaseinthefrequencyHurstexponent— quantified by Kendall’s tau rank correlation coefficient—before the blackout man- ifests itself in the grid powering down. Probably as the major result of this thesis, High Kendall’s tau of the frequency Hurst exponent is here proposed as an early- warning signal for blackout. xiii Introduction The smart grid concept was firstly introduced in 2005 [3] and the objective was to make the grid stronger, smarter and more secure. Simply speaking, the aim was to design a power grid that can cope with all new developments of renewable energy plants and to meet future energy demand such as the increase of deployed plug-in hybrid electric vehicle (PHEV). Aiming to design a smart grid that is efficient, reliable, resilient, stable, and secure requires both a fundamental understanding of its dynamics and an accurate mathematical framework. For instance, in order to have a stable power grid, the system should be able to regain the state of operating equilibrium after being subjected to a disturbance with most of the system variables remaining bounded [40]. The system stability can be classified into: voltage, frequency, and rotor angle stabilities. One fundamental challenge for ensuring voltage stability of the power grid is to deploy countermeasures to an imminent voltage collapse phenomenon in which the voltage drops to a low level due to a system disturbance [60]. The voltage collapse phenomenon has been studied extensively in the 80’s and 90’s. It was argued that the voltage collapse is either static or dynamic in nature; several researchers studied the voltage collapse from both aspects. In the static case, the voltage collapse was studied by looking at the feasibility of the load flow [28], the minimum singular value of the Jacobian matrix [62], and static xiv bifurcations of the load flow equations [41]. The static approach describes the load using the active and reactive powers, as shown in Eq. (0.1): P (V ) +jQ(V ) =k p (V ) np +jk q (V ) nq (0.1) where, k p and k q are the nominal active and reactive powers respectively. Voltage collapse was studied dynamically by investigating the interaction between the gen- erator and the load [16] and the interaction between the load and the On-Load Tap Changer (OLTC) [50]. The popular dynamic load model was introduced by Hill [33]byrepresentingtheactivepowerasafirst-ordernonlineardifferentialequation, as shown in Eq. (0.2): T p ˙ P d +P d =P s (V L ) +k p (V L ) ˙ V L (0.2) The dynamic model represents the load recovery with time constant (T p ) after a step change in the load voltage. The system frequency in this model is assumed to be constant. Even though many studies of the voltage collapse assumed no change in the system frequency during the voltage drop, there is still some debate on the need to investigate the dynamics of the frequency during this phenomenon. In [11], the author shows how the frequency dynamics can play a role in the voltage collapse, and explains how constant frequency assumption is not completely justifiable in case of a large time constant of the voltage. More recently, the research community and the government energy agency [52] advocate for the necessity of a better dynamic load modeling that can explain and describe the voltage and frequency variations over multiple time scales (short and long). xv This thesis follows up on those last footsteps, with the novel pitch of frac- tional dynamics and statistical PMU signal analysis with the ultimate objective of anticipating voltage collapse before it becomes catastrophic. xvi Chapter 1 Load Modeling 1.1 Introduction In the static analysis, active and reactive powers consumed by the loads are representedaspossiblynonintegerpowersoftheloadvoltagewithexponentsn p ,n q that depend on the load type as shown in Eq. (1.1), P L =K p V pv L , Q L =K q V qv L , (1.1) where K p and K q are constants. Classically, there are three types of static loads: constant power (p v = q v = 0), constant current (p v = q v = 1), and constant impedance (p v =q v = 2). Thedynamicloadmodel, whichwasproposedinHill[33],capturesthetransient load response to voltage step function. The model represents the active power as solution to a nonlinear differential equation. In this chapter, we use another load model, the Berg model [7], which in our opinion has been overlooked. Probably its most distinguishing feature is that it was derived in a real grid environment, which probably explains, in addition to noninteger exponents n p and n q of the load voltage, noninteger exponents of the frequency ω, as shown in Eq. (1.2), P (V L ,ω) =K p V pv L ω pω , Q(V L ,ω) =K q V qv L ω qω . (1.2) 1 Becauseofthisfrequencydependence, werefertothosemodelsas“inthestatic- dynamic gap." These models lend themselves to describing function “impedance" models, which are easily incorporated in multivariable feedback diagram and amenable to robust multivariable control techniques to reveal the possibility or the impossibility of voltage collapse. 1.2 Berg [7] Model versus Hill Model [33] The dynamical, transient-motivated model developed by Hill [33] reads T p ˙ P d +P d =P s (V L ) +k p (V L ) ˙ V L whereP d istheslowlytime-varyingactivepowerconsumedbytheloadandP s isthe steady-state component of the active power. P s (V L ) is certainly comparable with Model (1.1), valid in the specific cases of the filament lamp and the heater where the Berg modeling does not reveal frequency dependence. The major difference, however, is that ˙ P d and ˙ V L are the time derivatives of the active power and the voltage magnitude, resp., so that any describing function attempt to replace d/dt by jω sub , viz., P d (jω sub ,V L ) = αP s (V L ) +βjω sub k p (V L )V L 1 +jω sub T p would entail a subharmonic component ω sub , while the Berg model entail the har- monic (line frequency) ω. In a certain sense, the Berg and the Hill models complement each other, both having the similarity of a slow transient followed by recovery to steady-state (com- pare Berg [7, Fig. 8] with Hill [33, Fig. 1]), as shown in Fig.1.1. 2 The discrepancy, as shown in Fig.1.1, is that the Berg model ignores the tran- sient but captures a noninteger power of the frequency component in the steady- state whereas the Hill model captures the transient but ignores the frequency- dependence of the steady-state. (a) (b) Figure 1.1: Load Models: (a) Hill Model, (b) Berg Model Nonlinear load aggregation effect in the grid is a contributing factor to the non- integer exponent of the frequency, as the Berg model measured the dependency of the active power on voltage and frequency at the specific load, but in a real Scandinavian island grid. It could be argued that such exponents as w −2.1 in the induction motor could be the result of experimental inaccuracies, but such expo- nents as w 0.5 in the reduction furnace point to an inescapable reality, a fractional transfer function model of a high voltage transformer was already derived in [35]. The big difference between the Berg load model and the transformer model is that the fractional model of the transformer is mandated by matching the frequency response over a very large frequency sweep, exciting parasitic, distributed param- eter electromagnetic effects present in the wiring of the transformer. In our case, we deal with the frequency response over a very narrow frequency band around 50 cycles/sec. Under those circumstances, it is a bit difficult to see how distributed parameters modes could be excited, unless another hitherto unknown effect is 3 present. In Chapter 3, this other effect is conjectured to be some aggregation of the loads. 1.3 Static-Dynamic Gap Load Modeling Here, we fill the so-called “static-dynamic" gap by introducing describing func- tion models of such usual loads as reduction furnace, fluorescent lamps, etc. The static model was exlained in the introduction. The popular dynamical model of Hill [33] relies on a differential equation linking the instantaneous active power consumed by the load, the static power, and the voltage magnitude under sub- harmonic transients. On the other hand, the Berg model [7] adopted here ignores transients, but runs at the harmonic level and is hence amenable to robust multi- variabletechniques. Inacertainsense,wecompromisebetweenstaticanddynamic, with a purely harmonic model analytically extendable to include slightly damped behavior, and amenable to large scale multivariable analysis. In the Berg model [7], the active powerP L and the reactive powerQ L consumed arenonintegerexponentpowersoftheloadvoltageamplitudeV L andthefrequency w as shown in Eq. (1.3), P L =K p V pv L w pw , Q L =K q V qv L w qw , w> 0. (1.3) Intheabove,K p andK q areconstants. Theexponentsp v ,p w ,q v , andq w dependon the specific load. Table 1.1 shows the exponents for different loads (from [4, Table 6.1], [7,8]). 4 Load Type n p n q p f q f Filament lamp 1.6 0 0 0 Fluorescent lamp 1.2 3.0 -1.0 -2.8 Heater 2.00 0 0 0 Induction motor half load (HL) 0.2 1.60 1.5 -0.3 Induction motor full load (FL) 0.1 0.6 2.8 1.8 Reduction furnace 1.9 2.1 -0.5 0 Aluminum plant 1.8 2.2 -0.3 0.6 Table 1.1: Typical values of characteristic-load parameters The describing function of the load represents the equivalent gain from the load current I L to load voltage V L , Z L = V L I L . By multiplying the numerator and denominator by V L ∗ , we get Z L = V L V L ∗ I L V L ∗ = V 2 L S ∗ L . (1.4) The complex power is equal to S L =V L I L ∗ =P (V L ,w) +jQ(V L ,w). (1.5) Substituting Eq. (1.5) in Eq. (1.4) yields Z L = V 2 L P (V L ,w)−jQ(V L ,w) . (1.6) The describing function of different power grid loads can be derived by substituting the expression forP,Q of Eq. (1.3) in the describing function equation (1.6). The describing function of the load becomes: Z L = 1 K p V pv−2 L w pw −jK q V qv−2 L w qw , (1.7) 5 The describing functions of different types of loads are shown in Table 1.2. Load Type Describing Function Filament lamp 1 KpV −0.4 L −jKqV −2 L Fluorescent lamp 1 KpV −0.8 L w −1 −jKqV L w −2.8 Heater 1 Kp−jKqV −2 L Induction motor HL 1 KpV −1.8 L w 1.5 −jKqV −0.4 L w −0.3 Induction motor FL 1 KpV −1.9 L w 2.8 −jKqV −1.4 L w 1.8 Reduction furnace 1 KpV −0.1 L w −0.5 −jKqV 0.1 L Aluminum plant 1 KpV −0.2 L w −0.3 −jKqV 0.2 L w 0.6 Regulated aluminum plant 1 KpV 0.4 L w 0.4 −jKqV −0.4 L w 0.7 Table 1.2: Impedances of various loads For notational convenience, we switch to the admittance formulation Y L = 1/Z L = Lω p −jWω q . Next, we approximate Y L (V L ,ω) with a formal circuit- theoretic admittance, except for its amplitude dependency, Y L (V L ,ω)≈A(V L )× (jω) α −B(V L )× (jω) β , (1.8) where A, B are real-valued. The construction of the approximation is done as follows: First, write j α =a +jb,j β =c +jd. Then the approximant becomes A× (jω) α −B× (jω) β =aAω α −cBω β +j(bAω α −dBω β ). (1.9) The matching conditions then become aAω α −cBω β ≈Lω p , −bAω α +dBω β ≈Wω q . (1.10) 6 Matching the right and left-hand sides at ω = ω 0 = 2π50 rad/sec, along with the derivatives around ω = ω 0 , entails 4 conditions. This allows us to uniquely determine the 4 parameters A,α,B,β. The right-hand side of (1.8) has an obvious holomorphic extension, A× (σ + jω) α −B× (σ +jω) β that clearly satisfies the Cauchy-Riemann conditions. For ω = 0, the extension is obviously real and hence satisfies the “realness” condition of formal circuit theory [6]. 1.4 Power System Model The power system model consists of one generator source connected to a load (Z L )throughatransmissionadmittance(Y Line ), asshowninFig.1.2a. Thecircuit diagram in Fig. 1.2a can be redrawn to highlight the feedback structure with the line admittance (Y Line ) in the feedforward transfer function and the load (Z L ) is the feedback transfer function. This representation will be used through the whole analysis. The feedback model is shown in Fig. 1.2b. (a) (b) Figure 1.2: Power system model: (a) Circuit diagram, (b) Feedback model The generator will be modeled using a harmonic oscillator with a resonance frequencyw 0 = 377rad/s. Even though the generator model is a simplified one, it 7 still will give us some insight on how the voltage could collapse. A block diagram of the generator model is shown in Figure 1.3. Figure 1.3: Block Diagram of the harmonic oscillator The transmission line will be represented by an inductive admittance (Y Line ), assuming the purely inductive transmission line property can be justified as an approximation, that is, the resistance and the capacitive effects of the transmission line will be neglected. We will study different types of loads, one by one, to show conditions under which voltage collapse occurs. These models will be explained in details in the next chapter. 1.5 Summary The Berg model and the Hill model complement each other, on the common foundation of a transient followed by a recovery. It would be beneficial to develop a model that combines the Hill dynamics of the transient with the Berg frequency dependent static power. It would have the unique feature of a harmonic w and a subharmonic w sub frequency dependence. In Chapter 2, we integrate the load model with a simplified generator model, as shown in Fig.2.1, to investigate the voltage collapse phenomenon for different load types. 8 Chapter 2 Voltage Collapse in the Power System 2.1 Introduction The voltagecollapse has been explained using static and dynamic analysis. The static load has been used to study the voltage collapse by looking at the feasibility of the load flow [28], the minimum singular value of the Jacobian matrix [62], and static bifurcations of the load flow equations [41]. The dynamic load has been used to study the voltage collapse in the power system by looking at the interactions between the dynamics of the load and the generator or the load and the on-load tap changer as in [50] Voltage collapse is an intriguing phenomenon from a control perspective. Even though many theoretical “routes" to voltage collapse have already been proposed, it is not entirely clear what is really happening in the complex, large-scale envi- ronment of the power grid. It is our contention that voltage collapse does not only occur as a supply-demand imbalance, but could as well occur as a nonlinear phenomenon due to the complicated characteristics of the typical “messy" loads that the generators drive. The early suspicion that load characteristic is somehow the culprit in voltage collapse has led to deeper research on load modeling, but in a context where other phenomena (e.g., tap-changer effects) contribute to the collapse. In this chapter, we show that voltage collapse can occur—even with an 9 infinite bus, even with a simple one-generator, one-line, one-load configuration—as a phenomenon emanating from special effects in the aggregated load characteristic. This load characteristic is shown to be embedded in a hidden feedback structure that produces a nonlinear effect where, as the load voltage decreases, the damping of the oscillation increases. We call this phenomenon, the anti-Van der Pol effect. 2.2 Existence of a Voltage Collapse Solution Combining the generator model and the power system model, we will have the feedback model that is shown in Fig. 2.1. The model has a 2× 2 feedforward transfer matrix (G). The matrix G is an upper-triangular matrix, which has Y Line and w 0 2 on the diagonal. The feedback transfer matrix (F) is a 2× 2 diagonal matrix with−Z L and− 1 s 2 on the diagonal. Figure 2.1: Block diagram of the power system including the harmonic oscillator The output (y = [y 1 y 2 ] t ) can be written as function of the input (u = [u 1 u 2 ] t ), as follows: y = (I−GF ) −1 Gu, (2.1) 10 where G = Y Line Y Line 0 w 0 2 , F = −Z L 0 0 − 1 s 2 . It is important to note that Equation (2.1) should not be interpreted as a Laplace transform equation, but as an equation where y(jw) andu(jw) are “pha- sors" of the harmonic regime coswt in a sense formalized in [6]. We will further extendjw tos =σ +jw wheres is restricted to a small strip across the imaginary axis and around the harmonic frequency, where the analytic extension of the fre- quency response of the loads is guaranteed to exist, as will be proved in Sec. 2.2.1. Such an extension for stability analysis has been validated in [5]. 2.2.1 Analyticity Voltage collapse would mean showing that a purely harmonic solution to feed- backequationsundersomeV L wouldgotoasolutionoftheforme σt cos(wt),σ< 0, when V L drops below some nominal voltage. There is thus a need to validate the various describing functions for such signals, that is, do an analytical extension to w−jσ. The various functions are complicated and since we remain within some small neighborhood of 50 cycles/sec., analytical expansion in terms of small parameters is tempting. However, this can be accomplished only if such functions are analytic. The nonlinear impedance can clearly be written as Z L =R L +jX L and it is easily seen that bothR L andX L are of the formp(V L ,w)/q(V L ,w); some further manipulations allow us to rearrange terms so that both the numerator and the denominator of R L , X L involve positive (possibly noninteger) powers of w. Lemma 1. If p(V L ,w) and q(V L ,w) are real analytic in w around w 0 = 2π× 60, and if q(V L ,w 0 )6= 0, then p(V L ,w)/q(V L ,w) is real analytic at w 0 . 11 Proof. Define the function inv(x) = 1/x. It is real analytic at x 6= 0. Next, observe that 1/q(V L ,w) = (inv◦q)(V L ,w). Since the composition of two real analytic functions is real analytic [39, Prop. 1.4.2] (a corollary of the Faà di Bruno formula [39, Sec. 1.3]), it follows that 1/q(V L ,w) is real analytic. Finally, p(V L ,w)/q(V L ,w) is real analytic as the product of two real analytic functions having nonempty intersection of their domain of convergence [39, Prop. 1.1.7]. It remains to show that p(V L ,w) and q(V L ,w) are real analytic. Clearly, it suffices to show thatw r , wherer> 0 is a noninteger, possibly irrational exponent, is analytic. We simplify the exposition by assuming that r∈Q + . Lemma 2. w n/d , where n,d∈N is real analytic around w6= 0. Proof. Define the power function power r (x) = x r . For n∈ N, power n is clearly real analytic. Furthermore, power 1/d is real analytic by the real analytic inverse function theorem [39, Sec. 1.5]. Finally, observe that w n/d = (power 1/d ◦ power n ◦ 1)(w)andthelatterisrealanalyticbytherealanalyticpropertyofthecomposition of real analytic functions. 2.2.2 Voltage Collapse Condition The possibility for the power system in Fig. 1.2(a) to have oscillation is deter- mined by Eq. (2.1). The power system will have a solution for the output (y) when the input (u) is equal to zero iff the determinant of (I−GF ) is equal to zero. The determinant of (I−GF ) can be written as |I−GF| = (1 +Z L Y Line )(1 + w 0 2 s 2 ). (2.2) 12 Theorem 1. Consider the feedback interconnection of Fig. 2.1. Then, on the one hand, a harmonic solution cos(wt) always exists. On the other hand, a voltage collapsing solution e σt cos(wt) with σ< 0 exists if and only if 1+Z L (V L ,w−jσ)Y Line (w−jσ) = 0, for some σ< 0. Proof. In the particular case of one generator, one line, one load, the multivariable diagram of Fig. 2.1 yields y subject to collapse. By having the input u equal to zero in Eq. 2.1, existence of solutions is then given by det(I−GF ) = 0, which yields (1 +Z L Y L ) 1 + w 0 2 s 2 = 0. From here on, the result should be obvious. Assuming the load describing function is not equal to zero at any voltage (V L ) and frequency (w), the voltage collapse condition can be rewtitten as Y L (V L ,w−jσ) +Y Line (w−jσ) = 0, (2.3) where Y L =K p V pv−2 L ( w−jσ w 0 ) pw −jK q V qv−2 L ( w−jσ w 0 ) qw (2.4) and Y Line = K Line σ +jw , (2.5) where K Line is equal to Z base L . For simplification, we will substitute σ +jw by s, and multiply Eq.(2.3) by s. The voltage collapse condition becomes K p ( −j w 0 ) pw V pv−2 L (s) pw +1 −jK q ( −j w 0 ) qw V qv−2 L (s) qw +1 +K Line = 0. (2.6) 13 In the simple model of Fig. 2.1, the absence of “back-action" of the load to the generator results in two solutions only—the purely harmonic one imposed by the generator and the load-specific 1 +Z L Y Line = 0 solution. This architecture does not allow for a smooth transition from one to the other, so that some bifurcation needs to be invoked. But, as shown in Sec. 2.3, the message of this decoupled model is that the 1 +Z L Y Line = 0 solution reveals a behavior of the frequency concomitant with the voltage collapse that has been thus far derived in models where the swing equation of the generator appears to play a role [11]. Here, such behavior is created by the load. It thus appears that the load—properly modeled to include the frequency dependence as in the Berg model—might contribute much more to voltage collapse than has been assumed thus far. 2.3 The possibility of voltage collapse at differ- ent loads We will go over different groups of power system loads and identify the possi- bility of voltage collapse in each case. Firstly, we will go over synthetic loads with p v = q v and p w = q w . Then, we will use practical load models that were derived in [8]. Once we write the load impedance as a function of voltage and frequency, we will assume having an analytic extension of the imaginary axis interval (j(w 0 − δ),j(w 0 +δ)) in which Berg equations are valid to a neighboring strip comprising this interval and intersecting the right and left-half planes. That means that we will replace (w) by (w−jσ) in the equations of the load impedance (Z L ) and the line admittance (Y Line ). Recall that the damping ratio is ζ =−σ/w. 14 2.3.1 Synthetic loads with p v =q v ,p w =q w This group of loads have common voltage/frequency exponents for the active power and the reactive power. Even though these loads are more theoretical than realistic, however, they can help us to derive explicit expressions for the damping ratio (ζ) and frequency (w) as functions of the load voltage. Such solution can give us a bound on the voltage (p v =q v ) and frequency (p w =q w ) exponents of the loads that have a potential to collapse. Also, it can explain how other parameters, like active power coefficient (K p ), reactive power coefficient (K q ), and transmission line coefficient (K Line ), could support/mitigate voltage collapse. By having the coefficients of the active and reactive power equal to each other, Eq. (2.6) can be rewritten as, s =σ +jw = ( −k Line ( −j w 0 ) pw (K p −jK q ) ) 1 pw+1 | {z } α V L 2−p v p w + 1 | {z } β . (2.7) Thesystemwouldgotovoltagecollapseifthefollowingtwoconditionsaresatisfied: Condition 1:<(α)< 0 & =(α)> 0 Given that the damping ratio (ζ) and frequency w are positive during the voltage collapse,andsinceV L β isareal-valuedfunction,therealandimaginarycomponents of α should be negative and positive, resp. Condition 2: β < 0 A necessary condition for voltage collapse is that|σ| be inversely proportional to V L , which is satisfied when the exponent β of V L is negative. Our goal here is to find the ranges of voltage (p v ) and frequency (p w ) exponents that satisfy both of the conditions above. Using MATLAB, we found the region where Condition 1 is satisfied. We used the following numerical values:K p = 1, 15 K q = 0.328, and K Line = 1000. Condition 1 depends only on the frequency exponent (p w ) and it is satisfied when p w > 0.202. The second condition is satisfied if either (p v > 2 & p w >−1) or (p v < 2 & p w <−1). In Fig. 2.2, we plot the region where Condition 1 holds in red color and the region where Condition 2 in yellow color. The intersection of these regions (orange color) is the region that has potential for voltage collapse. Figure 2.2: Voltage collapse region To confirm the voltage collapse region, we plot σ and the frequency (w) versus the load voltage (V L ) at different load characteristics (p v i ,p w i ) chosen from the voltage collapse region of Fig. 2.2. We will sweep the voltage from 1 to 0.1. The curves in Fig. 2.4 for the three loads with exponents pairs ((2.5, 1), (3, 1.5), and (3.5, 2)) are showing increase in the damping ratio as the voltage decrease, which is a voltage collapse scenario. It is noteworthy that a similar pattern—abrupt decay of voltage concomitant with frequency increase–has already been observed in [11] using an approach that 16 involves the swing equation, as shown in Fig.2.3. Here, we show that a similar behavior can be obtained from the load characteristic only. (a) (b) (c) Figure 2.3: Voltage collapse phenomnon : (a) Voltage magnitude (V), (b) Phase angle of the generator bus (δ), (c) Frequency (ω) Figure 2.4: Sigma (σ) and frequency (w) for different synthetic loads 2.3.2 Berg load models Now, we will be more interested in applying the voltage collapse condition to the practical loads derived in Section 1.4 with load describing functions listed in Table 1.2. Since deriving an expression for the s = σ +jw solution is not easy 17 since the voltage collapse condition becomes more complicated with multiple non- integer exponents of s, we will use the ‘fsolve’ command in MATLAB to see how the damping ratio (ζ) and the frequency (w) depend on the load voltage (V L ). Based on the results, we can divide the loads into two groups: stable and unstable loads. Stable loads These loads are stable because they are not satisfying the voltage collapse condition. In case of Filament lamp, Heater, and Aluminum plant, the solution of the voltage collapse equation has negative frequency (w) values, which contradicts one of the assumptions (1.3). The solution of the voltage collapse equation does not exist in case of the Fluorescent lamp and the Reduction furnace. Induction motors (half/full load) have a solution for the voltage collapse equa- tion; however, the dependence of|σ| on the load voltage (V L ) is directly propor- tional, as shown in Figure 2.5. Figure 2.5: Sigma (σ) and frequency (w) for induction motor loads (half/full) 18 Unstable loads Regulated Aluminum plant load shows the possibility of having a voltage col- lapse as shown in Figure 2.6. This load is Aluminum Plant Load with voltage regulation that is accomplished by a transformer tap changer to maintain constant averageloadcurrent[7]. EventhoughtheAluminumplantbyitselfisastableload, the equivalent load for the load and the tap changer will make the load vulnerable to voltage collapse. Again the behavior of Fig. 2.6 has also been observed in [11] using the swing equation, while here this equation contribution does not appear to play a role. ‘ Figure 2.6: Sigma (σ) and frequency (w) for regulated aluminum plant The transmission line coefficient (K Line ), the active power coefficient (K p ), and the reactive power coefficient (K q ) have an impact on the proximity of the power system to voltage collapse. So, we will study how changing these variables could affect the damping ratio (ζ) at nominal load voltage (V L = 1). The closeness of the damping ratio (ζ) to zero at 1 p.u. means that a small bifurcation could throw the 19 voltage along the collapse route and is an indication of higher chance of voltage collapse. Figs. 2.7-2.9 show the dependency of the real partσ on transmission line coef- ficient, active power and reactive power coefficients. Since the transmission line coefficient (K Line ) is equal to Z base L Line and the maximum power transfer over trans- mission line is directly proportional to 1 L Line , the K Line is directly proportional to maximum power transfer. So, Figure 2.7 indicates that reducing the transmission line capacity could make the system closer to voltage collapse. Figure 2.7: The relationship between the transmission line coefficient (K Line ) and Sigma (σ) The impact of the active and reactive power coefficients on the voltage collapse is shown in Figs. 2.8-2.9. The higher these coefficients (K p andK q ) the higher the load demand, which could make the damping ratio (ζ) have a smaller value at 1 p.u., from where a small disturbance could take the system to voltage collapse. 20 Figure 2.8: The relationship between the active power coefficient (K p ) and Sigma (σ) Figure 2.9: The relationship between the reactive power coefficient (K q ) and Sigma (σ) 2.3.3 Correct initial conditions Since the results are essentially derived from the load, which, depending on the exponents in the model, would take the system from some initial conditions at 1 p.u. The correct initial conditions should be ζ(1p.u.) = 0 and w(1p.u.) = 21 2π60. The curves of Fig. 2.6 do not quite show the correct initial conditions; however, as shown by Fig. 2.7, manipulating the transmission line coefficient can bring us close to the correct initial condition ζ(1p.u.) = 0. The offset of the initial conditions could be created by complicated models extremely sensitive to exponents, themselves difficult to identify. However, despite this offset in initial conditions most probably created by numerical instability, the trend is correct as corroborated by [11]. 2.4 Summary We have shown a hitherto hidden route to voltage collapse. Its main feature is that it involves only a very simple model of the generator together with a model of the load realistically complicated by its strange frequency dependence, indicating that probably load characteristics have been overlooked in the general voltage collapse issue. Our finding of the possibility of increasing damping together with frequency increase/decrease as the voltage decreases is corroborated by the independent find- ing of [11], with the difference that our approach places more emphasis on the role of the load in this voltage collapse scenario. As far as the physically realistic frequency-dependent loads as modeled by the Berg paradigm are concerned, they can be classified into two groups: • Those loads for which the damping ratio ζ decreases as the load voltage decreases (induction motors). For those loads, voltage is stable, in a Van der Pol-like behavior. • Those loads that show an increase in damping ζ with decreasing voltage (regulated aluminum plant) together with a frequency disruption. These are 22 the “dangerous" loads, prone to voltage collapse, in an “anti-Van der Pol" behavior. The frequency-independent loads on the other hand do not create voltage col- lapse per the scheme unraveled here. What still needs further research is the transition from the harmonic regime (σ(1 p.u.) = 0) to the regime analyzed here characterized by increased damping under decreasing load voltage. Probably some couplings in Equations (2.1), (2.2) should give the key. 23 Chapter 3 Load Aggregation Effect in Power Grid 3.1 Introduction Load modeling is undoubtedly an important—still active—area of research in thepowergrid. Indeed, dependingontheiractiveandreactivepowerprofiles, loads could be the culprit in voltage collapse or other nuisances in the grid [11,55]. One issue that has never been made completely clear is what is the difference between, on the one hand, a model of a single isolated load and, on the other hand, a model of a load in a complex grid environment. Berg [7,8] is adamant that his load models represent the loads in the Scandinavian microgrid environment in which the experiments were done. Hill [33,50] on the other hand does “curve fitting" of the power disruption resulting from a voltage drop in a single load laboratory environment. A puzzling difference between the two models is that the Berg model involves the frequency in a fractional exponent in a narrow bandwidth around 50 cycles/sec, whereas a lingering issue in the Hill model is its lack of frequency dependence [33, p. 175], [52, p. 24]. The purpose is to offer evidence that his other effect is the aggregation of the loads. By “load aggregation," we mean that, after unraveling the hidden feedbacks in bus model, the resulting interconnected system is made up of one and only one 24 strongly connected component [9,10]. The latter means that under normal opera- tions loads are not “islanded." Quite to the contrary, every single load penetrates the whole grid. As an example [26], removing a load at a grid point in the Bay Area was observed by a micro-PMU 550 kilometers away in the Los Angeles area. 3.2 Berg load modeling and fractional dynamics Substituting s =σ +jω, the Laplace symbol, in the load describing function, Z L = 1 K p V pv−2 L ω pω −jK q V qv−2 L ω qω , (3.1) With this notation, the admittance takes the form Y L = I L V L =As α −Bs β . (3.2) Where A and B are real; otherwise the fractional differential models do not make sense. Reinterpreting the inverse Laplace transform of s α as the fractional differ- ential D α (∗) yields the fractional dynamics model: I L (t) =AD α ∗ V L (t)−BD β ∗ V L (t). (3.3) In our choice between the Riemann D α and the Caputo D α ∗ fractional differen- tial [18,29,42,43,47,48], we opted for the latter as its initial conditions involves only integer order derivatives of V L . Note that such a fractional order equation is known not to have periodic solutions [59]. It is believed that the fluctuations of the PMU variables observed in [53] betray their fractional dynamics. 25 3.3 Passivity of Berg loads In Eq. 3.3, the issue is whether the A and B coefficients are real when we go to the time domain fractional dynamics. It is closely related to that the admittance Y (V,jω +σ) is positive real, that is,<(Y (V,jω +σ))≥ 0 when σ≥ 0. Having the positive-real property means that the load has a passive real part that should not create any problems related to voltage stability. More precisely, this load would dissipate power only as expected (heater and filament lamp loads can obviously only dissipate power). However, when the sign of the real part changes, it means that the load could generate power under some conditions. Even though that might looks a little unusual as a load should dissipate power; this could be another sign of the load aggregation effect. In this section, we go over the positive-real property for two groups of loads, synthetic and berg loads. 3.3.1 Synthetic loads We study the positive-real property for synthetic loads (p v =q v and p ω =q ω ) by finding the condition under which the real part of the admittance stays positive whenσ≥ 0. Starting from the general form of the admittance describing function, Y (V,σ +jω) =K p V L pv−2 (ω−jσ) pω −jK q V L qv−2 (ω−jσ) qω (3.4) Assuming synthetic load with p v =q v and p ω =q ω , Y (V,σ +jω) =V L pv−2 (ω−jσ) pω (K p −jK q ) (3.5) The admittance can be written in a simplified form as, 26 Y (V,σ +jω) =K p ω pω V L pv−2 | {z } C (1−jσ/ω) pω (1−jK q /K p ) (3.6) Where C is positive and real. The synthetic load satisfies the positive-real propertyiftherealpartoftheadmittance (<(Y ))ispositivewhenσ≥0. Assuming σ, K p and K q are positive,<(Y ) is positive when p ω tan −1 (σ/ω) + tan −1 (K q /K p )≤π/2 (3.7) AssumingthatK q andK p valuesareclose, theresultindicatesthatthepositive- real property can be satisfied for small region around σ =0. 3.3.2 Berg loads In the section, we investigate the sign of the real part of the admittance for different Berg loads when voltage and sigma change . Then, we generate 2D and 3D plots for the real part of admittance as function of voltage and sigma. As shown in Fig. 3.1, we have 2D plots of the real part of the admittance as function ofsigma. The3DplotsoftherealpartoftheadmittanceisshowninFig.3.2. From the plots, it is clear that the sign of the real part of the admittance is positive for the region around σ=0. Asσ starts to move away from σ=0, the value of the real part of the admittance start to drop to negative values, as shown in Fig.3.2(d). Similar to the synthetic loads, the berg loads have positive-real property in the region around σ=0. 27 Figure 3.1: Real part of the admittance (Y) as function of sigma (σ) for different Berg loads: (a) Filament lamp (b) Fluorescent lamp (c) Induction motor- half load (d) Induction motor- full load (e) Reduction furnace (f) Aluminum Plant Figure 3.2: Real part of the admittance (Y) as function of voltage (V) and sigma (σ) for different Berg loads: (a) Filament lamp (b) Fluorescent lamp (c) Induction motor-halfload(d)Inductionmotor-fullload(e)Reductionfurnace(f)Aluminum Plant 28 3.4 Basic motif In the most basic power system (one generator with internal impedance, one line, one load), a feedback structure already appears in the light of its circuit theoretic model. The basic bus model of Fig. 3.3(a) has the 2-motif feedback structureofFig.3.3(b). ThefirstmotifreferstotheloopZ G Y 12 fedbythegenerator emfE G toproducethegeneratorvoltageV G , asshowninthetoppartofFig.3.3(b). The second motif is theZ L Y 12 loop that produces the load voltageV L . In the next section, we analyze those two basic motifs are repeated in a large scale architecture to yield the hidden feedback structure in the power grid. (a) (b) Figure 3.3: Simple 2-bus power system: (a) Circuit diagram (b) Feedback model 3.5 Basic motif interconnection Consider now the basic motif of Fig. 3.3(a) repeated in an interconnection pattern together with the basic feedback of Fig. 3.3(b) itself repeated in a matter consistentwiththebusmodel. Thelatterobviouslyleadstoacomplicatedfeedback interconnect. 29 D(U 1 ) D(U 2 ) D(U 3 ) D(U 4 ) Σ Σ Z Loads Y lines I L V L V l I l - E generators I other D(U) Figure 3.4: Hypothetical decomposition of a grid digraph into strongly connected components. Each strongly connected component (dotted circle) is itself a multi- variable version of the bottom motif of Fig. 3.3(b). Given a complicated feedforward-feedback path interconnection, the question as to whether closed-loop stability could be decided by examining closed-loop sta- bility of several “subsystems" was posed in a graph-theoretic setting by Callier and Desoer [9,10]. Given that transport of the various information/commodities occurs in a specific direction, the overall interconnection structure can be abstracted as a directed graph or digraphD = (V,E). A Strongly Connected Component (SCC), D(U),U⊆V, is a subset of vertices and oriented edges connectingU-vertices such that∀v i ,v j ∈D(U), there exists a directed path from v i tov j and a directed path from v j to v i , and furthermore the subset is maximal in the sense that adding a vertexv k 6∈U would destroy the strongly connected property. The SCCs are them- selves connected through the Interconnection Subsystem (IS), which, as shown in Fig. 3.4, has no oriented loops for otherwise some of theD(U i ) could not be max- imal. Then using the concept of essential set, each SCC component is partitioned into two parts: one that is the feedforward path and another that is the feedback part of the diagram, as illustrated in the dotted circle of Fig. 3.4. 30 We apply those concepts of [9,10] to the power grid. Instead of arguing in terms of active/reactive power flow, we argue in terms of such circuit theory commodities as generator electro-motive forces (E G ), line voltage drops (V ` ), line currents (I ` ), load voltages (V L ), and load currents (I L ). It is argued that the circuit theory interconnection has a great many feedbacks, hidden in the bus model, but revealed upon modeling the loads as impedances. Indeed, repeating the bottom motif of Fig. 3.3(b) yields the multivariable feedback shown in the dotted circle of Fig. 3.4, where the load impedances are in the feedforward path and the line admittances in the feedback path, although the reverse convention shown in the bottom motif Fig. 3.3(b) is fundamentally not different. Looking at the hypothetical decomposition of Fig. 3.4, it follows that, if a collapse occurs in one of theD(U i )’s, that part of the network upstream toD(U i ) won’t be affected. Now, assume there is adequate supply relative to the demand. Assume none of theD(U i )’s could go into the nonlinear feedback effect collapse as formulated in [55]. Then the whole interconnect would not go into collapse. Clearly, if such a decomposition as that of Fig. 3.4 were possible, checking for collapse would be immensely simplified and two loads in differentD(U i )’s won’t be aggregated in a feedback loop. The problem is that for any grid with the generator internal impedances prop- erly modeled, except in case of faults, the hidden feedback structure is strongly connected, cannot be decomposed in strongly connected component, and the var- ious loads are intertwined in an aggregating feedback. 31 3.6 Aggregation effect: one single Strongly con- nected component Before going through the detail, we provide an idea of the proof of Th. 2. Consider a 1-generator, 2-load interconnect, as shown in Fig. 3.5. Going from power flow to equivalent circuit theoretical formulation yields the multi-input, multi-output (MIMO) feedback of Fig. 3.6(a). By visual inspection, the latter interconnect looks like the MIMO feedback shown in the dotted circle of Fig. 3.4, with the extra difficulty that Fig. 3.6(a) involves generator internal impedance. The nodes of the feedback diagram of Fig. 3.6(a) are partitioned into voltage nodes and current nodes in the bipartite digraph shown in Fig. 3.6(b), obviously strongly connected. From this point on, it turns out that adding one generator will add another bipartite digraph, the two digraphs being connected in a way that allows us to go back and forth between the two of them. It turns out that this is a general fact; hence the hidden feedback graph is strongly connected. 3-bus power system In this case, we choose a 3-bus system with one generator at bus 2 feeding two loads at buses 1 and 3, as shown in Fig. 3.5. The feedback structure of this power system, depicted in Fig. 3.6(a), has three motifs, one for the generator and two for the loads. By looking at the feedback system, we can see that the generator is connected to the loadsZ L1 andZ L2 through line admittancesY 12 and Y 23 respectively. So, the system is strongly connected. 32 Figure 3.5: 3-bus power system: 1 generator, 2 loads The graph of this feedback system is a bipartite 2-color graph, as shown in Fig. 3.6(b). The green color is for the voltage nodes (V L1 ,V G , and V L2 ) and the orange color is for the current nodes (I L1 , I G , and I L2 ). (a) (b) Figure 3.6: 3-motif case (1 generator, 2 loads): (a) Feedback model, (b) Graph model Wecanfindthe“plant"transfermatrix ¯ Gandthefeedbackmatrix ¯ F bywriting the nodal equations for the power system, which yields V L1 V G V L2 | {z } η = −Z L1 0 0 0 −Z G 0 0 0 −Z L2 | {z } ¯ G I L1 I G I L2 | {z } e + 0 E G 0 | {z } V , (3.8) 33 I L1 I G I L2 | {z } e = Y 12 −Y 12 0 −Y 12 Y 12 +Y 23 −Y 23 0 −Y 23 Y 23 | {z } ¯ F V L1 V G V L2 | {z } η . (3.9) 6-bus power system In this case, we have a larger power system with 6 buses, 2 generators, and 4 loads as shown in Fig. 3.7 . The feedback model of the power system, as shown in Fig. 3.8(a), has 6 motifs, two for the generators and four for the loads 1, 2, 3, and 4. The feedback model is strongly connected. Figure 3.7: 6-bus power system The graph representation of this power system, depicted in Fig. 3.8, shows how generator and load motifs are connected. The graph is strongly connected, because the only way of decomposing the graph is by disconnecting two lines in the system, like lines 2-3 and 5-6. The plant matrix ¯ G and the feedback matrix ¯ F are derived below 34 V G1 V L1 V L2 V G2 V L4 V L3 | {z } η = −Z G1 0 0 0 0 0 0 −Z L1 0 0 0 0 0 0 −Z L2 0 0 0 0 0 0 −Z G2 0 0 0 0 0 0 −Z L4 0 0 0 0 0 0 −Z L3 | {z } ¯ G I G1 I L1 I L2 I G2 I L4 I L3 | {z } e + E G1 0 0 E G2 0 0 | {z } V , I G1 I L1 I L2 I G2 I L4 I L3 | {z } e = Y 12 +Y 16 −Y 12 0 0 0 −Y 16 −Y 12 Y 12 +Y 23 −Y 23 0 0 0 0 −Y 23 Y 23 +Y 43 −Y 43 0 0 0 0 −Y 43 Y 43 +Y 45 −Y 45 0 0 0 0 −Y 45 Y 45 +Y 56 −Y 56 −Y 16 0 0 0 −Y 56 Y 16 +Y 56 | {z } ¯ F V G1 V L1 V L2 V G2 V L4 V L3 | {z } η . 35 (a) (b) Figure 3.8: 6-motif case: (a) Feedback model, (b) Graph model IEEE 14-bus power system This 14-bus power system is one of the IEEE standard test cases which was provided online by [17] and depicted in Fig. 3.9 with its graph model shown in Fig. 3.10. 36 Figure 3.9: IEEE 14 bus power system Figure 3.10: Graph model in Small World configuration 37 3.7 The effect of contingencies on the power sys- tem graph In Section 3.6, we were able to show that any power system has a strongly con- nected feedback digraph. In this section, we show how the power system digraph canloseitsstrongconnectivityincaseofsingle-element (N−1)anddouble-element (N− 2) contingencies. There are different types of contingencies in the power system like, generator failure, transformer failure, transmission line disconnection, faults (single phase to ground, phase-to-phase, three-phase). Here, we will consider only transmission line tripping and three-phase fault and we will choose the 6-bus power system as our case study. 3.7.1 Single-element contingency In the graph shown in Fig. 3.8(b), each transmission line connects the two bipartite graphs by two edges in both directions. Disconnecting any transmission line randomly will not effect the strong connectivity of the graph. For example, if the line 5-6 is disconnected, the edgesV L3 -I L4 andV L4 -I L3 will disappear from the graph, but the remaining 2-3 line will ensure strong connectivity of the 2 bipartite graph interconnect. The weight of the edges connecting V L3 -I L3 and V L4 -I L4 will become Y 16 and Y 45 , respectively. 38 Figure 3.11: Graph model of 6-bus power system after disconnecting line 5-6 In case of a three-phase fault, the load impedance and the voltage at bus load become zero. Since the voltage at the load bus is zero, we can delete its vertex and all the edges connected to it. A fault at any of the loads will make the graph not strongly connected any more because there is no edge leaving the current vertex at the faulted bus. So, the graph will have two strongly connected components, one for the load current and the other one for the remaining vertices. For example, a three-phase fault at load 1 has a graph as the one shown in Fig.3.12. Figure 3.12: Graph model of 6-bus power system after a 3-phase fault at bus 2 3.7.2 Double-element contingency Disconnecting any two transmission lines from the power system will cause removing four edges from the graph, 2 for each line. This kind of contingencies 39 can make the system not strongly connected anymore and divide the graph into two smaller strongly connected ones. This kind of contingencies is very harmful to the power system, especially when the generation does not meet the demand in the new separated parts. If we disconnect the lines 2-3 and 5-6, the new graph will have 2 strongly connected components, each component with 1 generator and 2 loads, as shown in Fig.3.13. Figure 3.13: Graph model of 6-bus power system after disconnecting lines 2-3 and 5-6 In case of two three-phase faults, all the edges connected to the two voltage vertices will disappear as well as the two vertices themselves. The effect caused by these faults will be four strongly connected components, two of them have only one vertex and the other two have four vertices. If we choose for example load 1 and load 4 to have faults, the resulted graph model is shown in Fig.3.14. Figure 3.14: Graph model of 6-bus power system after faults at buses 2 and 5 40 The connected but not strongly connected cases of Figs. 3.12, 3.14 illustrate the necessity to restrict Theorem 2 to systems that have no faults. 3.8 Formal result In view of the above many examples, we can formulate the following major result: Theorem 2. Under the condition that the bus system has no signle-element or double contingencies like faults and line disconnections, if the bus model graph is connected, then it is strongly connected. Proof. The grid is composed of many “basic motifs," each motif being strongly connected. Given two motifs, say i and j, since the graph is connected, there are intermediate motifs, i + 1, i + 2, ..., j− 2, j− 1, connected with directed edges {i,i + 1},{i + 1,i + 2},...,{j− 1,j}, where{k,k + 1} denotes an edge without specification of its direction. Physically, each edge{k,k + 1} is an admittance Y flowing in one direction. But by the feedback created by this line admittance, there is another edge Y flowing in the other direction. Hence two consecutive motifs, say, k and k + 1 are connected by two edges: one k→ k + 1 and another one k + 1→k. Hence the graph is strongly connected. 3.9 Summary The main point is encapsulated in Theorem 2, saying that the hidden feedback structure of the power grid—essential for understanding subtle voltage collapse scenarios [55]—is strongly connected [9,10]. The latter is a robust property, as it survives some simple contingencies as shown in Sec. 3.7. This strongly connected 41 propertymeansthatthegreatmanyloadsinthepowergridareintertwined, sothat (P,Q) measurements made at one load do not reflect the nominal load behavior, but its behavior in the grid environment. At the limit of infinitely many loads, the nominal load will appear to have a distributed parameter frequency response, as the Berg model (1.3) shows. Parallel to this “first principles" approach, the fractal analysis of PMU sig- nals [53,54] points to the necessity of the Grünwald-Letnikov fractional dynam- ics [13] to explain the behavior of the signals. On the other hand, after some holomorphic extension of the Berg model, the inverse Laplace transform yields a Riemann-Caputo fractional dynamics [29], hence connecting the data driven and first principles approaches. Another corollary of the strongly connected component is that cascading failures are unlikely to be contained without human intervention. The Berg model [7] somehow misses the transient that leads to the perturbed steady state after voltage disruption, which is precisely taken care of by the Hill model [33]. Incorporating the Hill modeling of the transients in the Berg model would lead to the next generation of load models. Establishing the complete connection between the fractal PMU data-driven analysis and the first principles approach, relying on the Grünwald-Letnikov and the Riemann-Caputo fractional derivatives, respectively, calls for some reconcilia- tion between the two fractional dynamics approaches. 42 Chapter 4 Evidence of Long-Range Memory in Power System 4.1 Introduction A Phasor Measurement Unit (PMU) is a device that provides synchrophasor and system frequency estimates, as well as other optional information such as calculated megawatts (MW) and megavars (MVAR) [1]. The PMUs have been introducedinthe1990stosupportandovercomethedrawbacksoftheconventional Supervisory Control and Data Acquisition (SCADA) system. The drawbacks are related to the network security issues of the SCADA systems and the asynchronous data arrival due to the transmission delay and low sampling rate (one sample every 2-4 seconds). The PMUs provide secure data by having their own dedicated communication network. The data is measured at high sampling rate (30-50 samples per second). Moreover, the PMU resolves the issue of data delay by GPS time stamping of the data measurements. Aiming to build a more accurate dynamic load model, in this chapter, we perform an in-depth data driven statistical analysis of the power grid processes and identify a few fundamental mathematical characteristics that are essential for anticipating the voltage collapse phenomenon. Using real power system data, we demonstrate the existence of fractal and non-stationary behavior in the power grid 43 that justifies the need for capturing the frequency dependence as well as hints towards constructing long-range memory (dependence) models of the power grid dynamics. 4.2 Prior work and novel contribution The self-organized criticality has been discussed in the context of blackout impact size in power system [12],[14]. In [65], a long-range dependence in the elec- tricity prices has been measured. In the literature, the long-range dependence has not been applied to power system times series (voltage magnitude (V), frequency (f), and voltage phase angle (θ)). There are several methods to study the long- range dependence in time series. For stationary time series, we can use methods like, R/S, Variance, Absolute Moments, and Whittle.These methods do not give accurate results in case of non-stationary time series. However, Detrended Fluctu- ation Analysis (DFA) is a powerful and robust method to examine the long-range dependence in non-stationary time series. Our main contributions are: • We first investigate the statistical properties of the real power system data (volt- age magnitude (V), frequency (f), and phase angle (θ)) measured by several PMUs in Texas synchrophasor network. • Based on the observed non-stationarity, we employ the DFA method to calculate the scaling exponentα of time series (voltage magnitude (V), frequency (f), and phase angle (θ)) at different locations in the synchrophasor network. The scaling exponents indicate the type of mathematical models that should be constructed to describe the power system dynamics. 44 4.3 Statistical properties of the power system We will first describe the power system data we use in our analysis and then summarize our statistical analysis of several Phasor Measurement Units (PMUs) time series at various locations. 4.3.1 Description of the power system measurements We use real data from the power system obtained through several PMUs installed in the Texas synchrophasor network. Texas synchrophasor network has several PMUs distributed over several locations. The measured data of our inter- est are voltage magnitude (in p.u.), frequency (in Hz), and voltage phase angle (in degrees). These data are available online [30] for three PMUs at three locations: Baylor University, Harris Substation, and McDonald Observatory. There are four data sets of measurements at each location with one hour duration and at sampling rate of 30 samples/s. The data sets have been recorded on 5/25/2015 (6:00-7:00 PM and 10:00-11:00 PM), 05/27/15 (12:00-1:00 PM), and 05/30/15 (9:00-10:00 AM). 45 46 Figure 4.2: (a) The voltage magnitude at Baylor University on 05/27/15 (12:00- 1:00 PM) (b) The frequency at Baylor University on 05/27/15 (12:00-1:00 PM) (c) The voltage phase angle at Baylor University on 05/27/15 (12:00-1:00 PM) (d) The unwrapped phase angle (e) The empirical CDFs of the voltage magnitude at four intervals (900 seconds each) (f) The empirical CDFs of the frequency at four intervals (900 seconds each) (g) The empirical CDFs of the phase angle at four intervals (900 seconds each) (h) The empirical CDFs of the unwrapped phase angle at four intervals (900 seconds each) The plots of the voltage magnitude (V), frequency (f), and phase angle (θ) of data set 1 at Baylor university are shown in Figs. 4.2(a), 4.2(b), and 4.2(c) . We unwrapped the measurements of the phase angle, as shown in Fig. 4.2(d). The voltage magnitude (V) in Fig. 4.2(a) varies between 0.9632 and 0.9955 p.u., and the frequency (f) varies between 59.89 and 60.04 Hz, as shown in Fig. 4.2(b). There are noticeable events around the 40th minute of the voltage magnitude (V) plot and around 3rd minute in the frequency (f) plot. Due to the low sampling rate (30 samples/s) of the measured data, the frequency components of the data will be limited to 15 Hz and will not be useful to detect any high-frequency events. 4.3.2 Statistical analysis of PMU measurements Examining the stationarity of the PMU measurements is very important for choos- ing the most appropriate mathematical modeling framework. For instance, work- ing under the assumption that the processes are stationary and trying to construct compact mathematical models (with few parameters) can lead to misguiding con- clusions (e.g., applying a stationary method for quantifying long-range memory to a non-stationary time series can lead to misleading exponents). Consequently, in what follows, we investigate the stationarity of the voltage magnitude (V), 47 frequency (f), and phase angle (θ) by estimating their empirical cumulative dis- tribution functions (CDFs) over moving time intervals. A process is called second-order stationary if its mean and variance are constant over time, and the auto-covariance of the data does not depend on time. Moreover, the stationarity in the strict sense means the joint statistical distribution of any subset of the time series does not depend on time. To investigate the stationarity of the PMU data, we divide each data set into four intervals each with 900 seconds lengthandestimatetheircorrespondingempiricalCDFsforeachdataset. Identical empirical CDFs indicate that the processes are stationary; otherwise, they are non- stationary. The empirical CDFs of the voltage magnitude (V) in data set 3 are shown in Fig. 4.2(e). The empirical CDFs of the voltage magnitude (V) are not identical and empirical CDF of the first interval (black) is intersecting the empirical CDF of the second interval (blue). Also, the empirical CDFs of the frequency (f), shown in Fig. 4.2(f), are not identical and the empirical CDF of the first interval (black) is crossing the empirical CDFs of the three other intervals (blue, green, and red). The empirical CDFs of the phase angle (θ) and unwrapped phase angle, shown in Figs. 4.2(g) and 4.2(h), are crossing each other several times. So, the voltage magnitude (V), frequency (f), and phase angle (θ) in data set 3 are non-stationary. The non- stationarity has been confirmed for all the locations in the Texas synchrophasor network. 48 4.4 Existence of long-range dependence In this section, we first explain the DFA method, and then, study the long-range dependence in the time series of the voltage magnitude, frequency, and voltage phase angle. 4.4.1 Detrended Fluctuation Analysis (DFA) The measured data of the voltage magnitude (V), frequency (f), and phase angle (θ) in the power system are non-stationary, therefore we will not be able to use methods like, Variance, Absolute Moments, and R/S to confirm the long- dependency in the data. That is because these methods were derived based on the assumption that the data are stationary. The DFA method was introduced in [45] to study the mosaic organization of DNA nucleotides. This method has been proven as a robust method to show existence of long-dependency in time series. This method has been applied successfully in several non-stationary data like, heartbeat fluctuation [46], daily temperature [38], and wind speed [36]. For a data set y with N data points, the DFA analysis can be summarized in four steps: 1. Subtract the data set average y avg from each data point y(i) and integrate the data set using Eq. (4.1), y int (k) = k X i=1 (y(i)−y avg ) (4.1) 2. Divide the integrated data set into equal-sized boxes n and find the linear least square line y n inside each box. Then, subtract the least square linear fitting 49 y n (k) from the integrated data y int (k) to generate the detrended data y d , as shown in Eq. (4.2), y d (k) =y int (k)−y n (k) (4.2) 3. Find the root mean square (rms) fluctuation of the detrended data y d using Eq. (4.3), F (n) = v u u t 1 N N X k=1 (y d (k)) 2 (4.3) 4. The second and third steps are repeated at different box sizes n, and we find the least square line of the log-log plot of F (n) versus n. The slope of the least square linear fitting is called the scaling exponent α [51]. A linear relationship between log 10 (F (n)) and log 10 (n) means F (n)∝ n α . Based on the value of α, the long-range dependence in the data can be determined. A scaling exponent α between 0 and 0.5 means that the data are anti-correlated. The measured data is random if α is equal to 0.5. For α between 0.5 and 1, a long-range power-law dependence exists in the data. Finally, α > 1 indicates long-range dependence, but not in the power-law form. If the log-log plot of F (n) versusn has one linear trend, then the data is called mono-fractal. Some data sets exhibit a change in the slope of the least square line at a transition point. This phenomenon is called crossover and the data are considered to possess multi-fractal characteristics. 4.4.2 Analysis of the power system data using DFA As part of the DFA analysis, we have to find the root mean square (rms) of fluctuation at different box sizesn. We choose the box size (n) to change from 100 data points to 10,000 data points with logarithmic step size equal to 10 1 10 . 50 51 Figure 4.3: Log-log plots of the rms fluctuation function F (n) versus the box size n of V, f, and θ of different data sets: (a) Data set 1 at Baylor University (b) Data set 2 at Baylor University (c) Data set 3 at Baylor University (d) Data set 4 at Baylor University (e) Data set 1 at Harris Substation (f) Data set 2 at Harris Substation (g) Data set 3 at Harris Substation (h) Data set 4 at Harris Substation A logarithmic spacing is preferred because choosing a linear spacing will make the weight of the log-log plot higher by moving from 100 to 10,000. The range of the log 10(n) is from 2 to 4 with 0.1 step size. We applied the DFA method on four data sets of voltage magnitude (V), frequency (f), and phase angle (θ) at three different locations. So, our idea is to run the DFA analysis on data sets that are distant in time and space to make sure that results are robust and comparable. DFA analysis of the voltage magnitude The DFA analysis of the voltage magnitude (V) at the three different locations are shown in Table. 4.1. The scaling exponents of the four data sets at Baylor location are between 0.91 and 1.11. These values are close to the scaling exponent of the white noise (1/f) which has scaling exponent equal to 1. At Harris location, the scaling exponents are in the range of [0.81, 0.92] and they are also not far from the scaling exponent 1 (white noise). The McDonald location has scaling exponents between 1.30 and 1.37 which are clearly higher than the other two locations. The McDonald PMU was installed in the Texas synchrophasor network particularly at the McDonald Observatory because it is among several wind power plants in that part of the network [31]. Our explanation of the higher scaling exponents at this location could be related to the voltage fluctuation associated with the connected wind power plants in the area [15][44]. Figs. 4.3 (a-h) show the rms fluctuation function of the voltage magnitude (V) at Baylor and Harris using red-color data 52 points. Finally, since the scaling exponent at three locations varies between 0.81 and 1.37, long-range dependence does exist in the voltage magnitude (V). Data Baylor Harris McDonald Set V f θ V f θ V f θ #1 1.11 1.54 0.71 0.92 1.54 0.75 1.32 1.54 0.74 #2 1.11 1.53 0.66 0.81 1.53 0.63 1.30 1.53 0.64 #3 1.05 1.45 0.67 0.91 1.45 0.76 1.37 1.45 0.73 #4 0.91 1.49 0.63 0.89 1.49 0.64 1.32 1.49 0.64 Table 4.1: Scaling exponents of voltage magnitude, frequency, and phase angle for four data sets at three locations: Baylor University, Harris Substation, and McDonald Observatory DFA analysis of the frequency The scaling exponents of the frequency (f) of the four data sets at each location are shown in Table. 4.1. The three locations (Baylor, Harris, and McDonald) have similar scaling exponents between 1.45 and 1.54. All the data sets of the frequency (f)exhibitlong-rangedependence, butnotinapower-lowform. Thefrequency(f) has scaling exponents similar to Brownian noise which has 1.5 scaling exponent. The plots of the rms fluctuation function of the frequency (f) are shown in blue color in Figs. 4.3 (a-h) for Baylor and Harris locations. DFA analysis of the phase angle The scaling exponents of the phase angle (θ) data are shown in Table. 4.1. The scaling exponents of the phase angle (θ) varies between 0.63 and 0.76. Since the scaling exponents of the phase angle (θ) are between 0.5 and 1, that means the phase angle (θ) data have long-range power-law dependence. Figs. 4.3 (a-h) show the rms fluctuation function of the phase angle (θ) in green color. 53 Figure 4.4: Log-log plots of the rms fluctuation function F (n) versus the box size n of original and surrogate data at: (a) Baylor University (b) Harris Substation (c) McDonald Observatory 4.4.3 DFA analysis of the surrogate power system data To show the robustness of our result, we generated a randomly shuffled version of the PMU data (surrogate data), which has an identical CDF compared to the original data. Then, we run the DFA analysis on both the original and surrogate datatoidentifyanychangeinthecorrelation. Ifweobserveadecreaseinthescaling exponentα, the fractal behavior is due to the temporal structure of the measured data. Thistemporalstructureplaysafundamentalruleinconstructinganaccurate mathematical model. We generated log-log plots of the rms fluctuation function F (n) versus n for the original and surrogate on three data sets at three different locations, as shown in Fig. 4.4. The results show that the scaling exponents of voltage magnitude (V), frequency (f), and phase angle (θ) have changed to around 0.5 after the random shuffling of the data. Of note, time series characterized by 0.5 scaling exponent can be modeled through first order differential equations. However, theobservedtemporalstructurecharacterizedbyscalingexponenthigher 54 than 0.5 implies that the dynamics of the system variables should be characterized through fractional order differential equations. 4.5 Summary In this Chapter, we performed a statistical analysis of the PMU measurements, which demonstrates the existence of mono-fractality and non-stationarity as two main mathematical characteristics. In Chapter 5, these findings have recently been corroborated by European PMU data analysis; furthermore, the same analysis on India voltage collapse data seems to indicate that monitoring the PMU fractal dimension is able to anticipate imminent voltage collapse. The next challenge is to develop models able to reproduce the long range dependency of frequency, voltage and phase angle PMU data. Fractional dynamics models would reproduce such long range dependency. It is believed that the fractional powers of the frequency in the Berg model, properly re-interpreted in the time domain, would produce such models, as it has already been able to anticipate hitherto unknown voltage collapse scenarios [55]. 55 Chapter 5 Kendall’s Tau of Frequency Hurst Exponent as Blackout Proximity Margin 5.1 Introduction Due to the high sampling rate of the data measured by PMUs, extracting real- time useful information in a timely manner could be a challenge. The detection of events in the power grid using PMUs has been an active area of research [23,24,58]. Voltage collapse and more generally power system blackout—either accidental or malicious—are among the most severe events that cause power loss over a wide area of the power system. Anticipating such events is a high priority in smart grid research. In 2003, one of the largest blackouts in US history hit the Northeastern and Mid- western parts of the United States, and the Canadian province of Ontario. The blackout left 55 million without electricity with total economic cost between $7 and $10 billion. Two large consecutive blackouts occurred in the Northern part of the Indian power grid on July 30, 2012 and July 31, 2012. These two black- outs affected 300 million and 600 million, respectively. Unfortunately, the classical loading margin [57] is static and while it has been argued that a dynamic approach 56 is irrelevant [20], here we show that fractal dynamics is relevant and could provide another early warning of blackout. Wefirststudythescalingpropertiesof2079mediumvoltage(12kV)PMUdatasets (100,000 samples each) collected over four months from the École Polytechnique Fédérale de Lausanne (EPFL) campus network. We apply the DFA method on the voltage magnitude, frequency, and phase angle data sets to show the strong consistency of the scaling exponents of 120V and 69 kV PMU data collected from the Texas Synchrophasor Network [53] and the EPFL campus network. From our earlier work [53], it appears that scaling properties of the frequency do not depend on the voltage. Secondly, we introduce a new method based on the change of fractal characteristics of the frequency data before the power system goes to blackout. We compare the Hurst exponents of frequency data collected from the EPFL campus network and frequency data collected before the 2012 Indian blackout. The increase in the Hurst exponent of the frequency time series can be used as a new early-warning signal of the proximity of the power system to a blackout. 5.2 Related work and novel contribution Complexsystemsusuallyhavecriticalthresholdsbeforetheygointosuddenchange from one state to another. Because of the complex structure of these systems, it is certainly not easy to build models that describe their dynamics accurately. However, such systems show a critical slowing down phenomenon as early-warning signal prior to the critical transition [66]. Critical slowing down means that the system needs more time to recover from a perturbation before critical transition. 57 The critical slowing down manifests itself as a sudden increase in the first auto- regressive (AR(1)) coefficient before the occurrence of the transition [34]. Critical slowing down has been used as an early-warning signal for critical tran- sition in climate [66], environment [22], ecosystems [25], and power system [19]. In [19], the slowing down has been investigated by calculating the AR(1) coefficient of the frequency data collected before the Western Interconnect Blackout in 1996. Our contribution consists in using another frequency parameter as an indication of an imminent blackout. Specifically, our contributions are: • First, we show mono-fractality of the PMU data (V,f, andθ) collected in EPFL campus during normal operation. We run the DFA over a PMU big data set collected over four months (January, April, June, and October) in 2015. • Second, we introduce a new measure of the proximity to blackout in the power system by comparing the Kendall’s tau of the Hurst exponent of the frequency data collected from the EPFL campus (normal condition) and before the 2012 Indian blackouts. 5.3 Power system data In this section, we give descriptions of the PMU data (V,f, andθ) collected in the EPFLcampusandreviewhowtheHurstexponentofthesePMUdataiscomputed. 5.3.1 Description of the PMU data The EPFL campus network has five PMUs installed throughout the campus to collect several data measurement of the power system variables. These variables are voltage magnitude (V), frequency (f), phase angle (θ), active power (P), and 58 reactive power (Q). The PMU data measurements for several months in 2015 are available online [49]. The PMUs are installed at the medium voltage side with nominal voltage 12 KV and nominal frequency 50 Hz. The three time series of voltage magnitude, frequency, phase angle measured on 1/15/2015 (8:00-8:33 AM) are shown in Figs.5.1 (a-c). Figure 5.1: Analysis of PMU#4 data recorded at EPFL on 01/15/2015 (8:00- 8:33 AM) during normal operation. First row: raw data recorded; second row: corresponding log-log plot of rms fluctuation function versus the box size n. First (red) column: voltage magnitude; second (blue) column: frequency; third (green) column: phase angle In our analysis, we have chosen four months (January, April, June, and October) of 2015 to represent the different seasons in the year. Then, we picked four time series (100,000 samples each) per day of each of the power system variables (V, f, and θ). Two of the four time series were measured during day time (8:00-10:00 AM) and the other two measured during night time (8:00-10:00 PM). Based on the 59 availability of data, we will analyze 693 time series of each of voltage magnitude (V), frequency (f), and phase angle (θ) with a total of 207.9 million data samples. 5.3.2 PMU data fractality In [53], we have shown the existence of non-stationarity and mono-fractality of the PMU data sets collected from the Texas Synchrophasor Network. Here, we analyze a larger number of PMU data from the EPFL campus to compare the results with PMU data from Texas. Also, the large size of the data gives some confidence in the calculated scaling exponents of the data. Using the DFA method, we found the scaling exponents of the three time series shown in Figs. 5.1 (a-c). The plots of the rms fluctuation function (F (n)) versus the window size (n) for each of the three power system variables are shown in Figs. 5.1 (d-f). The scaling exponent of the voltage magnitude (V) in Fig. 5.1 (d) is 1.20. The time series of frequency (f) is similar to Brownian noise with scaling exponent equal to 1.55, as shown in Fig.5.1 (e). The scaling exponent of the phase angle time series is around 1.27 as shown in Fig. 5.1 (f). The Hurst exponent histograms of the time series collected from EPFL are shown in Figs. 5.2 (a-c). The Hurst exponent histograms of the voltage magnitude, frequency, and phase angle are shown in red, blue, and green, respectively. Each of the Hurst exponent histograms represents the Hurst exponents of one of the power system variables (V, f, or θ) collected during the four months (January, April, June, or October). The Hurst exponent histogram of voltage magnitude, shown in Fig. 5.2 (a), has mean 1.23 with standard deviation 0.11. The histogram of the frequency time series has a mean Hurst exponent 1.51 with standard deviation 0.05. Finally, the phase angle Hurst exponent has mean and standard deviation 1.21 and 0.08, respectively. 60 Figure 5.2: Hurst exponent histograms of EPFL PMU#4 data during all 4 seasons (January, April, June, and October): (a) Voltage magnitude (b) Frequency (c) Phase angle 5.4 Kendall’s tau as Blackout Proximity Margin The power blackout usually starts as instability in the voltage, frequency, or phase angle that leads to a major blackout. The power system blackout is traditionally explained using bifurcation theory [21], where the system at bifurcation point move from the stable region to the unstable one. To test the possibility of existence of critical slowing down phenomena in the power system, in Sec. 5.4.1, we investigate the change in the AR(1) coefficient (short- range correlation) in the frequency time series before the 2012 Indian blackout. In Sec. 5.4.2, we introduce our new method to anticipate an imminent power system blackout based on the change of the Hurst exponent (long-range correlation) before the blackout. Then, in Sec. 5.4.3, we investigate the change in AR(1) coefficient and Hurst exponent over large data set of frequency time series collected from the EPFL campus during normal conditions. 61 5.4.1 Kendall’s tau of AR(1) The frequency time series, collected before the 2012 Indian blackout, is shown in blue color in Fig. 5.3(a). The length of the time series is 167,600 samples (∼56 minutes) and it is measured at sampling rate of 50 samples/second. The frequency time series is non-stationary [53]; therefore, we should remove the trends in the time series before calculating the AR(1) coefficient. Figure 5.3: (a) Raw (blue) and Gaussian kernel smoothed (red) frequency time seriesbefore2012Indianblackout(b)Autocorrelationcoefficientatlag1ofresidual differences between raw and smoothed frequency time series (c) Hurst exponent of the raw frequency time series We remove the trends in the frequency time series using Gaussian convolution kernel[64]. Wefirstfindthesmoothedfunctionofthetimeseriesandthensubtract the smoothed function from the original time series to get the residual. Choosing the bandwidth (∼ 2.7σ) of the Gaussian kernel correctly is a critical step to avoid underfitting (large bandwidth) and overfitting (small bandwidth). The bandwidth thathasbestsmoothingforthefrequencytimeseriesis18,000samples(6minutes). Using the command ksmooth in R software [61], we calculated the Gaussian kernel smoothed function of the frequency time series before the blackout as shown in red color (over blue frequency color) in Fig. 5.3(a). We have chosen the normal kernel with 18,000 samples bandwidth. 62 Now, we find the change in the AR(1) in the residual time series by calculat- ing the AR(1) coefficient over a window (grey shaded area) of 100,000 samples (∼33 minutes). Then, we move the window 100 samples (2 seconds) to the right (toward the blackout) and calculate the AR(1) coefficient over the next window. Fig. 5.3(b) shows the change in the AR(1) coefficient of the residuals time series before the blackout. The AR(1) coefficient was calculated using ar.ols command in R software. The AR(1) coefficient before the blackout shows an increase from 0.89 to 0.97. The increase starts around 33 minutes before the blackout, which could be an early- warning for power system blackout. The increase in the AR(1) coefficient can be quantified using Kendall’s tau [37] [63] [2] [56] [32]. Kendall’s tau is a rank correlation coefficient that is used to measure the ordinal association between two quantities. Assuming that we have n pairs of x and y ((x 1 ,y 1 ), (x 2 ,y 2 ),..., (x n ,y n )). Kendall’s tau is defined as: τ = # of concordant pairs− # of discordant pairs n(n− 1)/2 (5.1) The pair is concordant if x i > x j & y i > y j or x i < x j & y i < y j . On the other hand, the pair is discordant if x i >x j & y i <y j or x i <x j & y i >y j . The range of Kendall’s tau is between -1 and 1. TheKendalltauoftheAR(1)coefficientcouldgoashighas0.92closetoablackout. With 685 data points, the Null Hypothesis of no trend in the AR(1) is rejected withp = 2.2×10 −16 . The very high confidence level means that the AR(1) increase is symptomatic of a dynamical shift in the grid. However, the observed increase in the autocorrelation could not be exclusively happening before the power system blackout. Indeed, in Sec. 5.4.3, we will show that such high confidence changes in 63 the AR(1) coefficient are indeed happening over large data set of frequency time series collected from the EPFL campus network under normal conditions. 5.4.2 Kendall’s tau of the Hurst exponent as dynamic proximity margin to blackout In Sec. 5.3.2, we have shown the existence of long-range memory in the frequency time series in the power system by analyzing a large number of data sets collected from the EPFL campus network. Our novel method to predict the power system blackout will be based on the change in long-range correlation (Hurst exponent) of frequency time series instead of the short range correlation at lag 1. In Fig. 5.3(a), we have the frequency time series before the 2012 Indian blackouts for 56 minutes. We study the change in Hurst exponent of a moving window of length 110,000 samples (∼37 minutes) and a shift of 900 samples (18 seconds). We calculated the Hurst exponent of each window using the Detrended Fluctuation Analysis (DFA) method, which is robust against the presence of trends and non- stationarity in the time series. Using dfa command in R software, the Hurst expo- nents of the frequency before the blackout were calculated as shown in Fig. 5.3(c). To calculate the Hurst exponent, we have used 21 window sizes with range from 100 to 10,000 samples. The Hurst exponent before the blackout increases from 1.55 to 1.72. Kendall’s tau of the Hurst exponent is 0.86. With 65 data points, the Null Hypothesis of no trend is rejected with p = 2.2× 10 −16 , a very high confidence level. This increase in the Hurst exponent before the blackout can be a sign of the proximity of the power system to blackout. Next, we need to show that the increase in the AR(1) coefficient and Hurst expo- nent before the 2012 Indian blackout is unique and can be used as reliable measures 64 for the proximity of the power system to blackout. This goal will be achieved by calculating the change in the AR(1) coefficient and Hurst exponent of several fre- quency data sets collected from EPFL campus under normal conditions. 5.4.3 Kendall’s tau of AR(1) coefficient versus Kendall’s tau of the Hurst exponent under normal conditions In this section, we use the Kendall’s tau to quantify the change in the AR(1) coef- ficient and Hurst exponent of 230 frequency time series (180,000 samples each). These time series were collected during normal conditions from the EPFL campus. Kendall’s tau close to 1 means that AR(1) coefficient or Hurst exponent has con- sistently increasing trend. However, Kendall’s tau close to -1 indicates consistently decreasing trend. Then, we will compare the Kendall’s tau of the frequency time series before blackout with the ones that are collected during normal conditions. The AR(1) coefficient of 230 frequency data set is calculated over 100,000 samples window with 100 samples shift between consecutive windows. The Kendall’s tau distributions of frequency AR(1) coefficient during each of the four months are shown in Fig. 5.4(a). The histogram of Kendall’s tau of the frequency AR(1) coefficient for all the month is shown in Fig. 5.4(b). The histogram has -0.08 mean (black line) and 0.57 standard deviation (orange line) with range between -0.97 and 0.96. The histogram is very close to a uniform distribution over the range of the Kendall’s tau from -1 to 1. Since the Kendall’s tau of the frequency AR(1) coefficient before the Indian blackout is 0.92 (red line) and inside the range of the histogram of the EPFL data, there is no significant difference between blackout and normal data. The Hurst exponent of each frequency data set is calculated over a 100,000 samples window with 1,000 samples shift between consecutive windows. Then, we calculate 65 Figure 5.4: (a) Box-plots of Kendall’s tau of frequency AR(1) coefficient for data collected in EPFL during January, April, June, and October (b) Histogram of Kendall’stauoffrequencyAR(1)coefficientduringthesefourmonths(c)Box-plots of the Kendall’s tau of frequency Hurst exponent for data collected in EPFL during January, April, June, and October (d) Histogram of Kendall’s tau of frequency Hurst exponent during these four months. Due to blackout data scarcity, a red histogram and a test for random draws from different distributions couldn’t be obtained. the Kendall’s tau of the Hurst exponents for each of 230 frequency data sets. The distributions of the Kendall’s tau of the frequency Hurst exponent for each of the four months are shown in Fig. 5.4(c). The histogram of the Kendall’s tau for all the frequency data sets is shown in Fig. 5.4(d). The mean of the histogram is -0.04 (black line) and the standard deviation is 0.39 (orange color). It is clear that the histogram of the Kendall’s tau of frequency Hurst exponent is centered around 0 withrangebetween−0.77and 0.77. TheKendall’stauoffrequencyHurstexponent before the Indian blackout is 0.86 (red line) and it is outside the histogram range 66 of the normal EPFL data. That means there is a difference between the blackout and normal data and hence the Kendall’s tau of frequency Hurst exponent is a good measure of the proximity to blackout. 5.5 Summary We first investigated the long-range memory in a large PMU data sets from the EPFL campus network. All the PMU data (V,f, andθ) showed a long-range cor- relation with average Hurst exponents approximately 1.23, 1.51, and 1.21, respec- tively. Next, we studied the existence of critical slowing down phenomenon in the frequency time series before the power system blackout. We provided an evidence that the increase in the autocorrelation coefficient at lag 1 before the blackout is not specific to blackout; however, the increase in the Hurst exponent of the frequency, more specifically having a high Kendall’s tau of the frequency Hurst exponent, can be a better early-warning signal for power system blackout. Should the Hurst exponent increase, the pressing question is whether it is global (catastrophic) or area-local (manageable) and what has caused the increase. Assuming that the area of Hurst exponent increase and the root cause are identi- fied, the action would be to act as quickly as possible using FACTS, at the primary (millisecond)layerFig. 5of[27], consistentlywiththefastsamplingrateofPMUs. Should the problem not be rectified at this layer, we would attempt to rectify it by reactive power management at the secondary (seconds) layer. Finally, if the prob- lem is not yet resolved at the secondary layer, the last attempt would be to resolve it at the tertiary (minutes) layer, by for example islanding, or load shedding, which unfortunately would affect some consumers. The overall process is depicted in the tree of Fig. 5.5. 67 Figure 5.5: Decision tree acting on the multi-layer decomposition of Fig. 5 of [27] as response to increase of Hurst exponent 68 Chapter 6 Conclusion It is our hope that this thesis has provided yet another hitherto hidden "road to voltage collapse." Theoretically, part of its foundation rests on fractional dynamics and the other part on graph theory–the load aggregation effect. Next to its theo- retical foundation, it also has a more practical aspect: the objective analysis of the data, showing long range dependence. While intuitively the latter appears to be related to the theoretical foundation, the bridge between fractional dynamics and long range dependency of the data still remains to be reinforced. Also in deed to be further explored, one will mention the connection among the various fractional dif- ferential models and the Auto-Regressive Fractionally Integrated Moving-Average (ARFIMA) model. Even though some of these mathematical development still need further exploration, the reward already at this stage appears enormous–the possibility to anticipate blackout by statistical signal analysis before the blackout develops to a major catastrophic event. Finally, it is hoped that PMU signals under normal conditions and under Cyber attacks will show discrepancies that can be caught by the same fractal PMU signal analysis that has allowed us to antici- pate blackout. This conjecture is supported by some Cyber atatck scheme aiming at creating voltage collapse. Regarding your dissertation, I believe an important issue is whether the A, B coefficient when you go to the time domain fractional dynamics are real. Closely related, is the Y (V,jω +σ) admittance positive real, that is,<(Y (V,jω +σ))≥ 0forσ≥ 0? 69 Reference List [1] IEEE Standard for Synchrophasor Data Transfer for Power Systems. IEEE Std C37.118.2-2011 (Revision of IEEE Std C37.118-2005), pages 1–53, 2011. [2] H. Abdi. The Kendall rank correlation coefficient. In N. Salkind, editor, Ency- clopedia of measurements and Statistics. Sage, Thousand Oaks, CA, 2007. [3] S. M. Amin and B. F. Wollenberg. Toward a smart grid: power delivery for the 21st century. IEEE Power and Energy Magazine, 3:34–41, September 2005. [4] J. Arrillaga and C. P. Arnold. Computer Analysis of Power Systems. John Wiley and Sons, Chichester, New York, Brisbane, Toronto, Singapore, 1990. [5] D. P. Atherton. Stability of Nonlinear Systems. Research Studies Press, New York, 1981. [6] V. Belevitch. Classical Network Theory. Holden-Day, San Francisco, 1968. [7] G. J. Berg. System and load behavior following loss of generation— experimental results and evaluation. Proc. IEE, 119(10):1483–1486, October 1972. [8] G. J. Berg. Power system load representation. Proc. IEE, 120(3):344–348, March 1973. [9] F. M. Callier, W. S. Chan, and C. A. Desoer. Input-output stability theory of interconnected systems using decomposition techniques. IEEE Transactions on Circuits and Systems, CAS-23(12):714–729, December 1976. [10] F. M. Callier, W. S. Chan, and C. A. Desoer. Input-output stability of inter- connected systems using decompositions: An improved formulation. IEEE Transactions on Automatic Control, AC-23(2):150–163, April 1978. [11] Claudio A. Canizares. On bifurcation, voltage collapse and load modeling. IEEE Transactions on Power Systems, 10(1):512–522, 1995. 70 [12] B. Carreras, D. Newman, I. Dobson, and A. B. Poole. Evidence for self- organized criticality in a time series of electric power system blackouts. IEEE Transactions on Circuits and Systems, 51(9), September 2004. [13] M. Chakraborty, D. Maiti, A. Konar, and R. Janarthanan. A study of the Grunwald-Letnikov definition for minimizing the effects of random noise on fractional order differential equations. In Information and Automation for Sustainability, 2008. ICIAFS 2008. 4th International Conference on, pages 449–456, Dec 2008. [14] J. Chen, J. Thorp, and M. Parashar. Analysis of electric power system dis- turbance data. In Proceedings of the 34th Hawaii International Conference on System Sciences, 2001. [15] Z. Chen. Issues of connecting wind farms into power systems. IEEE/PES Transmission and Distribution Conference and Exhibition: Asia and Pacific, 2005. [16] H.-D. Chiang, I. Dobson, R. J. Thomas, and J. S. Thorp. On voltage collapse inelectricpowersystems. IEEE transactions on Power Systems, 5(2):601–611, 1990. [17] R. Christie. Uw power system test case archive. Available: http://www.ee.washington.edu/research/pstca/. [18] F. Cipriano, H. Ouerdiane, and R. Vilela Mendes. Stochastic solution of a KPP-type nonlinear fractional differential equation. Fractional Calculus and Applied Analysis, 12(1):47–56, 2009. [19] E. Cotilla-Sanchez, P. Hines, and C. M. Danforth. Predicting critical tran- sitions from time series synchrophasor data. IEEE Transactions on Smart Grid, 3(4), December 2012. [20] I. Dobson. The irrelevance of electric power system dynamics for the loading margin to voltage collapse and its sensitivities. NOLTA, IEICE, 2(3):263–280. [21] I. Dobson. Observations on the geometry of saddle node bifurcation and voltage collapse in electrical power systems. IEEE Transactions on Circuits and Systems I: Fundamental Theory and Applications, 39(3):240 – 243, March 1992. [22] J. M. Drake and B. D. Griffen. Early warning signals of extinction in deteri- orating environments. Nature, 467, September 2010. [23] D. I. Kim et al. Wavelet-based event detection method using pmu data. IEEE Transactions on Smart Grid, October 2015. 71 [24] JA Jiang et al. An adaptive PMU based fault detection/location technique for transmission lines part i: Theory and algorithms. IEEE Transactions on Power Delivery, 15(2), April 2000. [25] V. Dakos et al. Slowing down in spatially patterned ecosystems at the brink of collapse. The American Naturalist, 177(6), June 2011. [26] P. Fairley. Sniffing out grid attacks. IEEE Spectrum, pages 13–15, Sept. 2016. [27] B. Fardanesh. Future trends in power systems control. IEEE Computer Appli- cations in Power (CAP), 15(3):24–31, July 2002. [28] F.D. Galiana. Load flow feasibility and the voltage collapse problem. IEEE Proceedings of 23rd Conference on Control and Decision, pages 485–487, December 1984. [29] F. Gorenflo and F. Mainardi. Fractional calculus: Integral and differential equations of fractional order. In A. Carpinteri and F. Mainardi, editors, Fractals and Fractional Calculus in Continuum Mechanics, pages 223–276. Springer, New York, 1997. [30] W. M. Grady. Synchrophasor event screening examples. http://web.ecs.baylor.edu/faculty/grady/. Accessed: 2015-11-6. [31] W. M. Grady and D. Costello. Implementation and application of an inde- pendent texas synchrophasor network. 63rd Annual Conference for Protective Relay Engineers, 2010. [32] K. H. Hamed. The distribution of Kendall’s tau for testing the significance of cross-correlation in persistent data. Hydrological Sciences Journal, 56(5):841– 853, 2011. doi: 10.1080/02626667.2011.586948. [33] D. J. Hill. Nonlinear dynamic load models with recovery for voltage stability studies. IEEE Transactions on Power Systems, 8(1):166–176, February 1993. [34] A.R.Ives. Measuringresilienceinstochasticsystems. Ecological Monographs,, 65(2):217–233, May 1995. [35] A . K. Kamath, J. R. Gandhi, A. S. Bohra, A. V. Goel, D. U. Patil, O. V. Kulkarni, and J. O. Chandle. Modeling of transformer characterisitics using fractional order transfer functions. In IEEE International Conference on Con- trol and Automation, pages 2123–2128, Christchurch, New Zealand, December 2009. 72 [36] R. G. Kavasseri and R. Nagarajan. Evidence of crossover phenomena in wind- speed data. IEEE Transactions on Circuits and Systems, 51(11), November 2004. [37] M. G. Kendall. A new measure of rank correlation. Biometrika, 30(1-2):81–93, 1938. [38] A. Király and I. M. Jánosi. Detrended fluctuation analysis of daily tem- perature records: Geographic dependence over australia. Meteorology and Atmospheric Physics, 88:119–128, April 2005. [39] S. G. Krantz and H. R. Parks. A Primer of Real Analytic Functions (Second Edition). Birkhäuser Advanced Texts, Boston, Basel, Berlin, 2002. [40] P. Kundur et al. Definition and classification of power system stability. IEEE Transactions on Power Systems, 19(3):1287–1401, August 2004. [41] H.G. Kwatny, A.K. Pasrija, and L.Y. Bahar. Static bifurcations in electric power networks: Loss of steady-state stability, and voltage collapse. IEEE Transactions on Circuits and Systems, CAS-33(10):981–991, October 1986. [42] F. Mainardi, Y. Luchko, and G. Pagnini. The fundamental solution of the space-time fractional diffusion equation. Fractional Calculus and Applied Analysis, 4:153–192, 2001. [43] R. Vilela Mendes and L. Vazquez. The dynamical nature of a backlash system with and without fluid friction. Nonlinear Dynamics, 47:363–366, 2007. [44] Joaquin Mur-Amada and JesÞs SallÃąn-Arasanz. Power fluctuation in a wind farm compared to single turbine. In Gesche Krause, editor, From turbine to wind farms—Technical requirements and spin-off products, chapter 6, pages 101–132. Intech, April 2011. ISBN 978-953-307-237-1; DOI: 10.5772/15957. [45] C.-K. Peng et al. Mosaic organization of dna nucleotides. Phys Rev E, 49(2), 1994. [46] C.-K. Peng, S. Havlin, H. E. Stanley, and A. L. Goldberger. Quantification of scaling exponents and crossover phenomena in nonstationary heartbeat time series. Chaos, 5, 1995. [47] V. Pipiras and M. S. Taqqu. Fractional calculus and its connections to frac- tional brownian motion. In P. Doukhan, G. Oppenheim, and M. S. Taqqu, editors, Theory and Applications of Long-Range Dependence, pages 165–201. BirkhÃďuser, Boston, 2003. 73 [48] I. Podlubny. Geometric and physical interpretation of fractional integra- tion and fractional differentiation. Fractional Calculus and Applied Analysis, 5:367–386, 2002. [49] The École polytechnique fédérale de Lausanne. PMU data from EPFL cam- pus. http://nanotera-stg2.epfl.ch/, Online; accessed March 8, 2016. [50] D. Popović, I. A. Hiskens, and D. J. Hill. Investigation of load-tap changer interaction. Electrical Power & Energy Systems, 18(2):81–97, 1996. [51] G. Samorodnitsky and M. Taqqu. Stable Non-Gaussina Random Processes– Stochastic Models with Infinite Variances. Chapman& Hall, CRC, 1994. [52] M. Schwalbe. Mathematical Science Research Challenges for the Next- Generation Electric Grid: Summary of a Workshop. The National Academies Press, Washington, DC, 2015. [53] L. Shalalfeh, P. Bogdan, and E. Jonckheere. Evidence of long-range depen- dence in power grid. In Power and Energy Society General meeting (PESGM), Boston, MA, July 2016. IEEE Power & Energy Society. [54] L. Shalalfeh, P. Bogdan, and E. Jonckheere. Kendall’s tau of frequency Hurst exponent as blackout proximity margin. In IEEE International Conference on Smart Grid Communications, Sydney, Australia, November 2016. [55] L. Shalalfeh and E. Jonckheere. The existence of a voltage collapse solution in the static-dynamic gap. In American Control Conference (ACC), pages 4126–4131, Boston, MA, July 2016. [56] G. P. Sillito. The distribution of Kendall’s τ coefficient of rank correlation in rankings containing tie. Biometrika, 34(1/2):36–40, January 1947. [57] A. Sode-Yome, N. Mithulananthan, and K. Y. Lee. A maximum loading mar- gin method for static voltage stabilityu in power systems. IEEE Transactions on Power Systems, 21(2):799–808, May 2006. [58] J. E. Tate and T. J. Overbye. Line outage detection using phasor angle mea- surements. IEEE TRANSACTIONS ON POWER SYSTEMS, 23(4), Novem- ber 2008. [59] M. S. Tavazoei and M. Haeri. A proof for non existence of periodic solutions in time invariant fractional order systems. Automatica, 46:1886–1890, 2009. [60] C.W. Taylor. Power System Voltage Stability. McGraw-Hill, 1993. 74 [61] R Development Core Team. R: A language and environment for statistical computing. https://www.r-project.org, Online; accessed May 22, 2016. [62] A.TiranuchitandR.J.Thomas. Aposturingstrategyagainstvoltageinstabil- itiesinelectricpowersystems. IEEE Transactions on Power Systems, 3(1):87– 93, 1988. [63] P. D. Valz and A. I. McLeod. A tremendously simpli- fied derivation of the variance of Kendall’s τ. Available at http://www.stats.uwo.ca/faculty/aim/vita/ps/kendall.pdf. [64] G. S. Watson. Smooth regression analysis. Sankhya A, 26(4):359–372, 1964. [65] R. Weron. Measuring long-range dependence in electricity prices. in H. Takayasu ed., "Empirical Science of Financial Fluctuations" (Springer-Verlag Tokyo, 2002), pages 110–119. [66] C. Wissel. A universal law of the characteristic return time near thresholds. Oecologia, 65:101–107, 1984. 75
Abstract (if available)
Abstract
Voltage collapse, one of the critical phenomena that threatens the power infrastructure, manifests itself by a sudden and fast collapse of the system voltage. The simplest bus system—one generator, one line, and one load—reveals a feedback structure, which, depending on the nature of the load, could go in a voltage collapse resembling an “anti van der Pol"" behavior. A static-dynamic compromise for the load model is provided by an “impedance"" describing function that has the advantage of being derived from experimental data collected in a real grid environment. This reveals, among other things, a noninteger exponent of the frequency requiring an extension for the model to be usable for a collapsing solution away from the imaginary axis. ❧ We argue that an aggregation effect of all of the loads in a grid endows a single load with grid-characteristics properties in addition to the usual load-specific properties. Topologically, the hidden feedback structure of a bus model reveals that the resulting digraph is strongly connected, meaning that all loads are intertwined in a single system that cannot be decomposed into islands. ❧ To build more accurate dynamic load models in the power system, we examine the mathematical characteristics of the power system variables. Toward this end, we demonstrate that the voltage magnitude, frequency, and voltage phase angle of Phasor Measurement Unit data exhibits long-range dependence. Our findings call for the development of a fractal modeling approach of the smart grid. ❧ Last but not least, we introduce a novel method to assess the distance to blackout or other instability of the smart grid. Based on the existence of long-range correlation in the PMU data, we exhibit an increase in the frequency Hurst exponent—quantified by Kendall’s tau rank correlation coefficient—before the blackout manifests itself in the grid powering down. Probably as the major result of this thesis, High Kendall’s tau of the frequency Hurst exponent is here proposed as an early warning signal for blackout.
Linked assets
University of Southern California Dissertations and Theses
Conceptually similar
PDF
Probabilistic data-driven predictive models for energy applications
PDF
System stability effect of large scale of EV and renewable energy deployment
PDF
Electric vehicle integration into the distribution grid: impact, control and forecast
PDF
Understanding dynamics of cyber-physical systems: mathematical models, control algorithms and hardware incarnations
PDF
Curvature analysis of the power grid: congestion, congestion management, and line tripping cascade anticipation
PDF
Predictive real-time assessment for power grid monitoring and situational awareness using synchrophasors with partial PMU observability
PDF
Data-driven methods for increasing real-time observability in smart distribution grids
PDF
Prediction models for dynamic decision making in smart grid
PDF
Utility – customer interface strategy and user application for electric vehicles and distributed generation management
PDF
Model-driven situational awareness in large-scale, complex systems
PDF
Novel and efficient schemes for security and privacy issues in smart grids
PDF
Nonlinear control of flexible rotating system with varying velocity
PDF
A function-based methodology for evaluating resilience in smart grids
PDF
Defending industrial control systems: an end-to-end approach for managing cyber-physical risk
PDF
Theoretical foundations for modeling, analysis and optimization of cyber-physical-human systems
PDF
Reinforcement learning in hybrid electric vehicles (HEVs) / electric vehicles (EVs)
PDF
High-accuracy adaptive vibrational control for uncertain systems
PDF
Improving mobility in urban environments using intelligent transportation technologies
PDF
In-situ digital power measurement technique using circuit analysis
PDF
Dynamic graph analytics for cyber systems security applications
Asset Metadata
Creator
Shalalfeh, Laith
(author)
Core Title
Data-driven aggregated fractal dynamic load modeling for early warning of voltage collapse
School
Viterbi School of Engineering
Degree
Doctor of Philosophy
Degree Program
Electrical Engineering
Publication Date
02/15/2017
Defense Date
01/11/2017
Publisher
University of Southern California
(original),
University of Southern California. Libraries
(digital)
Tag
load modeling,long-range dependence,OAI-PMH Harvest,phasor measurement units,smart grid,voltage stability
Language
English
Contributor
Electronically uploaded by the author
(provenance)
Advisor
Bogdan, Paul (
committee chair
), Jonckheere, Edmond (
committee chair
), Wang, Chunming (
committee member
)
Creator Email
alshalalfeh@hotmail.com,shalalfe@usc.edu
Permanent Link (DOI)
https://doi.org/10.25549/usctheses-c40-337563
Unique identifier
UC11258383
Identifier
etd-ShalalfehL-5058.pdf (filename),usctheses-c40-337563 (legacy record id)
Legacy Identifier
etd-ShalalfehL-5058.pdf
Dmrecord
337563
Document Type
Dissertation
Rights
Shalalfeh, Laith
Type
texts
Source
University of Southern California
(contributing entity),
University of Southern California Dissertations and Theses
(collection)
Access Conditions
The author retains rights to his/her dissertation, thesis or other graduate work according to U.S. copyright law. Electronic access is being provided by the USC Libraries in agreement with the a...
Repository Name
University of Southern California Digital Library
Repository Location
USC Digital Library, University of Southern California, University Park Campus MC 2810, 3434 South Grand Avenue, 2nd Floor, Los Angeles, California 90089-2810, USA
Tags
load modeling
long-range dependence
phasor measurement units
smart grid
voltage stability