Transcript ppt presentation - The Missouri S&T Power Laboratory

Presentation at the Missouri University of Science and
Technology
PSERC
The Hilbert Transform: Applications in
the Analysis of Power Engineering
Dynamics
G. T. Heydt
Arizona State University
October, 2009
Я
2фӘ-Я/‫ﮎ‬Ǽ4
It is obvious that …
PSERC
Ǽ5
Я*
μφ21 = 2π±σg+Q=buuY(t) >> 12.1028746*x(z)
ξΩ
5
Except
πΔΨ = Ω+Ξ/0.24r ∩∏ ≠≈Ǿ*Ǽ *ф/Ө-1 \ 1.1111+Єζ/ξΩ√ќ
when n is
q34-atan(acos(z))Λcosh(q)/Ж2=0.2957352957Юq
odd or n =
8ф21/Ө-1 = ‫ﻅ‬
3


 12i 
i 
 1279gk  r    1
k 22

 PFq
q 104769
b
cosh(a tan(sinh(t * z )))

 i
i 


 0
0
0

0
 1
  * *  
   0
0
0


    g k  r    1 0
k  22
y 0


 Fq
0 

q 104769

Ξ2-3Ω8>1
57Юq
0

 1 0.111

 1.1t / 2   
0
 * 202  1.1t
  1.1

 0    1.1

0



22.5  2
   2*

  
 22
 12k 
k qw 5

    g k  r  
k  22

 Fq
q 104769
VERY IMPORTANT FOR
Z = 2₣‫שׂ‬Я
Note
1
  16.7 qw 2 * q  y
a b
r


8 qq(    )  14v * * ô/44 - 1
cd
5vô/2 - 1
2

qr
6
3
4

7
s2
8
 
