Data driven fractal modeling for blackout 

and malicious threat detection 

E. Jonckheere

(joint work with L. Shalalfeh and P. Bogdan), 

Dept. of Electrical Engineering

University of Southern California

Los Angeles, CA 90089

jonckhee@usc.edu

http://eudoxus2.usc.edu
Plan of Action

 Fractal PMU signal analysis

 Texas & EPFL (Switzerland) normal PMU data

 Why are the PMU signals fractal???

 Fractional dynamics load modeling 

 Hidden feedbacks in power grid

 Strong connectivity of power grid graph, 

aggregating all loads

 Early warning of imminent blackout

 Indian blackout PMU data

 Shift in AR(1) coefficient and Hurst exponent.
Long-Range Dependence or Memory

(in PMU data)

Long-range memory is one of the characteristics of

fractal patterns. It relates to slow decay of the

correlation as the lag between samples increase.

Autocorrelation

Lag

Short-memory series

Autocorrelation

Lag

Long-memory series

d = 0 d = 0.3

exponential decay non exponential decay
 Number of incrementation or differentiation steps (d):

 Power Spectral Density exponent (β):

 Hurst exponent (α):

It relates to the autocorrelation of time series and the rate at 

which these decrease as the lag increases. 

Long-Range Dependence or Memory



S( f ) 

1

 f

 There are several parameters that quantify the severity of 

the fractal behavior in a time series:



ARFIMA : 1 iL

i

i1

p













1 L

d

Xt 1  iL

i

i1

q













 t

, 1  AR(1)
 Steps:

1. Subtract average and integrate the data set:

Detrended Fluctuation Analysis (DFA)

Sample Number

Data

Sample Number

Integrated Data
2. Divide the data into n equal-sized boxes and find the Linear

Least Squares (LLS) line inside each box.

3. Subtract the LLS fitting from the integrated data to generate the

detrended data:

Detrended Fluctuation Analysis (DFA)

Sample Number

Integrated Data

Detrended Data

Sample Number

yint(k)-yn(k)=yd(k)
4. Find the Root Mean Square (RMS) fluctuation of the detrended

data:

4. The second and third steps

are repeated

at different box sizes:

Detrended Fluctuation Analysis (DFA)

log10(n)

log10(F(n))



F(n) 

1

n

yd  (k)

2

k1

n





  limn

log10F(n)

log10n
Texas Synchrophasor Network

 Several PMUs are installed at 120V and 69KV over

several locations:

 Baylor University (Waco),

 Harris Substation, and

 McDonald Observatory.

 The data we analyzed here are

 voltage magnitude,

 frequency, and

 phase angle.

 The sampling rate of the data is

30 samples/second.
PMU Time Series (Texas)
Hurst Exponent (Texas)

α = 1.11 

α = 1.54 

α = 0.71 

α = 0.92 

α = 1.54 

α = 0.75 



0.5   1: long range with power law

  > 1: long range but no power law
Hurst Exponent (Texas)

 Frequency and angle data are consistent across the 3 stations. 

 Voltage definitely has higher Hurst exponent at McDonald... Why???

 Proximity of wind farm?

 Is the Hurst exponent of voltage a sign of penetration of 

renewables in the larger grid?
PMU-Based Monitoring in EPFL

(Ecole Polytechnique Fédérale de Lausanne)

PMUs installed in EPFL campus perform real time 

monitoring of the EPFL pilot smart grid.

The PMUs were

installed on medium

voltage buses (12KV)

The sampling rate is

50 samples/second
PMU Time Series (EPFL)
Hurst Exponents (EPFL)

Voltage magnitude Frequency Phase angle 



 1.20



 1.55



 1.27

Amazing consistency between the frequency α in Texas (1.54) and Switzerland (1.55)
Hurst Exponent Histograms (EPFL)

Voltage magnitude Frequency Phase angle 

mean(α)=1.23 mean(α)=1.51 mean(α)=1.23
Plan of Action

 Fractal PMU signal analysis

 Texas & EPFL (Switzerland) normal PMU data

 Why are the PMU signals fractal???

 Fractional dynamics load modeling 

 Hidden feedbacks in power grid

 Strong connectivity of power grid graph, 