2
Outline
PSERC
• Why look to transform theory for any help in
power system dynamic analysis?
• The Hilbert transform
• Some interesting mathematics
• Modal analysis, damping and stability
• Some complications
• Summary, conclusions, recommendations,
possible venues for new work in power
engineering
3
Objectives
PSERC
• To introduce the Hilbert transform in a comprehensible
way
• To discuss applications in power engineering
• To give a capsule summary of challenges in the area
4
Il existe de nombreuses façons d'afficher une
image
PSERC
transform
What is an ‘image’
•A way to see something
•A view not easily interpreted otherwise
•Trans = across Form = manifestation
TRANSFORMATION
A mapping from one
space to another
5
Il existe de nombreuses façons d'afficher une
image
PSERC
• The concept is to make calculations easier in the
transformed domain
• And not to waste too much time in transforming and
untransforming
TRANSFORMATION
A mapping from one
space to another
6
Issues in power signal identification
PSERC
Transforms are often useful
for these applications
POWER
SYSTEM
Measurements
IDENTIFICATION
Take corrective control
action, alarms, PSS
signals
The main contenders
• Fourier analysis
• Prony analysis
• Hilbert analysis
• Various control theory approaches such as observer design
7
Why use transformations?
Fourier transform
Laplace transform
Hartley transform
To convert a differential equation to an
algebraic equation
To convert the convolution integral into
something that is more easily calculated
To convert a signal with a wide frequency
bandwidth into something that has a narrow
bandwidth in the transformed domain
Fourier transform
Laplace transform
Hartley transform
Discrete Hartley and Fourier
transforms
To get rid of unbalanced three phase quantities
To make calculations easier
And to conform with widely used notation
PSERC
Walsh transform
Symmetrical
components,
Clarke’s
components
8
David Hilbert
PSERC
1862 – 1943
born in Königsberg, East Prussia
algebraic forms
algebraic number theory
foundations of geometry
Dirichlet's principle
calculus of variations
integral equations
theoretical physics and dynamics
foundations of mathematics
the Hilbert transform
9
The Hilbert transform
PSERC
1
H [ x( t )]  X ( t ) 
* x( t )
t
Some points of interest
The transformed variable is still t
The convolution integral is best performed by taking the FT of
both sides – and use the convolution property of the FT
Recall that the FT of the 1/t term is –jsgn(ω)
This can be verified by the reciprocity theorem: if f(t) and F(jω)
are transform pairs, then f(jω) and F(t) are also transform pairs
10
The Hilbert transform
PSERC
1
H [ x( t )]  X ( t ) 
* x( t )
t
The FT of the 1/πt term is –jsgn(ω)
This can be verified by the reciprocity theorem: if f(t) and F(jω)
are transform pairs, then f(jω) and F(t) are also transform pairs
These are FT
transform
pairs
11
The Hilbert transform
PSERC
1
H [ x( t )]  X ( t ) 
* x( t )
t
Therefore, one way to obtain the HT is to MULTIPLY the FT of
1/πt (namely –sgn(ω)) with the FT of x(t).
But that is easy – just reverse the signs of all the terms of the FT
of x(t) over negative values of ω. Then take the IFT if you really
need X(t).
12
Some rather interesting Hilbert
transforms
x(t)
Aet cos(d t   )
PSERC
X(t)
Special interest in
dynamic studies of
all kinds of linear
systems
Aet sin(d t   )
13
Some rather interesting properties of the
Hilbert transform
Linearity
H(ax(t))=aX(t)
PSERC
H(x(t)+y(t))=X(t)+Y(t)
Double application
When the HT is applied twice to x(t), the result
is –x(t). This is also called anti-involution.
Inverse HT
H-1 = -H
Differentiation
H(dx/dt) = d[H(t)]/dt
Convolution
H(x*y) = X*y = x*Y
The analytic
function
XA(t) = x(t) + j H[x(t)]
14
The analytic function of the decaying
sinusoid
Aet cos(d t   )
Aet sin(d t   )
The analytic function
PSERC
HT pair
XA(t) = x(t) + j H[x(t)]
Therefore
XA(t) = Aeσtcos(ωdt+φ)+j Aeσtsin(ωdt+φ)
|XA(t)|= Aeσt
This property is useful in calculating system damping on line – and
potentially in calculating PSS signals and signals that might be used
to separate systems that will break apart in uncontrolled separation.
15
For example
A 0.270 Hz decaying
sinusoid, damping factor 0.1
PSERC
The HT of this signal
The magnitude of the analytic function – plotted
on a log scale
16
For example
PSERC
Observations
•The slope of the log of the |XA| function is the value of σ, namely the negative of the
damping factor, 0.1 in this case
•The plot is obtained numerically, and only the near end values of the plot lie off the line y
= mx+b. This is due to end effects of the DFT calculation of the HT from a finite sample.
•Since the HT is in the time domain, if the damping changes at time to, the slope of the log
plot will simply change at time to.
•Since the DFT is used, as measured data become available, the oldest datum is simply
dropped out of the DFT calculation, and the new datum is brought in – in the fashion of a
sliding window.
The magnitude of the analytic function – plotted on a log
scale
17
The phase of the analytic function
PSERC
XA(t) = Aeσtcos(ωdt+φ)+j Aeσtsin(ωdt+φ)
Arg(XA(t)) = ωdt+φ
This is ωd
18
A synthetic example
A synthetic example – corrupted
by noise (‘S5’ with SNR = 2, ‘S6’
with SNR = 5). The base signal S5
is augmented with a second mode
at 0.6 Hz, unity amplitude, time
constant 8 s in S6.
PSERC
S5
f(t) et/10sin( 2π 0.27 t )  noise
Prony analysis
Signal
Frequency Identified
(Hz)
Component
(Hz)
Attenuation factor
identified (s)
0 – 50s
0 – 10s
0 – 50s
0 – 10s
S5
0.27
0.271
0.272
9.3
7.2
S6
0.27
0.60
0.272
0.604
N/A
0.566
5.1
6.4
N/A
2.4
Hilbert analysis
Signal
Component
(Hz)
Frequency
Identified (Hz)
Attenuation factor
identified (s)
0 – 50s
0 – 10s
0 – 50s
0 – 10s
S5
0.27
N/A
N/A
N/A
N/A
S6
0.27
0.60
N/A
N/A
0.272
0.602
N/A
N/A
8
10
19
Actual signal taken in a power system
after a large disturbance
M1 A measured signal
PSERC
•Successively zoomed traces
•Prony ‘sees’ potentially spurious modes – the
number is selected by the user
•Hilbert ‘beats’ Prony in computational speed
•Hilbert can identify changes in modes as the
event unfolds
•Prony assumes stationarity in the signal
•Prony has been programmed in
commercially available packages – readily
used
•Accuracy is similar between Prony and
Hilbert
(Hz)
Damping
ratio
Prony
1
2
3
4
0.23
0.30
0.49
0.77
0.001
0.72
0.012
0.023
Comment
on
amplitude
Dominant
Minor
Minor
Negligible
Hilbert
1
0.23
-0.003
Sole
Frequency
Method Component
20
Bases of assessing the tools used for
power system signal processing
•
•
•
•
•
•
•
•
PSERC
Multiple modes and modes that are near each other
Noise in the measurements
Missing measurements
Finite sample of the time domain signal (finite time window)
Three phase issues
Speed of the identification – can it be done in real time?
Suitability for control action
Accuracy of the identification
21
Execution speed
PSERC
A second measured signal: M2
These are tie line flows
Zoomed traces
•100 identifications
•On-line capability
•Hilbert generally ‘beats’ Prony in speed
•Accuracy in synthetic signals appears to
be about the same
•Does not include preprocessing
22
Time domain windowing
Time domain windowing will impact both Prony
and Hilbert analysis. The impact on Prony can
not be corrected, but there is potential for
correction in the Hilbert domain.
1
PSERC
Windowing may be
viewed as multiplication
by a rectangular pulse
p(t). Thus the signal
measured is not x(t), but
p(t)x(t)
0
23
Bedrosian’s theorem
1
PSERC
Sample length T
0
Period of oscillation To
•The signal x(t) is known only in a finite time window [0,T ]
•The Hilbert transform is x(t)p(t) where p(t) is a rectangular pulse that captures
that window
•The Hilbert transform is H[x(t)p(t)] ≈ p(t)H[x(t)] = p(t)X(t) for pulse widths that
are significant relative to the period of oscillation of x(t) , To << T
•This approximation is Bedrosian’s theorem and it is a consequence of a narrow
band model
•Under the narrow band model, X(t) changes from cosine forms to sine forms, and
the angle of the analytic function of x(t) is calculated accurately from the
arg(XA(t))
24
Bedrosian’s theorem
1
It would be nice to
reduce T as much
as possible. This
can be done via
several routes
Sample length T
0
Period of oscillation To
Preprocess
data
Capture data
?
Remove the
assumptions of
Bedrosian's
theorem
PSERC
Reduce the BW
of the signal
Modulate the signal
with a sweeping
frequency
Noise filters
Separate even and odd parts of
x(t)
Splines to envelope
Work with moving time window
the signal
and process only changes in
Remove high frequency
X(t)
signals and process
Combine with wavelet analysis
separately
25
Hilbert Huang method
PSERC
•Effectively reduces the BW of x(t) and allows high speed processing of individual
component bands of frequencies
•Programmed in a commercial prototype, and proven in a range of applications
•Although not based on Bedrosian’s theorem, the HH method breaks the signal
x(t) into component band limited signals, and processes those separately. The
HHT method uses splines and time domain ‘sifting’. These are similar to
demodulation. The preprocessing is in the time domain.
Preprocess
data
Capture data
Use peaks to demodulate the signal
Modulate the signal
with a sweeping
frequency
Noise filters
Splines to envelope
the signal
Remove high frequency
signals and process
separately
26
Hilbert Huang method
PSERC
SPLINES
•The basic idea is to develop a series of splines that span time intervals 1, 2, …,
k, … such that the signal is stationary within the spline horizon.
•Then subtract a projected modal function within each spline horizon, m1(t) =
x(t) – h1(t), m11(t) = h1(t)-h11(t), …
•Stop subtracting estimated modal functions when the Cauchy convergence test
is satisfied, and repeat over all splines
C is sufficiently small as
set by the user. This is
effectively a nonlinear low
pass filter
27
The ‘challenges’
PSERC
Fully exploit Bedrosian’s theorem
Bedrosian’s theorem seemingly allows the use of shortened time windows of data
if the product p(t)x(t) accounts for the rectangular pulse p(t) . There have been
published ways to handle products such as this – but no one has fully exploited the
results. It is possible that much shorter clips of data would be useable in obtaining
intrinsic power system modes. A side benefit: if Bedrosian’s theorem is applied to
band limited signals, the convolution property results – but it is in the time
domain.
Combine the Prony and Hilbert methods
The Prony method has been programmed, commercialized and widely used for
many years – and there are many proponents of the method. The Hilbert method
may be viewed as a ‘competitor’ by some. But there are real possibilities to use the
time specific properties of Hilbert to size the sample window for Prony, or to
obtain accurate results for nearly collocated modes, or to simply obtain a second
estimate which may be a sanity check.
28
The ‘challenges’
PSERC
Apply Titchmarsh's theorem
Titchmarsh’s theorem: if f(t) is square integrable over the real axis, then any one of the
following implies the other two:
1. The FT, F, is 0 for negative time
2. In the FT, replacing ω by x+jy, results in a function that is analytic in the complex plane
and its integral is bounded.
3. The real and imaginary parts of F(x+jy) are the HTs of each other
This theorem may allow one to calculate the HT very rapidly by construction of the analytic
function of f(t). Also, there are some consequences of autofiltering of f(t) working in the
Hilbert domain.
Solve the Riemann – Hilbert problem for this application
Form an analytic function from the even and odd parts of a signal: fe(t) and fo(t) namely
M(t)=fe+jfo. Then consider two additional functions a(t) and b(t) such that afe-bfo = c. The
question is to find a and b such that the even part of M(z) [where z replaces t and z is a
complex number] is the HT of c(t). This may allow the selection of functions a and b that
are band limited and this will allow rapid calculation of the HT of f. And this may allow
extraction of the component modes of f(t).
29
The ‘challenges’
PSERC
Perform an error analysis for the HT and HT to quantify the
accuracy of the methods
All practical uses of the HT actually use the
discrete HT. The DHT is obtained from the
DFT. For an n-point implementation, there
is a known error introduced in the DHT
calculation. This implies that some kind of
error correction may be possible. The DFT
calculation error for one type of signal is
shown in red.
Apply the method for large scale, high profile applications
The HT method has been applied in ‘laboratory’ controlled circumstances. The need to is
to apply the idea in large systems with many intrinsic modes. And implement the
calculation alongside a Prony calculation. And also to make the HT calculation an option in
commercial software.
30
Some additional potential applications of the HT
Hilbert Transformers
This is a pair of digital filters that generates outputs
u(t) and v(t) given an input x(t) where u and v are in
quadrature – that is, their joint integral is zero.
x(t)
All pass H1(jf)
PSERC
u(t)
All pass H2(jf) v(t)
Application
In real time generate a voltage that is the q-axis component of a three phase signal, and
use a power electronic amplifier to generate a signal –xq(t) which is injected in series with
the supply for ‘power conditioning’.
Hilbert Phase Modulation
x(t)
HT
Electronic
waveform
generator
Phase modulation occurs inadvertently in bus
voltages, at low frequencies, due to power swings.
The HT of a phase modulated signal is of the
Bessel function
form Jn(β)ej2πnat , Jn is a Bessel function, 2πna is
look up table
the frequency of the phase modulation. This HT
can be calculated easily in real time, and it may
Power system stabilizer
be possible to inject a signal into the
transmission system to cancel interarea
Application: a power system stabilizer 31
oscillations.
Contributors to the Hilbert method of signal
analysis
Georg Friedrich
Bernhard Riemann
1826 – 1866
Germany
黃鍔
Norden E. Huang
1942 - …
Taiwan
PSERC
Edward Charles
Titchmarsh
1899 – 1963
U. K.
Edward Bedrosian
1922 - …
U. S. A.
32
More information
PSERC
•N. E. Huang, Z. She, S. R. Long, M. C. Wu, S. S. Shih, Q. Zheng, N.-C. Yen, C. C. Tung,
H. H. Liu, “The empirical mode decomposition and the Hilbert spectrum for nonlinear
and non-stationary time series analysis,” Proc. Royal Society of London, vol. 454, pp. 903995, 1998.
•S. L. Hahn, Hilbert Transforms in Signal Processing, Boston, Artech House, 1996.
•J. Hauer, D. Trudnowski, G. Rogers, B. Mittelstadt, W. Litzenberger , J. Johnson,
“Keeping an eye on power system dynamics,” IEEE Computer Applications in Power, vol.
10, No. 4, pp. 50-54, Oct. 1997.
•A. R. Messina, V. Vittal, D. Ruiz-Vega, G. Enríquez-Harper, “Interpretation and
visualization of wide-area PMU measurements using Hilbert analysis,” IEEE
Transactions on Power Systems, vol. 21, No. 4, pp. 1763-1771, Nov. 2006.
•Timothy Browne, V. Vittal, G. T. Heydt, Arturo R. Messina, “A real time application of
Hilbert transform techniques in identifying inter-area oscillations,” Chapter 4, Interarea
Oscillations in Power Systems, Springer, New York NY, 2009, pp. 101 – 125
33
PSERC
Questions?
Comments?
34