aggregating all loads

 Early warning of imminent blackout

 Indian blackout PMU data

 Shift in AR(1) coefficient and Hurst exponent.
Static versus Dynamic Load Models

• Static load model:

 Constant Power 

 Constant Current 

 Constant Impedance

• Dynamic load model (Hill):



PL

KpVL

pv



QL

KqVL

qv



v p v q  0



v p v q 1



v p v q  2



TP Pd



 Pd  Ps(VL )  kp (VL )V



L



Should be 

PL VL  ,
Berg Data-Driven Load Modeling Experiment 

in a real microgrid

Electric Reduction Furnace (Berg) Mechanical Plant (Berg)

Active 

power in a load (Hill) 

voltage and reactive 

power 

in a 

P load (Hill) 

Q 

V

Q measurement

P measurement

!

! 
Load Type pv pω qv qω

Filament lamp

Fluorescent lamp

Heater

Induction motor (HL)

Induction motor (FL)

Reduction furnace

Aluminum plant

Regulated aluminum plant

1.6

1.2

2.0

0.2

0.1

1.9

1.8

2.4

0

-1.0

0

1.5

2.8

-0.5

-0.3

0.4

0

3.0

0

1.6

0.6

2.1

2.2

1.6

0

-2.8

0

-0.3

1.8

0

0.6

0.7

L L L S  P  jQ





PL KPVL

pvp



QL KQVL

qvq

Berg load model involves frequency to a 

noninteger exponent
Load Type Describing Function

Filament lamp

Fluorescent lamp

Heater

Induction motor (HL)

Induction motor (FL)

Reduction furnace

Aluminum plant

Regulated aluminum plant



1

(KpVL

0.4 q

jK VL

2

)



1

(KpVL

0.8 1

  q

jK VL

2.8

 )



1

(Kp q

jK VL

2

)



1

(KpVL

1.8 1.5

  q

jK VL

0.4 0.3

 )



1

(KpVL

1.9 2.8

  q

jK VL

1.4 1.8

 )



1

(KpVL

0.1 0.5

  q

jK VL

0.1

)



1

(KpVL

0.2 0.3

  q

jK VL

0.2 0.6

 )



1

(K pVL

0.4 0.4

  q

jK VL

0.4 0.7

 )

Impedance Describing Function
Analytic Extension of Describing Function

Crude way: 

Leaves some coefficients complex, not completely in line with formal circuit theory

Better way:

Coefficients are kept real, in line with formal circuit theory;

However, positive realness does not hold unless the load is a heater



   j

           

where (.) and (.) are real valued.

extension

A B

Y A V j B V j A V s B V s L L L L L

   

        

   

q

L

p

L

L

L L V jW V

Z

Y     

1
Can we replace s by d

dt

???

Yes, but subject to correct interpretation:

 Caputo, D*(initial conditions in terms of integer derivatives)

 Riemann-Liouville, D (initial conditions in terms of fractional derivatives)

 Grunwald-Leitnikov, d (close to ARFIMA model)



   

 

  





  



Generation

Line

1

Transmission network

Line

(load voltages)

State

1

Distribution network

(*) (*)

1 (*) 1 (*) 1

1 1

0

0

























































































m

n n

G

G

n n n n V

V

B Y

V

V

A Y

V

V

a D b D

a D b D

 

 

related
Hidden Feedback in Power Systems

 

 

Deliberately simplified 

model of the generator... 

to put the load

in the spotlight.
Feedback Model of Power System



G 

YLine YLine

0  0

2











F 

ZL 0

0 s

2











y  (I  GF)

1Gu G

F



u  u1 u2 

t



y  y1 y  2 

t 

 







Simplification: 

No back-action of the load to the generator
Towards more Complicated 

Feedback Models of Power System
Time to Conceptualize

Nominal impedance, line Zi, Yii

Connecting lines 1, 2, ..., i, ..., m

Are we sure that

    i FiGi det I  FG   det I 

:

.

:

.

:

.

:

.
D(U1

)

D(U2)

D(U3) D(U4

)

 ZLoads 

Ylines

I

L

VL

Vl

I

l

-

Egenerators

Iother D(U)

Decomposition of Digraph into 

Strongly Connected Components D(Ui)

Feedback connections, if any, are lumped 

into strongly connected subsets

No large scale feedback connections 

at the large scale of the structure graph
Bus model

Circuit model

Graph model
Bus model

Circuit model

Graph model
Bus model

Circuit model

Graph model

No loss of strong connectivity! 

Effect of Single Contingency 

Single transmission line 5-6 tripping:
Effect of Single Contingency 

Three-phase fault at Load 1:

Loss of strong connectivity: two strongly connected components! 
Loss of connectivity: two connected components! 

Effect of Double Contingency 

Double transmission line 5-6, 2-3 tripping:
Effect of Double Contingency 

Two three-phase faults at Loads 1 and 4:

Loss of strong connectivity: four strongly connected components! 
Main Theorem

Theorem: Under the conditions that 

 the bus system is connected, 

 all generators have nonvanishing internal impedance, 

and the contingencies are restricted to

 single transmission line tripping,

the graph model is strongly connected.
Plan of Action

 Fractal PMU signal analysis

 Texas & EPFL (Switzerland) normal PMU data

 Why are the PMU signals fractal???

 Fractional dynamics load modeling 

 Hidden feedbacks in power grid

 Strong connectivity of power grid graph, 

aggregating all loads

 Early warning of imminent blackout

 Indian blackout PMU data

 Shift in AR(1) coefficient and Hurst exponent.
Critical Transition in

Harvested Population
blackout!

2012 Indian Blackout

 The blackout occurred on July 30, 2012 and affected

more than 300 million people living in Northern India.
Frequency (Hz) Residuals AR(1) coefficient

window

Increase in Autoregressive Coefficient 

before Blackout 
Frequency (Hz) Hurst exponent

window

Increase in Hurst Exponent 

before Blackout 
 Concordant pair xi > xj & yi > yj or xi < xj & yi < yj

 Discordant pair xi > xj & yi < yj or xi < xj & yi > yj

t =

# of concordant pairs - # of discordant pairs

n(n-1) / 2

Kendall’s tau

 Kendall’s tau is a rank correlation coefficient that is

used to measure—in a statistically meaningful sense—

the ordinal association between two datasets, {(t

i

,αi)}.

 Assuming that we have n pairs of x and y data

 ((x1,y1); (x2,y2); ...; (xn,yn)),

 Kendall’s tau is defined as
AR(1) coefficient

time (minutes)

Hurst exponent

time (minutes)

τ = 0.92 τ = 0.86

Kendall’s Tau of

AR(1) Coefficient versus Hurst Exponent
AR(1) versus Hurst Exponent

Sample Distributions

Normal frequency data Frequency data before blackout
Conclusions

 The frequency appears to be the most relevant data to

anticipate blackout.

 The fundamental observation is that the τ(AR(1)) and the τ(α) of

the frequency blackout data point are shifted to the right of the

empirical distributions of the Kendall tau of AR(1) and Hurst

exponent of normal frequency data.

 The shift is more pronounced for the Hurst data.

 The Hurst exponent of the frequency data appears the best

bet to anticipate blackout.

 $1,000,000 question: Could it anticipate malware?

There is hope to achieve this as it was shown that

during Distributed Denial of Service (DDoS) and UDP

flooding attacks the Akaike/Kolmogorov informational

statistics of the link utilization signals changed!
Future Work

 With more blackout data points, we hope to demonstrate—with

enough confidence—that the empirical distributions of the

normal and blackout Hurst frequency data are random draws

from different distributions.

Kendall tau of Hurst exponent of frequency

normal blackout
Thank you!

Questions?

jonckhee@usc.edu

Kendall tau of Hurst exponent of frequency

normal blackout
Conclusions

 The fractal behavior of the PMU signals is puzzling...

 Its potential for anticipating black-out and/or cyber attacks 

has been demonstrated.

 So, it is of paramount importance to understand why the PMU 

signals are fractal.

The Berg load models provide a clue with their fractional 

exponents of .

In the Berg experiment, the load is modeled in its 

microgrid environment.

The aggregation of the loads combines a great many 

lumped parameter circuit elements to make distributed 

parameter elements.