ABRUPT CHANGE DETECTION IN AUTOMATIC DISTURBANCE RECOGNITION IN ELECTRICAL POWER SYSTEMS
by
Abhisek Ukil
Submitted in partial fulfilment of the requirements for the degree
DOCTOR TECHNOLOGIAE
in the
Department of Mathematical Technology FACULTY OF NATURAL SCIENCES TSHWANE UNIVERSITY OF TECHNOLOGY
Supervisor: Prof R Zivanovic Co-supervisor: Prof SV Joubert
October 2005
DECLARATION BY CANDIDATE
“I hereby declare that the thesis submitted for the degree DTech, at Tshwane University of Technology, is my own original work and has not previously been submitted to any other institution of higher education. I further declare that all sources cited or quoted are indicated and acknowledged by means of a comprehensive list of references.”
Abhisek Ukil
Copyright© Tshwane University of Technology 2005
ii
I would like to dedicate this work to
my parents, Mrs Ranjana Ukil and Mr Ajit Kumar Ukil, my constant source of inspiration, my aunty, Mrs Sadhana Dutta, for her inspiration, and my wife, Sangita, for her selfless support and encouragement.
iii
ACKNOWLEDGEMENTS My sincere thanks to various people who played a role in this interesting study and project. Additionally, I would like to convey my earnest gratitude and appreciation to the following people and organisations:
9 Prof Zivanovic, my supervisor, and Prof Joubert, my co-supervisor, for being such supportive mentors, both technically and personally
9 Tshwane University of Technology and the National Research Foundation (NRF) for providing the financial support to make this study possible
9 Eskom Transmission, Germiston, South Africa, for their important technical support, and for providing all the real disturbance signal recordings for the study
9 My friends, Jaco, Willy, Predrag, Paul and Christian for their constant help, multidimensional support and friendship.
iv
ABSTRACT The analysis of faults and disturbances in power systems is a fundamental foundation for a secure and reliable electrical power supply. The Automatic Disturbance Recognition Project focuses on automatic disturbance recognition and analysis based on the disturbance signals obtained from the digital fault recorders in the power transmission network in South Africa.
In an automated recognition-oriented task, the first step is to apply the abrupt change detection algorithms to segment the fault recordings into different segments, such as, the prefault segment, fault, after the circuit-breaker opening, after auto-reclosure of the circuitbreakers, and so forth. The next step involves constructing the appropriate feature vectors for the different segments. Finally, pattern-matching algorithms are applied using these feature vectors to accomplish the disturbance recognition and analysis tasks. The main focus of this thesis is to take the first step towards an automatic disturbance analysis, namely abrupt change detection-based segmentation. The aim is to accurately estimate the time-instants of the changes in the signal model parameters during different events of the disturbance and accordingly segment the signal. This is critical for improving the fault recognition rate and automatic analysis quality. It also provides scope for analysing certain kinds of disturbances directly from the segmented recordings. This work proposes and establishes various new and customised techniques of detecting the abrupt changes in the electrical power systems disturbance signals. Using those techniques, event-specific automatic signal segmentation and various applications like synchronisation, relay performance monitoring and such like are performed. Commercial implementation of the project as Application Service Provider (ASP) solution is also proposed.
v
CONTENTS
PAGE ACKNOWLEDGEMENTS ..........................................................................................................................IV ABSTRACT .................................................................................................................................................... V LIST OF FIGURES.......................................................................................................................................XI LIST OF TABLES....................................................................................................................................... XV GLOSSARY................................................................................................................................................XVI CHAPTER 1
INTRODUCTION.............................................................................................................. 1
1.1
OBJECTIVES OF THE THESIS ........................................................................................................ 2
1.2
MAJOR CONTRIBUTIONS............................................................................................................. 4
1.3
ORGANISATION OF THE THESIS................................................................................................... 6
CHAPTER 2
PROJECT BACKGROUND.............................................................................................. 9
2.1
BACKGROUND AND MOTIVATION............................................................................................... 9
2.2
DISTURBANCE ANALYSIS SCHEME AT ESKOM, SOUTH AFRICA ............................................. 10
2.3
EXISTING DISTURBANCE ANALYSIS SYSTEM ........................................................................... 13
2.4
PROPOSED AUTOMATED ANALYSIS SYSTEM ............................................................................ 15
2.5
RESEARCH TOWARDS AUTOMATED ANALYSIS SYSTEM ........................................................... 17
2.6
OVERVIEW OF THE AUTOMATED ANALYSIS SYSTEM ............................................................... 18
2.6.1
Disturbance Recordings from Digital Fault Recorders........................................................... 19
2.6.2
Abrupt Change Detection-based Segmentation ....................................................................... 22
2.6.3
Feature Vector Construction ................................................................................................... 23
2.6.4
Pattern Matching and Disturbance Recognition ..................................................................... 25
2.7
SUMMARY ................................................................................................................................ 27
CHAPTER 3 3.1
ABRUPT CHANGE DETECTION................................................................................. 28 INTRODUCTION ......................................................................................................................... 28
vi
3.2
VARIOUS TECHNOLOGIES ......................................................................................................... 28
3.3
PROBLEM STATEMENT ............................................................................................................. 30
3.4
ABRUPT CHANGE DETECTION: SIMPLE METHODS .................................................................... 31
3.4.1
Superimposed Current Quantities ........................................................................................... 31
3.4.2
The Linear Prediction Error Filter.......................................................................................... 33
3.4.3
The Adaptive Whitening Filter................................................................................................. 36
3.5
ABRUPT CHANGE DETECTION: THE LINEAR MODEL BASED APPROACH ................................... 40
3.5.1
Additive Spectral Changes....................................................................................................... 40
3.5.2
Autoregressive (AR) Modelling and Joint Segmentation......................................................... 42
3.5.3
State-Space Modelling............................................................................................................. 45
3.6
ABRUPT CHANGE DETECTION: THE MODEL-FREE APPROACH .................................................. 47
3.6.1
Support Vector Machines ........................................................................................................ 47
3.7
ABRUPT CHANGE DETECTION: NONPARAMETRIC APPROACH................................................... 50
3.7.1
The Fourier Transform............................................................................................................ 50
3.7.2
The Wavelet Transform ........................................................................................................... 51
3.8
SUMMARY ................................................................................................................................ 54
CHAPTER 4
THE RECURSIVE IDENTIFICATION METHOD ..................................................... 55
4.1
INTRODUCTION ......................................................................................................................... 55
4.2
RECURSIVE IDENTIFICATION .................................................................................................... 57
4.3
IMPLEMENTATION STEPS .......................................................................................................... 62
4.4
APPLICATION RESULTS ............................................................................................................ 63
4.5
COMMENTS ON APPLICATION RESULTS .................................................................................... 66
4.6
SUMMARY ................................................................................................................................ 67
CHAPTER 5
THE WAVELET TRANSFORM METHOD................................................................. 69
5.1
INTRODUCTION ......................................................................................................................... 69
5.2
WAVELET TRANSFORM ANALYSIS ........................................................................................... 70
5.2.1
The Continuous Wavelet Transform ........................................................................................ 70
5.2.2
The Discrete Wavelet Transform ............................................................................................. 71
5.2.3
Multiresolution Signal Decomposition and Quadrature Mirror Filter ................................... 72
vii
5.3
SIGNAL DECOMPOSITION.......................................................................................................... 73
5.4
APPLICATION OF THE THRESHOLD METHOD............................................................................. 76
5.5
THE APPLICATION RESULTS ..................................................................................................... 77
5.6
COMMENTS ON APPLICATION RESULTS .................................................................................... 81
5.7
SUMMARY ................................................................................................................................ 83
CHAPTER 6
ADAPTIVE ABRUPT CHANGE DETECTION ........................................................... 84
6.1
THE ADAPTIVE WHITENING FILTER ......................................................................................... 84
6.2
APPLICATION OF THE ADAPTIVE WHITENING FILTER ............................................................... 86
6.3
THE ADJUSTED HAAR WAVELET .............................................................................................. 89
6.3.1
Overview of the Haar Wavelet................................................................................................. 89
6.3.2
Adjustment to the Haar Wavelet .............................................................................................. 90
6.3.3
Compact Support ..................................................................................................................... 91
6.3.4
Orthogonality .......................................................................................................................... 91
6.3.5
Perfect Reconstruction ............................................................................................................ 92
6.3.6
The Adjusted Scaling Function................................................................................................ 93
6.3.7
The Adjusted Wavelet Function ............................................................................................... 96
6.4
APPLICATION OF THE ADJUSTED HAAR WAVELET ................................................................. 100
6.4.1
Comments on the Application Results ................................................................................... 103
6.5
SUMMARY .............................................................................................................................. 104
CHAPTER 7
THE COMPLETE ALGORITHM................................................................................ 106
7.1
INTRODUCTION ....................................................................................................................... 106
7.2
THE DISTURBANCE SIGNAL READ MODULE ........................................................................... 107
7.3
THE SIGNAL REPRESENTATION ALGORITHMS ........................................................................ 108
7.3.1
The Recursive Identification method ..................................................................................... 108
7.3.2
Wavelet Transform method.................................................................................................... 108
7.3.3
The Adaptive Whitening Filter method.................................................................................. 109
7.3.4
The Adjusted Haar Wavelet method ...................................................................................... 109
7.4
THE THRESHOLD CHECKING ALGORITHM .............................................................................. 110
7.5
HEURISTIC SMOOTHING FILTERING ........................................................................................ 110
viii
7.6
THE DECISION-MAKING ALGORITHM..................................................................................... 112
7.7
SUMMARY .............................................................................................................................. 118
CHAPTER 8
APPLICATIONS ............................................................................................................ 119
8.1
ABRUPT CHANGE DETECTION-BASED SEGMENTATION .......................................................... 119
8.2
SYNCHRONISATION ................................................................................................................ 121
8.3
RELAY PERFORMANCE MONITORING ..................................................................................... 126
8.3.1
Fastest Relay Operating Time ............................................................................................... 126
8.3.2
Auto-reclosing of the Circuit-Breakers.................................................................................. 127
8.3.3
Main-1 and Main-2 Relay Operation .................................................................................... 128
8.4
DISTURBANCE ANALYSIS EXAMPLES ..................................................................................... 129
8.4.1
Cleared Single Phase Fault ................................................................................................... 129
8.4.2
Uncleared Single Phase Fault ............................................................................................... 130
8.4.3
Circuit-Breaker Restrike........................................................................................................ 131
8.4.4
Reactor Ring Down ............................................................................................................... 132
8.4.5
Capacitive Voltage Transformer Transient Behaviour.......................................................... 134
8.4.6
Energising of a Transformer ................................................................................................. 135
8.5
ANALYSIS OF POWER SIGNALS FROM THE MEXICAN POWER NETWORK ................................ 136
8.6
SUMMARY .............................................................................................................................. 139
CHAPTER 9
COMMERCIAL IMPLEMENTATION ...................................................................... 140
9.1
INTRODUCTION ....................................................................................................................... 140
9.2
OUTCOMES OF COMMERCIAL IMPLEMENTATION .................................................................... 141
9.2.1
First Objective Outcome: Automatic Fault Analysis System ................................................. 141
9.2.2
Second Objective Outcome: Application Service Provider ................................................... 142
9.3
THE APPLICATION SERVICE PROVIDER (ASP)........................................................................ 143
9.4
OVERVIEW OF ASP................................................................................................................. 143
9.4.1
ASP Categories...................................................................................................................... 143
9.4.2
Benefits .................................................................................................................................. 144
9.4.3
Key Features.......................................................................................................................... 145
9.4.4
Key Technologies................................................................................................................... 147
ix
9.5
WEB SERVICES ....................................................................................................................... 148
9.5.1
Introduction ........................................................................................................................... 148
9.5.2
Web Service Technology........................................................................................................ 149
9.5.3
Web Service Architecture ...................................................................................................... 149
9.5.4
Web Service Building Blocks ................................................................................................. 151
9.5.5
Web Service Development ..................................................................................................... 152
9.6
THIN CLIENT COMPUTING ...................................................................................................... 153
9.6.1
Tarantella® ............................................................................................................................ 153
9.6.2
Windows® Terminal Services................................................................................................. 155
9.6.3
Citrix® .................................................................................................................................... 155
9.7
WINCONNECT® SERVER XP™................................................................................................. 156
9.8
APPLICATION IN COMMERCIAL IMPLEMENTATION ................................................................. 158
9.9
HIGHLIGHTS AND ADVANTAGES OF COMMERCIAL IMPLEMENTATION ................................... 160
9.10
SUMMARY .............................................................................................................................. 161
CHAPTER 10
CONCLUSION........................................................................................................... 162
10.1
AUTOMATIC DISTURBANCE RECOGNITION IN ELECTRICAL POWER SYSTEMS ........................ 162
10.2
ABRUPT CHANGE DETECTION BASED SIGNAL SEGMENTATION .............................................. 164
10.2.1
Abrupt Change Detection Algorithms............................................................................... 164
10.2.2
Applications of Abrupt Change Detection Algorithms...................................................... 165
10.3
FUTURE WORK ....................................................................................................................... 166
10.3.1
On-line Abrupt Change Detection in Electrical Power Systems....................................... 167
10.3.2
Early Disturbance Prediction and Prevention.................................................................. 167
10.3.3
Application in other Domains........................................................................................... 167
BIBLIOGRAPHY ....................................................................................................................................... 168 APPENDIX A: COMTRADE ................................................................................................................... 180 APPENDIX B: MATLAB® IMPLEMENTATION CODES .................................................................. 183
x
LIST OF FIGURES PAGE FIGURE 2.1: Disturbance Analysis Scheme at National Control, ESKOM, showing the existing systems in solid-line blocks and proposed systems in dotted-line blocks. ................................................. 12 FIGURE 2.2: Architecture of the proposed Automated Disturbance Analysis System................................... 18 FIGURE 2.3: Digital Fault Recorder data file format. .................................................................................... 21 FIGURE 2.4: Feature Vector construction using the Semi-Parametric Approach. ......................................... 25 FIGURE 2.5: Algorithm Flowchart of the SVM-based Fault Classifier. ........................................................ 27 FIGURE 3.1: The Superimposed Current Quantities technique which calculates the difference of two sample points, exactly one cycle (50 Hz system) apart......................................................................... 32 FIGURE 3.2: Fault Detection Logic using Superimposed Current Quantities................................................ 33 FIGURE 3.3: LPE-based Transient Monitoring of a Current Signal............................................................... 36 FIGURE 3.4: Adaptive Whitening Filter-based Transient Monitoring. .......................................................... 39 FIGURE 3.5: Method and Model of generating a Quasi-stationary Signal..................................................... 43 FIGURE 3.6: Practical Application Result for Joint Segmentation using the EM Algorithm......................... 45 FIGURE 3.7: A Geometrical depiction of the SVM-based Abrupt Change Detection Algorithm.................. 49 FIGURE 3.8: An Application Result (Music Signals) for SVM-based Abrupt Change Detection. ................ 50 FIGURE 3.9: Basis Functions for Fourier Transform (sine wave) and Wavelet Transform (db10). .............. 52 FIGURE 3.10: Time, Frequency, STFT and Wavelet views of Signal Analysis............................................. 53 FIGURE 4.1: Segmentation of RED-Phase Voltage Signal using the Recursive Identification Method. ....... 65 FIGURE 4.2: Segmentation of RED-Phase Current Signal using the Recursive Identification Method......... 66 FIGURE 5.1: Multiresolution Signal Decomposition realised by Quadrature Mirror Filter banks................. 73 FIGURE 5.2: Daubechies 1 (Haar) wavelet: (a) Scaling function and (b) Wavelet function.......................... 74 FIGURE 5.3: Daubechies 4 wavelet: (a) Scaling function and (b) Wavelet function. .................................... 76 FIGURE 5.4: Segmentation of the RED-Phase Current Signal using the Wavelet Transform Method. ......... 79 FIGURE 5.5: Segmentation of the RED-Phase Voltage Signal using the Wavelet Transform Method. ........ 80 FIGURE 6.1: Architecture of Abrupt Change Detection based on the Adaptive Whitening Filter................. 85 FIGURE 6.2: Frequency Response Plots of the Adaptive Whitening Filter. .................................................. 86
xi
FIGURE 6.3: Segmentation of RED-Phase Voltage Signal using the Adaptive Whitening Filter and the Wavelet Transform Method. ..................................................................................................... 87 FIGURE 6.4: Pole-Zero plot of the Adjusted Haar Wavelet Scaling Filter, for n=0, which corresponds to the Original Haar Wavelet. ............................................................................................................. 95 FIGURE 6.5: Pole-Zero plot of the Adjusted Haar Wavelet Scaling Filter with Adjusting Zeroes, for n=1, which introduces one pair of Complex Conjugate Zeroes. ....................................................... 96 FIGURE 6.6: Pole-Zero plot of the Adjusted Haar Wavelet Scaling Filter with Adjusting Zeroes, for n=2, which introduces two pairs of Complex Conjugate Zeroes. ..................................................... 96 FIGURE 6.7: The Fourier Spectrum of the Adjusted Haar Wavelet, for n=0, which corresponds to the Original Haar Wavelet. ............................................................................................................. 99 FIGURE 6.8: The Fourier Spectrum of the Adjusted Haar Wavelet, for n=1, which decreases the strong ripples of the Original Haar Wavelet. ....................................................................................... 99 FIGURE 6.9: The Fourier Spectrum of the Adjusted Haar Wavelet, for n=2, which further decreases the strong ripples of the Original Haar Wavelet. .......................................................................... 100 FIGURE 6.10: Segmentation of the Voltage Waveform during a Phase-to-Ground Fault. Plot (a) shows the Result using the Haar Wavelet, while Plot (b) shows the Result using the Adjusted Haar Wavelet. .................................................................................................................................. 101 FIGURE 6.11: Segmentation of the Voltage Waveform during a Phase-to-Ground Fault. Plot (a) shows the Result using the Haar Wavelet, while Plot (b) shows the Result using the Adjusted Haar Wavelet. .................................................................................................................................. 102 FIGURE 6.12: Segmentation of the Voltage Waveform during a Phase-to-Phase Fault. Plot (a) shows the Result using the Haar Wavelet, while Plot (b) shows the Result using the Adjusted Haar Wavelet. .................................................................................................................................. 102 FIGURE 6.13: Segmentation of the Voltage Waveform during a Phase-to-Phase Fault. Plot (a) shows the Result using the Haar Wavelet, while Plot (b) shows the Result using the Adjusted Haar Wavelet. .................................................................................................................................. 103 FIGURE 7.1: Top-down Architecture of the Complete Segmentation Algorithm. ....................................... 106 FIGURE 7.2: The Decision-Making Algorithm Flowchart........................................................................... 114 FIGURE 8.1: Abrupt Change Detection-based Segmentation of the BLUE-Phase Voltage Recording during a Phase-to-Ground Fault. ........................................................................................................... 120
xii
FIGURE 8.2: Three Segmented but Unsynchronised DFR Voltage Recordings for the same Phase-to-Ground Fault. ....................................................................................................................................... 122 FIGURE 8.3: Three Segmented and Synchronised DFR Voltage Recordings for the same Phase-to-Ground Fault. ....................................................................................................................................... 124 FIGURE 8.4: Synchronised Analogue Voltage Recording Plot and Binary Plots for the Fault Duration and Auto-Reclosing of the Circuit-Breaker. .................................................................................. 125 FIGURE 8.5: Segmented Analogue Voltage Recording during a Phase-to-Ground Fault Successfully Cleared.................................................................................................................................... 129 FIGURE 8.6: Segmented Analogue Current Recording during a Phase-to-Ground Fault not Successfully Cleared.................................................................................................................................... 130 FIGURE 8.7: Voltage (Upper Plot) and Current (Lower Plot) Recordings during a Circuit-Breaker Restrike. ................................................................................................................................................ 132 FIGURE 8.8: Schematic Oscillating Circuit Diagram of the Reactor Ring Down phenomenon. ................. 133 FIGURE 8.9: Voltage Recording during a Reactor Ring Down. .................................................................. 133 FIGURE 8.10: Segmented Voltage Recording for the Transient Behaviour of the Capacitive Voltage Transformer (CVT)................................................................................................................. 134 FIGURE 8.11: Segmented Current Recording reflecting the Energising of a Transformer, which is associated with High Magnetising Inrush Currents.................................................................................. 135 FIGURE 8.12: Active Power Flow Oscillations in the Transmission Line of the Mexican Interconnected System (MZD-DGD). ............................................................................................................. 136 FIGURE 8.13: Selected Time windows for Linear Spectral Analysis. ......................................................... 137 FIGURE 8.14: Segmentation of the Power Oscillation Signal from the Mexican Power Network, using the Adjusted Haar Wavelet Method.............................................................................................. 137 FIGURE 8.15: Segmentation of the Power Oscillation Signal from the Mexican Power Network, using the Wavelet Method (db4 mother wavelet)................................................................................... 138 FIGURE 8.16: Segmentation of the Power Oscillation Signal from the Mexican Power Network, using the Wavelet Method (db4 mother wavelet)................................................................................... 138 FIGURE 9.1: Web-Service Architecture: XML-based Interface to Back-end Services................................ 150 FIGURE 9.2: Simplified Service-Client Request-Response Mechanism. ..................................................... 151 FIGURE 9.3: Tarantella® Working Principle................................................................................................ 155
xiii
FIGURE 9.4: WinConnect® Server XP™ Working Principle. ....................................................................... 158 FIGURE 9.5: Commercial Implementation: ASP of Automatic Disturbance Analysis in Electrical Power Systems. .................................................................................................................................. 159
xiv
LIST OF TABLES PAGE TABLE 8.1: Fault Inception Time of the DFR Recordings shown in Figure 8.2.......................................... 123 TABLE 8.2: Change Time-instants of the Synchronised, Segmented Analogue DFR Recordings shown in Figure 8.3 ................................................................................................................................ 126 TABLE 8.3: Fault Duration in the Synchronised, Segmented Analogue DFR Recordings shown in Figure 8.3 ................................................................................................................................................ 126 TABLE 8.4: Circuit-Breaker Auto-Reclosing time of the Synchronised, Segmented Analogue DFR Recordings shown in Figure 8.3 ............................................................................................. 128
xv
GLOSSARY AAS
automated analysis system
ADAMO
adaptive whitening filter
ADAU
analogue data acquisition unit
ADSL
asymmetric digital subscriber line
AFMM
adaptive forgetting through multiple models
API
application programming interface
AR
autoregressive
ARMA
autoregressive moving average
ASP
application service provider
BDAU
binary data acquisition unit
CPU
central processing unit
CUSUM
cumulative sum (algorithm)
CVT
capacitive voltage transformer
CWT
continuous wavelet transform
DAU
data acquisition unit
DFR
digital fault recorder
DFT
discrete Fourier transform
DOS
disk operating system
DSL
digital subscriber line
xvi
DWT
discrete wavelet transform
ebXML
electronic business extensible markup language
EKF
extended Kalman filter
EM
expectation-maximisation
FDI
failure detection and isolation
FFT
fast Fourier transform
FIR
finite impulse response
FTP
file transfer protocol
GLR
generalised likelihood ratio
GUI
graphical user interface
HPF
highpass filter
HTML
hypertext markup language
HTTP
hypertext transfer protocol
IDE
integrated development environment
IDL
interface description language
iff
if and only if
ISDN
integrated services digital network
ISP
Internet service provider
ISV
independent software vendor
IT
information technology
LAN
local area network
xvii
LMS
least mean square
LPE
linear prediction error
LPF
lowpass filter
MAP
maximum a posteriori
MoM
mobile multicast protocol
MOV
metal oxide varistor
MSD
multiresolution signal decomposition
MZD-DGD
mexican interconnected system
NC
national control
OASIS
Organisation for the Advancement of Structured Information Standards consortium
PC
personal computer
PDA
personal digital assistant
PDF
probability density function
QMF
quadrature mirror filter
RAM
random access memory
RDC
remote desktop connection
RDP
remote desktop protocol
RPC
remote procedure call
SCADA
supervisory control and data acquisition
SLA
service-level agreement
xviii
SLT
statistical learning theory
SMS
short messaging service
SOAP
simple object access protocol
STFT
short-time Fourier transform
SV
support vector
SVC
static var compensator
SVM
support vector machine
TCP-IP
transfer control protocol –Internet protocol
TIPPS
transmission integrated plant performance system
UDDI
universal description, discovery and integration
viz.
namely
VPN
virtual private network
W3C
World Wide Web consortium
WAN
wide area network
WML
wireless markup language
WSDL
Web service description language
WT
the wavelet transform
WTS
windows terminal services
XML
extensible markup language
xix
xx
CHAPTER 1: INTRODUCTION
1
CHAPTER 1
INTRODUCTION
The analysis of faults and disturbances has always been a fundamental foundation for a secure and reliable electrical power supply. The introduction of digital recording technology opened up a new dimension in the quantity and quality of fault and disturbance data acquisition, resulting in the availability of a huge amount of new information to power systems engineers. Information from the analysis of digital records can provide muchneeded insight into the behaviour of the power system as well as the performance of protection equipment. Manual analysis of these records, however, is both time-consuming and complex. Today the challenge is to automatically convert data into knowledge, which frees the human resources to implement corrective or preventive action. The purpose of the Automatic Disturbance Recognition Project is to augment the existing semi-automatic fault analysis system at the National Control, ESKOM, South Africa, with more robust and accurate algorithms and techniques to make it fully automated. In the direction towards a recognition-oriented task, the abrupt change detection algorithms are first applied to segment the fault recordings from the digital fault recorders (DFRs) into different event-specific segments such as, the pre-fault segment, after initiation of the fault, after circuit-breaker opening, and after auto-reclosure of the circuit-breakers. Appropriate feature vectors are then constructed for the different segments, and finally, the patternmatching algorithms are applied using these feature vectors to accomplish the fault recognition and disturbance analysis tasks (Ukil & Zivanovic, 2005a).
2
CHAPTER 1: INTRODUCTION
The main focus of this thesis is taking the first step towards the proposed automatic fault analysis, namely abrupt change detection-based segmentation. The aim is to estimate exactly the time-instants of the changes in the signal model parameters during the pre-fault condition, after initiation of fault, after circuit-breaker opening and auto-reclosure. This is critical to improving the fault recognition rate and automatic analysis quality. It also provides huge scope for analysing certain kinds of disturbances directly from the segmented recordings and monitoring the performances of the protective relays before conforming with any further significant and complex feature vector analysis (Ukil & Zivanovic, 2005d). In this work, we propose and establish various new and customised techniques of detecting the abrupt changes in the electrical power systems disturbance signals and performing event-specific automatic signal segmentation. The aim of this chapter is to introduce the methods developed in this thesis. For this purpose we present the problems associated with automatic disturbance recognition in electrical power systems. We introduce the adopted approach and present the organisation of this thesis.
1.1
Objectives of the Thesis
The detection of abrupt changes in signal characteristics is a widely studied subject with many different approaches. It has a significant role to play in failure detection and isolation (FDI) systems. The aim of this thesis is to introduce abrupt change detection technology as an effective segmentation tool for use in the power systems disturbance analysis. A possible approach to recognition-oriented signal processing entails using an automatic segmentation of the signal based on abrupt change detection as the first processing step. A
CHAPTER 1: INTRODUCTION
3
segmentation algorithm splits the signal into homogeneous segments, the lengths of which are adapted to the local characteristics of the analysed signal which can be achieved either on-line or off-line. In power systems applications, traditionally this step has not received much attention. Only extremely simple technologies such as peak value detection and superimposed current quantities (Stokes-Waller, 1998) are used to identify certain sections. Explicit event-specific segmentation of the power systems disturbance signal is critical in the disturbance analysis. After effective segmentation it is possible focus on further signal processing and analysis for the specific segments, not the whole signal, improving both the speed of operation and accuracy. For the Automatic Disturbance Recognition in Electrical Power Systems Project, this is even more important. Various new algorithms and approaches for abrupt change detection-based segmentation of the power systems disturbance signals, presented in this thesis, suggest an effective approach for the automatic disturbance recognition and analysis. The main goals of this thesis are as follows: (1) to study the abrupt change detection technologies and various technologies associated with it (2) to develop algorithms for effective application of the techniques in the power systems domain (3) to implement the different algorithms and test them extensively on the real disturbance signals, obtained from the digital fault recorders of the electrical power network of South Africa
4
CHAPTER 1: INTRODUCTION
(4) to compare the results of the different segmentation algorithms based on the practical tests and other factors (5) to develop an optimised overall algorithm, considering all the algorithms for best possible segmentation for all kinds of disturbance signals (6) to examine the scope of additional processing of the disturbance signals for some effective inferences directly from the segmentation step (7) to develop a framework for the commercial implementation of the overall Automatic Disturbance Recognition in Electrical Power Systems Project, including the abrupt change detection-based segmentation for effective technology transfer to the power utilities, such as ESKOM, South Africa, which will be the principal users of the outcome of the research.
1.2
Major Contributions
Several new methods or developments are introduced in this thesis. The major contributions are as follows: •
the development of a new approach for automatic disturbance recognition in electrical power systems by conducting event-specific segmentation based on abrupt change detection (Ukil & Zivanovic, 2005a)
•
the introduction and evaluation of an abrupt change detection algorithm based on the recursive identification technique (Ukil & Zivanovic, 2005f)
CHAPTER 1: INTRODUCTION
•
5
the introduction and evaluation of an abrupt change detection algorithm using the wavelet transform and threshold method (Ukil & Zivanovic, 2005b)
•
the introduction and evaluation of an improvement technique for the abrupt change detection algorithm using the adaptive whitening filter (Ukil & Zivanovic, 2005e)
•
the introduction of a new general adjustment technique to the standard mother (Haar) wavelet. This generates a new series of adjusted mother (Haar) wavelets and also improves the accuracy and effectiveness of the wavelet transform-based algorithm for the applications addressed in this thesis (Ukil & Zivanovic, 2005g)
•
the introduction and evaluation of different direct applications of the abrupt change detection based segmentation, such as, relay performance monitoring (Ukil & Zivanovic, 2005c) and the analysis of certain kinds of disturbances (Ukil & Zivanovic, 2005d)
•
the study and development of commercial implementation of the project as Application Service Provider (ASP) solution (Ukil & Zivanovic, 2005h).
The secondary contributions are as follows: •
the study and evaluation of several methods for abrupt change detection
•
a comparison of the developed algorithms of abrupt change detection-based segmentation by means of extensive testing of the implementations on real disturbance signals
•
the development of an optimised overall algorithm, considering all the individual
6
CHAPTER 1: INTRODUCTION
algorithms developed •
the development of an effective heuristic smoothing filtering operation for the fine accomplishment of abrupt change detection-based segmentation.
1.3
Organisation of the Thesis
This thesis is organised in ten chapters: •
Chapter 1: We review the problem treated in this thesis and define the goals of the work. Finally, the organisation of the thesis is presented.
•
Chapter 2: We present the background information of the Automatic Disturbance Recognition in Electrical Power Systems Project. The background and motivation of the disturbance recognition in electrical power systems is presented. Detailed discussions about the existing semi-automatic system, the problems associated with it and the reasons behind the current research to develop a fully automatic system are provided. We define the functional specification of the proposed automated analysis system, and discuss briefly about the four major modules of the automated analysis system, abrupt change detection-based segmentation being the second step.
•
Chapter 3: We present a comparative detailed discussion of the existing research works on abrupt change detection relevant to our project. It explains the theories necessary to understand the project and its necessity on the basis of existing research.
CHAPTER 1: INTRODUCTION
•
7
Chapter 4: We look into the recursive identification technique, and develop an algorithm for abrupt change detection-based segmentation on the basis of that. The application results of the recursive identification-based algorithm on the real disturbance recordings are also presented along with comments.
•
Chapter 5: We develop a new algorithm for abrupt change detection-based segmentation based on the wavelet transform and threshold method. The application results of the wavelet transform-based algorithm on the real disturbance recordings are also presented along with comments.
•
Chapter 6: We discuss about the adaptive abrupt change detection techniques for some special disturbance signals, not showing distinct abrupt changes in the signal parameters. Two techniques are proposed. The first method is an improvement technique based on the adaptive whitening filter. In the second method, a general adjusting technique to the standard Haar wavelet is proposed. This generates a series of new adjusted Haar wavelet which can be successfully used as the basis for the wavelet transform algorithm for accomplishing the adaptive abrupt change detection task.
•
Chapter 7: We develop the overall, optimised algorithm considering all the developed algorithms for abrupt change detection based segmentation. The architectural structure of the software implementation and the smoothing filter are also discussed.
•
Chapter 8: We present the various applications of abrupt change detection-based segmentation on real disturbance signals. This includes event-specific segmentation for further signal processing towards automatic disturbance recognition, synchroni-
8
CHAPTER 1: INTRODUCTION
sation, relay performance monitoring and analysis of certain kinds of disturbances directly from the segmentation steps. •
Chapter 9: We describe the idea of commercialisation of the whole project as an application service provider (ASP) solution, to be used using the Internet. Following a general discussion about the ASP, key technologies (e.g. Web services, thinclient computing etc) for the ASP implementation are discussed in details. We also discuss the utilisation of the ASP technologies for the commercial implementation of the project (direct and indirect mode).
•
Chapter 10: We present the general conclusions of this thesis and propose possible improvements and directions for future research.
CHAPTER 2: PROJECT BACKGROUND
9
CHAPTER 2 PROJECT BACKGROUND
Automatic disturbance recognition and analysis from the recordings of the digital fault recorders plays a significant role in fast fault clearance, helping to achieve a secure and reliable electrical power supply. This chapter provides background information on the Automatic Disturbance Recognition in Electrical Power Systems Project. Application of abrupt change detection in various dimensions, the central theme of this thesis, is a critical module of the overall project. Before going into further detail on abrupt change detection, it is necessary to discuss the overall project more broadly in order to provide a clear picture of the whole project. This will help the reader to understand the importance of the abrupt change detection module, which is the focus of this thesis, in the overall project.
2.1
Background and Motivation
The first disturbance recording systems in electrical power systems were based on electromechanical and analogue technology (Barth, Ludwig & Schegner, 2003). The measurements were recorded on metallic or photosensitive papers or magnetic tapes. The accuracy, the number of signals, the recording time and the number of data records per paper roll or magnetic tape were highly restricted. Data handling and interpretation were also extremely complex processes. With the introduction of the digital recording technology, nowadays, engineers in power utilities have more data than can be processed and assimilated in the time available. Users
10
CHAPTER 2: PROJECT BACKGROUND
of knowledge extracted from recorded data are engineers and technicians working in operations, maintenance and protection departments. It is the task of operating personnel to return to service as much of the electric system as practical in the shortest time possible. The following questions concern them: Where and what is the problem? Did the line reclose and stay in operation? What equipment was operating? Was everything working correctly? If so, can it be returned to service? If not, what needs to be isolated? A complete analysis should be available for use within five minutes after the conclusion of the event. Maintenance personnel are charged with repairing and returning outage equipments to service. They require information on what is damaged or operating outside the normal parameters. The time requirement to notify the maintenance personnel is in the region of two hours or less. Protection personnel are responsible for the final assessment of the correctness of any electrical system response to a given disturbance condition. Normally the protection engineer is given adequate time to collect all the data necessary for a complete evaluation of a disturbance. Key questions that need to be answered are: Did the right thing respond in the right way? Did the wrong thing respond in the wrong way? Did the right thing respond in the wrong way? The “thing” in these questions might be a relay, relaying system, circuit breaker, or switch.
2.2
Disturbance Analysis Scheme at ESKOM, South Africa
In 1998, a 24 hour support service for National Control, ESKOM (South African power utility) was introduced to supply operating personnel with information from disturbance records analysis with a view to enhancing operational decisions, immediately identifying potential risk and improving the response time to latent weaknesses or failures (Bartylak, 2002). Since 1993 ESKOM has been installing centralised per substation disturbance re-
CHAPTER 2: PROJECT BACKGROUND
11
corders on the transmission lines (currently installed at 136 substations). The DFRs are installed on the feeder bays with an additional few installed on the Static Var Compensators (SVCs) (Stokes-Waller, 1998). The recorders monitor all voltages and currents and 30 selected protection operation outputs per line. The recorders are triggered by protection operation pulses, and some of them by dv/dt (change of phase voltage). The recording duration varies from two to 12 seconds with scanning frequency of 2.5 kHz. The implementation of the X.25 communication facility was initiated in 1995, and all recorders are at present remotely accessible (Stokes-Waller & Keller, 1998). In line with the commissioning of new recorders, the amount of mathematical analysis of recorded signals was continuously growing. Storage and retrieval of disturbance data became a real challenge (presently approximately 2,000 records a year). Semi-automatic software, OSCOP (Siemens, 1999) for data acquisition and storage in performance databases is currently available. The results from the analysis of the digital fault records are also captured in a database called TIPPS (Transmission Integrated Plant Performance System) (Keller, Henze & Zivanovic, 2005). In the fully developed automatic system, mathematical analysis of all (presently known) patterns of incorrect behaviour (disturbances) should be performed automatically and only a short message should be issued for the controllers (operating personnel), summarising the required information and knowledge. Figure 2.1 (obtained from ESKOM) shows the disturbance analysis scheme employed at the National Control, ESKOM, South Africa. In Figure 2.1, the blocks in solid lines indicate the existing systems, whereas the blocks in dotted lines indicate the proposed automatic disturbance recognition and analysis systems.
12
CHAPTER 2: PROJECT BACKGROUND
FIGURE 2.1: Disturbance Analysis Scheme at National Control, ESKOM, showing the existing systems in solid-line blocks and proposed systems in dotted-line blocks.
As indicated in Figure 2.1, when an event (Tx-event) on the power transmission network occurs, National Control (NC) will be informed via a SCADA (Supervisory Control And Data Acquisition) network and the Digital Fault Recorders (DFRs) will also be triggered. The protection engineers provide 24-hour support through the analysis of these records. It is possible to download the fault record via the X.25 communication network. The protection engineer has to conduct a manual analysis to create a report for National Control. The manual analysis is cumbersome and time-consuming, typically from one hour to 10 hours or more, depending on the complexity and severity of the event. Additionally an automatic scanning PC downloads the records which need to be saved manually at a file server. Results from the analysis of the digital fault records are also captured in the TIPPS database.
CHAPTER 2: PROJECT BACKGROUND
13
Owing to large amount of historical data available trend analysis on the performance of primary and secondary equipment can be done. An automated analysis system (shown by dotted lines in Figure 2.1) would be much faster, and could be placed in two different ways. The more efficient way is to conduct the automated analysis directly at the substation and only transmit a report to NC. This is the fastest way, but additional PCs and the communication infrastructure (for example, a LAN) are needed. The other way is to run the automated analysis on a central PC. The disadvantage is that it takes longer to download the whole record than to send a report. The main advantage is the cost factor because only one PC is needed and the communication network is much simpler. The automated analysis system should also be able to write the analysis results to the TIPPS database.
2.3
Existing Disturbance Analysis System
The majority of digital fault recorders at ESKOM are from Siemens: SIMEAS R and P531 recorders (Siemens, 1999). The abbreviation ‘SIMEAS R’ stands for SIemens MEASurement Recorder (SIEMENS new generation fault recorder), while OSCILLOSTORE P5xx stands for the earlier generation of fault acquisition systems. Both generation types are referred to by the generic name “quality recorder” (recorders used to measure the quality of a power network and contribute to maintaining that quality at a high level). The DFRs trigger because of, say, power network fault conditions, protection operations, breaker operation, and the like. Each DFR recording typically consists of 32 points binary information and analogue information in the form of voltages and currents per phase as well as the neutral current (Keller, 2004). Following the IEEE standard for Common For-
14
CHAPTER 2: PROJECT BACKGROUND
mat for Transient Data Exchange (COMTRADE) as specified by IEEE STANDARD C37.111 (1991), the DFR recordings are provided as input into the existing semi-automatic fault analysis software OSCOP (Siemens, 1999) which uses discrete Fourier analysis and superimposed current quantities (Stokes-Waller, 1998). All the disturbance analysis is done manually which is quite a complex and time-consuming process (Stokes-Waller & Keller, 1998). Also, some of the information from the OSCOP, especially that on the circuit breaker binaries are often incomplete, which is why it is not that trustworthy and needs to be handled with caution. An automated fault analysis system based only on the binaries should be avoided. The semi-automatic analysis of OSCOP provides a report that includes (Henze, 2005) the following: •
the values of all analogue channels during pre-fault, fault and post-fault period
•
the fault start and fault end
•
the fault type and involved phases
•
the fault resistance
•
the fault location
•
power values.
However, the values of the automated analysis tool of OSCOP are not that accurate as determined by the ESKOM engineers (Henze, 2005). To solve this problem, the proposed independent automated analysis of these records can provide additional value because it is far more robust, trustworthy and quicker and gives the engineer more time for detailed analysis. This is the driving point of the Automatic Disturbance Recognition Project.
CHAPTER 2: PROJECT BACKGROUND
2.4
15
Proposed Automated Analysis System
The functional specifications of the proposed automated analysis system (AAS) are described in the work of Keller, Henze and Zivanovic (2005). The first requirement for the AAS is that it must be able to read a COMTRADE file (IEEE STANDARD C37.111, 1991). The AAS must be able to automatically import the COMTRADE file from a user-defined directory, do the analysis and produce the results in a format which can be viewed by any text editor. The distribution of the results should be done through any of the following media: fax, e-mail, SMS (short messaging service), print or Web (Keller, Henze & Zivanovic, 2005). The AAS must extract from the COMTRADE file the following information (Keller, Henze & Zivanovic, 2005): •
Faulted phase(s)
•
Fault type
•
Total fault duration
•
Main 1 protection operating time
•
Main 2 protection operating time
•
Fault location
•
Fault resistance
•
DC offset
•
Breaker operating time
•
Auto reclose time.
16
CHAPTER 2: PROJECT BACKGROUND
In addition, the AAS must also determine and report on the following as described by Keller, Henze and Zivanovic (2005): •
Was the fault on the specific feeder?
•
Did the main 1 relay operate?
•
Did the main 1 relay operate correctly?
•
Did the main 2 relay operate?
•
Did the main 2 relay operate correctly?
•
Was the main 1 permissive carrier signal sent, and was it correct?
•
Was the main 2 permissive carrier signal sent, and was it correct?
•
Was the main 1 permissive carrier received, and was it correct?
•
Was the main 2 permissive carrier received, and was it correct?
•
How did the breaker operate – 1 or 3 poles?
•
Did the breaker auto reclose (ARC)?
•
What was the magnitude of the fault current?
•
What was the magnitude of the neutral current?
•
What was the magnitude of the healthy phase currents during the fault?
•
What was the depression in voltage on the faulted phase/s?
•
What was the depression in voltage on the healthy phase/s?
•
What were the dominant frequencies before and during the fault?
•
Did any of the breaker poles restrike?
CHAPTER 2: PROJECT BACKGROUND
2.5
17
Research towards Automated Analysis System
As described above, the practical problem identified by ESKOM (South Africa) and many other power utilities worldwide is as follows: How is it possible to automate mathematical analysis of data recorded during disturbances in an electrical network? To solve this practical problem we formulate the following subproblems: •
Modelling: Define mathematical models (structure and range of parameters) for a selected set of disturbances.
•
Feature extraction: Identify the disturbance model parameters (features) using recorded data.
•
Classification and analysis: Based on estimated features, automatically recognise the disturbance type and perform qualitative analysis.
The aims of the research can be summarised as follows:
to develop mathematical models for a selected set of disturbances
to develop the parameter identification technique for the specified model structure
to develop automatic disturbance recognition algorithm(s)
to develop various techniques for qualitative analysis based on the disturbance model. Examples include circuit-breaker restrike, secondary arc and reclosing failure, various resonance conditions with trapped energy, transformer inrush current condition, instrument transformer transients, fault location, performance analysis of
18
CHAPTER 2: PROJECT BACKGROUND
controllers and protection systems and analysis of electromechanical oscillations, etc
to create the software package for the automatic disturbance recognition and analysis.
2.6
Overview of the Automated Analysis System
In working towards achieving an automatic recognition-oriented task, we would first apply the abrupt change detection algorithms to the fault recordings. This would segment the fault recordings into different event-specific segments, namely pre-fault segment, after circuit-breaker opening and after auto-reclosure of the circuit-breakers. The appropriate feature vectors for the different segments would then be constructed. Finally, the patternmatching algorithm would be applied using those feature vectors to accomplish the fault recognition and associated tasks. The complete architecture of the proposed automated analysis system is shown in Figure 2.2.
FIGURE 2.2: Architecture of the proposed Automated Disturbance Analysis System.
CHAPTER 2: PROJECT BACKGROUND
19
In Figure 2.2, the architecture of the automated analysis system is a sequential top-down block diagram. Each of the blocks is described in the following subsections.
2.6.1 Disturbance Recordings from Digital Fault Recorders Digital Fault Recorders (DFRs) are highly accurate recording instruments providing sampled waveform and contact data using a relatively high sampling rate, typically above 5 kHz (a sample every 0.2 milliseconds). A DFR has all recording functions, as well as some sensing, detection and calculation capabilities. The main characteristics of DFR are as follows: •
the capability of remote interrogation for data analysis and manipulation
•
low-cost mass storage.
Because of the different nature and processing techniques of the input signals – analogue and binary – two types of data acquisition units (DAU) are designed, namely: •
Analogue (ADAU)
•
Binary (BDAU).
As discussed in Section 2.3, the majority of digital fault recorders in ESKOM are from Siemens, namely SIMEAS R and OSCILLOSTORE P531 recorders (Siemens, 1999). The following section briefly discusses main important characteristics of the individual recorders. SIMEAS R comes with 16-bit resolution and a 12.8 kHz maximum scanning frequency per channel. In SIMEAS R, the functions of a fault recorder, digital recorder, power and frequency recorder, diagnosis system and message printer are combined in one instrument.
20
CHAPTER 2: PROJECT BACKGROUND
There are two types of central unit: one with 8 analogue and 16 binary channels, and one with 32 analogue and 64 binary channels (Siemens, 1999). OSCILLOSTORE P531 comes with 8- and 12-bit resolution and a 5 kHz maximum scanning frequency. It has 31 Data Acquisition Units (DAUs) per central unit, that is, 124 analogue or 992 binary data acquisition channels (Siemens, 1999). Every feeder uses three data acquisition unit (DAU): 2 x ADAU (analogue data acquisition units used for 4 x voltage and 4 x current signals), and 1 x BDAU (binary data acquisition unit used for up to 32 digital signals). One recorder is able to monitor all values of 10 feeders. Per feeder, 8 analogue channels (3 phase voltages, 1 phase-to-Phase voltage at the busbar, 3 phase currents and the neutral current) and 32 binaries (carrier receive and send, protection operation, status of breakers…) can be recorded. The DFR measures continuously, but only a predefined part of the signal will be stored when the DFR is triggered. Usually one second of pre-fault and two seconds after the recorder has been triggered are recorded, that is, the whole record is about three seconds long (Henze, 2005). The 32 binary values are either stored as a 0 or a 1 and indicates the status of a contact, for example, breaker auxiliary contact (Keller, 2004). Binaries are divided into the following four groups: •
Main 1 distance relay contacts
•
Main 2 (back up) distance relay contacts
•
Other relays
•
Circuit breaker and other necessary signals.
CHAPTER 2: PROJECT BACKGROUND
21
Analogue values indicate the magnitude of an analogue signal (voltage or current) measured at a specific point in time (Keller, 2004). The analogue information consists of •
Voltages per phase
•
Currents per phase as well as the neutral current.
As per the IEEE COMTRADE standard, each event will have three types of files associated with it (Henze, 2005). Each of the three types carries a different class of information, namely header (*.HDR), configuration (*.CFG) and data (*.DAT). A more detailed discussion of the COMTRADE structure is provided in Appendix A. Each record (row) in the data file (*.DAT) obtained from the ESKOM DFRs has 42 data items (columns): the first column is the sample number, the second column the time in microseconds from the beginning of the record; the next eight columns are the four analogue voltage and four analogue current recordings, next 32 columns are the 32 binary points representing the contact changes. Figure 2.3 shows a sample record from a data file.
FIGURE 2.3: Digital Fault Recorder data file format.
22
CHAPTER 2: PROJECT BACKGROUND
2.6.2 Abrupt Change Detection-based Segmentation Detection of abrupt changes in signal characteristics is a much-studied subject with many different approaches. It has a significant role to play in failure detection and isolation (FDI) systems and segmentation of signals in recognition-oriented signal processing (Basseville & Nikoforov, 1993). Many of these signals are quasi-stationary in nature, that is, they are composed of segments of stationary behaviour with abrupt changes in their characteristics in the transitions between different segments. It is important to find the time instants when the changes occur and to develop models for the different segments during which the system does not change (Ukil & Zivanovic, 2005f). Segmentation of the fault recordings by detecting the abrupt changes in the characteristics of the fault recordings, obtained from the DFRs of the power network in South Africa, is the first step towards automatic disturbance recognition and analysis. Abrupt change detection algorithms collectively segment the fault recordings into different event-specific segments, namely, pre-fault segment, after initiation of fault, after circuit-breaker opening and after auto-reclosure of the circuit-breakers. Ukil and Zivanovic (2005a) conducted a comparative study of the different technologies for abrupt change detection and categorised the techniques are broadly as simple methods, linear model-based approaches, model-free approaches and nonparametric approaches. The following chapters provide detail descriptions of the techniques. Utilising the segmented recording and working on the specific segments such as, pre-fault, fault, after circuit-breaker opening, etc, is critical for improving the fault recognition rate and automatic analysis quality.
CHAPTER 2: PROJECT BACKGROUND
23
Besides facilitating further complex feature vector analysis and pattern recognition leading to automated disturbance recognition and analysis, the abrupt change detection-based segmentation procedure followed by the synchronisation step can be directly used to monitor the relay performance as studied by Ukil and Zivanovic (2005c). Abrupt change detection based segmentation also provides huge scope for analysing certain kinds of disturbances directly from the segmented recordings before conforming with any further significant and complex feature vector analysis (Ukil & Zivanovic, 2005d).
2.6.3 Feature Vector Construction The main objectives of this step are as follows: •
to develop the parameter identification techniques for various models
•
to develop a parameter identification technique in the presence of nonlinearity, such as, the case with slow time-varying amplitude and frequency
•
to extract features from sources other than analogue signals.
A practical recorded signal is corrupted with noise and nuisance components (higher frequencies due to wave transients and measurement system bias). A link between the system model and the signal model could be established through the state space theory. To filter noise and bias, an extended version of the structured total least squares method (Lemmerling, 1999) will be used. To estimate the parameters, state space method proposed by Kung, Arun and Rao (1983), will be applied. Additional modelling techniques are required to include network nonlinearities. Typical network nonlinearity is of the saturation type and is exited only during extreme conditions
24
CHAPTER 2: PROJECT BACKGROUND
(e.g. saturation of power transformer, instrument current transformer, coil or MOV protection of series capacitors). The system can be modelled in all these instances, using a system of linear differential equations perturbed with nonlinearity (a quasi-linear system). In the modelling process, we use the Averaging approximation technique (only the first order approximation) referred to in engineering as ‘Describing Function’ method (Geleb and Velde, 1968). The signal model has the same structure as described above but the parameters are slow time-varying. The time-varying version of the structured total least squares method should be adopted here as discussed by Zivanovic (2001) and Jordaan and Zivanovic (2002). Higher harmonics associated with a nonlinearity will be treated as nuisance components and filtered out in the estimation process. A link between the time-varying signal model and the system model will be established using the averaging theory (Geleb and Velde, 1968). Other sources for the feature extraction process are binary disturbance records (opening and closing of the various contacts associated with the recorded event) and expert knowledge about disturbances provided by ESKOM specialists. Zivanovic (1999) also proposed a frequency estimation algorithm based on a local polynomial approximation technique for use in the feature vector construction. The semiparametric estimation approach for the feature vector construction is depicted in Figure 2.4. The input data in Figure 2.4 are the segmented disturbance signal as explained in Section 2.6.2.
CHAPTER 2: PROJECT BACKGROUND
25
FIGURE 2.4: Feature Vector construction using the Semi-Parametric Approach.
2.6.4 Pattern Matching and Disturbance Recognition The main objectives of this step are as follows: •
The best set of features (obtained from the parameter identification process and other sources) should be selected to represent disturbances.
•
The classification algorithm and training using recorded data and human expertise should be adopted.
•
The simple demonstration prototypes of the following analysis algorithms will be developed: fault location algorithm, an analysis for resonant conditions in the system with trapped energy, protection and control system performance analysis, electromechanical oscillation and system stability analysis.
26
CHAPTER 2: PROJECT BACKGROUND
Development in this step will rely on the available techniques from the statistical learning field (Vapnik, 1995). Unsupervised learning methods, such as, principal components and curves (Hastie, Tibsirani & Friedman, 2001) will be used for the selection and analysis of features. The classification task should deal with mixed data type (continuous values, breaker open and close, expert knowledge, etc). Keller, Henze and Zivanovic (2005) proposed a Support Vector Machine (SVM) based fault classifier. SVM is a training algorithm for learning classification and regression rules from data, first introduced by Vapnik in the 1960s (Vapnik, 1995). The SVM-based fault classifier separates ground faults and non-ground faults. Detailed mathematical description of the fault classification system can be found in the practical training report by Henze (2005). The input feature space is two-dimensional and consists of the magnitude of the 50 Hz component of the neutral current (X1) and the zero sequence voltage (X2). Typically 10 ground faults and 10 phase-to-phase faults are used as training data (Keller, Henze & Zivanovic, 2005). For the output, Yi = 1 means a ground fault and Yi = –1 a non-ground fault. After training the SVM classifier, another data set is used to test it. Typical test data set consists also of 10 ground and 10 non-ground faults (Keller, Henze & Zivanovic, 2005). Figure 2.5 shows the algorithm of the SVM fault classifier as a flowchart. The typical faults in power systems are three-phase short circuit (L-L-L), three-phase-toground (L-L-L-G), line-line short circuit (L-L), single-line-to-ground (L-G) and double line-to-ground (L-L-G) faults, where the terms ‘L’ and ‘G’ refer to ‘Line’ and ‘Ground’ respectively. The fault classifier system is directed mainly towards ground faults and singlephase-to- ground faults. It further classifies phase-to-phase, three-phase, phase-to-phase-
CHAPTER 2: PROJECT BACKGROUND
27
to-ground faults (See Figure 2.5). Voltage and current recordings for all the three phases as well as the neutral current are utilised for that.
FIGURE 2.5: Algorithm Flowchart of the SVM-based Fault Classifier.
2.7
Summary
In this chapter we have provided the necessary background information on the entire Automatic Disturbance Recognition in Electrical Power Systems Project. Detailed discussions about the background and motivation, the existing system, its shortcomings, and a brief overview of the proposed automated system were presented. We also provided introductory information about the subsystem blocks of the proposed automated systems. These include disturbance recordings from the digital fault recorders, abrupt change detection-based segmentation, feature vector construction, pattern matching and disturbance recognition.
28
CHAPTER 3: ABRUPT CHANGE DETECTION
CHAPTER 3 ABRUPT CHANGE DETECTION
This chapter forms the literature section of the thesis and contains the theoretical content necessary to understand the project and the need for it with reference to existing research. Following the work of Ukil and Zivanovic (2005a), a comparative detailed discussion of the existing and research studies on abrupt change detection relevant to this project and thesis is presented in this chapter which forms the basis of the rest of topics in the thesis.
3.1
Introduction
The detection of abrupt changes in signal characteristics is a widely studied topic with various approaches (Basseville, 1988; Basseville & Nikoforov, 1993). It plays a significant role in Failure Detection and Isolation (FDI) systems and the segmentation of signals in recognition-oriented signal processing. A possible approach to recognition-oriented signal processing involves the use of an automatic segmentation of the signal based on abrupt change detection as the first processing step. A segmentation algorithm splits the signal into homogeneous segments, the lengths of which are adapted to the local characteristics of the analysed signal which can be achieved either on-line or off-line (Basseville & Nikoforov, 1993).
3.2
Various Technologies
There are various algorithms to perform abrupt change detection, for example:
CHAPTER 3: ABRUPT CHANGE DETECTION
•
Limit Checking Detectors and Shewhart Control Charts (Shewhart, 1931)
•
Geometric Moving Average Control Charts (Roberts, 1959)
•
Finite Moving Average Control Charts (Laï, 1974)
•
Filtered Derivative Algorithms (Basseville et al., 1981)
•
Generalised Likelihood Ratio (GLR) Algorithm (Lorden, 1971)
•
Cumulative Sum (CUSUM) Algorithm (Page, 1954)
•
Bayesian-theory type detection (Girshick & Rubin, 1952).
29
The focus of this study is on power system fault and disturbance signals. Many of these signals are quasi-stationary, that is, the signals are composed of segments of stationary behaviour with abrupt changes in their characteristics in the transitions between different segments. It is imperative to find the time-instants in which the changes occur and to develop models for the different segments during which the system does not change. Of the various approaches, some depend on the signal statistical model, for example, the Generalised Likelihood Ratio (GLR) techniques and the Bayesian-theory type technique (Basseville & Nikoforov, 1993) where the signal under consideration can be assumed to be associated with a parametric system model. In other cases, where it is substantially difficult to estimate an accurate model, model-free approaches are proposed namely, machinelearning-based algorithms (Vapnik, 1995) and support-vector-based techniques (Gretton & Desobry, 2003). To accomplish the abrupt change detection, hence segmentation of the fault and disturbance signals, the following categorised approaches are considered. •
Simple methods o Superimposed Current Quantities o Linear Prediction Error Filter
30
CHAPTER 3: ABRUPT CHANGE DETECTION
o Adaptive Whitening Filter, •
Linear Model-based approach o Additive Spectral Changes o Autoregressive (AR) Modelling and Joint Segmentation o State-Space Modelling,
•
Model-free approach o Support Vector Machines,
•
Nonparametric approach o Discrete Fourier Transform o Wavelet Transform.
3.3
Problem Statement
Assuming a parametric system model, that is, when some inherent signal model parameter is changing, we can consider a quasi-stationary sequence of k independent observations x , with a d-dimensional parameter vector θ which describes the properties of the observations. Before the unknown change time t 0 , the parameter θ is equal to θ 0 , while after the change, it is equal to θ1 ≠ θ 0 . At this stage, two tasks are necessary: to detect the change time-instant t 0 and to estimate the corresponding parameter vectors θ 0 and θ1 . With the primary focus on detecting the change time-instant t 0 , it is useful to consider t 0 a random unknown variable with unknown distribution (Basseville & Nikoforov, 1993). In the Bayesian approach, the parameter vector θ is considered to be a random variable with a certain prior probability distribution, that is, before the observations have been
CHAPTER 3: ABRUPT CHANGE DETECTION
31
made (Ljung, 1999). The output and input data, y t and u t will be correlated with the θ [the superscript indicates the data record up to time t ]. At time t , we want to determine the posterior probability density function (PDF) for θ , that is, p (θ | y t , u t ) . From the posterior PDF different estimates of θ , θˆ(t ) can be estimated (Ljung & Söderström,
1986). For example: •
the mean of the posterior distribution, that is, the conditional expectation:
θˆ(t ) = E (θ | y t , u t ) , •
(3.1)
the value for which the PDF attains its maximum, in other words, the most likely value, which is known as maximum a posteriori (MAP) estimate (Ljung, 1999).
3.4
Abrupt Change Detection: Simple Methods
In this section we describe the simple methods for abrupt change detection in the power system fault and disturbance signals. The methods discussed here are: •
Superimposed Current Quantities
•
The Linear Prediction Error Filter
•
The Adaptive Whitening Filter.
3.4.1 Superimposed Current Quantities The superimposed current quantities technique (Stokes-Waller, 1998) is used to detect the signal peak values in the existing ESKOM semi-automated disturbance analysis software. Superimposed current quantities (∆I ) are created by calculating the difference of two
sample points, exactly one cycle (50 Hz system) apart. Constant system conditions will
32
CHAPTER 3: ABRUPT CHANGE DETECTION
produce no superimposed current output while any changes in network conditions, as on fault occurrence, produces a superimposed current quantity (Stokes-Waller, 1998). Figure 3.1 illustrates the principle of superimposed current quantities.
FIGURE 3.1: The Superimposed Current Quantities technique which calculates the difference of two sample points, exactly one cycle (50 Hz system) apart.
Figure 3.2 (Stokes-Waller, 1998) illustrates the way superimposed current quantities are used in the existing semi-automated analysis software for fault detection at ESKOM. Variables K1 ,....., K 6 shown in Figure 3.2 are user defined.
CHAPTER 3: ABRUPT CHANGE DETECTION
33
FIGURE 3.2: Fault Detection Logic using Superimposed Current Quantities.
3.4.2 The Linear Prediction Error Filter The optimal Linear Prediction Error (LPE) filter for transient monitoring perfectly decorrelates the signal, leaving only white noise (a whitening filter). The finite impulse response (FIR) whitening filter can be defined by, A( z −1 ) = 1 − z − C ,
(3.2)
where A indicates the finite impulse response of the whitening filter, z indicates the Ztransform and C is the rounded number of samples per cycle at the nominal network frequency. f C = round S f fund
,
(3.3)
where f S is the sampling frequency and f fund is the fundamental frequency (C = 50 at 50 Hz for a sampling frequency of 2.5 kHz). This filter compares the current signal value xk with the value approximately one cycle before, xk −C . The zeroes of this filter are the harmonics of the fundamental signal, up to
34
CHAPTER 3: ABRUPT CHANGE DETECTION
the Nyquist frequency, they are evenly spaced on the unit cycle. This whitening filter is known as the Fourier filter (Philippot, 1996). When the true fundamental frequency is close to f S C , but not exactly equal to this quantity, the zeroes of the whitening filter move away from the true harmonics. However, with additional computation, it is still possible to force a filter zero at the estimated fundamental frequency. This leads to the adjusted Fourier filter in case of a noninteger frequency ratio (Philippot, 1996). Using the same definition of C as before, let
A( z −1 ) = 1 + aC −1 z − C +1 + aC z −C .
(3.4)
From the condition: z C + aC −1 z + aC = 0 for z = e ± jω0T ,
(3.5)
we obtain the following formulae for computing aC −1 and aC .
aC −1 =
− sin(ω0TS C ) , sin(ω0TS )
aC = − cos(ω 0TS C ) − aC −1 cos(ω 0TS ) .
(3.6)
(3.7)
In the above cases, we consider the linear prediction. For a sampled signal {x k } , we can construct an error prediction {ε k } based on an extrapolation of x,
ε k = xk − xˆ k .
(3.8)
CHAPTER 3: ABRUPT CHANGE DETECTION
35
After launching the filter, an estimate of the error variance can be constructed by processing S successive samples,
σˆ S (k ) =
1 S −1 2 ∑ε k− j . S j =0
(3.9)
However, the described Fourier filters operate best in steady-state conditions. Trained
LPE filters are far more effective in transient conditions (Philippot, 1996). Given a data window covering a typical system state and containing the transients of interest, we can construct the trained LPE filters of any length by least squares fitting (e.g. by using the forward-backward technique). Hence the suggested selection of LPE filters as proposed in (Philippot, 1996) is
to use the Fourier filter in steady-state conditions
to use a short trained whitening filter just after a trigger
to switch from the short filter to a one-cycle trained whitening filter as soon as the error variance estimate of the second filter has been initialised.
LPE-based transient monitoring of a current signal is shown in Figure 3.3.
36
CHAPTER 3: ABRUPT CHANGE DETECTION
FIGURE 3.3: LPE-based Transient Monitoring of a Current Signal.
3.4.3 The Adaptive Whitening Filter The ideal prediction filter for transient monitoring will perfectly decorrelate the signal, leaving only white noise (whitening filter). In practice, mathematical simplicity and computation speed are key factors in the choice of a prediction filter. In the previous section we described the classical Fourier filter and the adjusted Fourier filter for this purpose. The frequency dependent coefficients of the adjusted Fourier filter ( aC −1 and aC ) are given by (3.6) and (3.7). The adaptive whitening filter is based on the adjusted Fourier filter taking
CHAPTER 3: ABRUPT CHANGE DETECTION
37
into account that the output of the filter is at a minimum when its coefficients are well adapted; a least mean square (LMS) estimate of aC −1 and aC is performed to minimise the output; with the constraint of exactly filtering the dc component of the filter (Wiot, 2004). For a sampled signal {x k } , we can construct an error prediction {ε k } so that
ε k = x k − α k x k − N +1 − β k x k − N ,
(3.10)
where α and β are equal to aC −1 and aC , and N is equal to C as given by (3.3). An LMS method is then used to estimate both α and β coefficients with the constraint that the sum of all coefficients is null. In other words, the dc component is perfectly filtered out (Wiot, 2004). This gives n
n
n
i =0
i =0
i =0
∑ ε k2 = ∑ ( xk −i − α k xk − N +1−i − β k xk − N −i ) 2 − ∑ λ (1 −α k − β k ) 2 .
(3.11)
The minimisation conditions are given by d n 2 ∑ε k = 0 , dα k i =0
(3.12)
d n 2 ∑ε k = 0 , dβ k i =0
(3.13)
d n 2 ∑ε k = 0 . dλ i =0
(3.14)
From (3.11) we obtain the solution of the linear equation system as
38
CHAPTER 3: ABRUPT CHANGE DETECTION
n
αk =
∑ ( xk −i − xk − N −i )( xk − N +1−i − xk − N −i ) i =0
n
∑ (x i =0
k − N +1−i
− x k − N −i )
2
n
=
∑I i =0
k −i
n
∑K i =0
K k −i , 2 k −i
(3.15)
βk = 1−αk . The recursive form of the adaptive whitening filter is given as follows: Starting from the general relation we get n
αk =
∑I i =0
n
k −i
∑K i =0
K k −i 2 k −i
=
γk . σk
(3.16)
Both numerator and denominator are easily updated
γ k +1 = γ k + I k +1 K k +1 − I k −n K k −n σ k +1 = σ k + K k2+1 − K k2−n
.
(3.17)
The computation of the β coefficient can be avoided by replacing it by (1- α ).
ε k = xk − α k xk − N +1 − β k xk − N ε k = xk − α k xk − N +1 − (1 − α k ) xk − N = ( xk − xk − N ) − α k ( xk − N +1 − xk − N ) = I k − α k K k
(3.18)
Hence the complete recursive equations of the Adaptive Whitening filter are as follows:
CHAPTER 3: ABRUPT CHANGE DETECTION
39
γ k = γ k −1 + I k K k − I k −n −1 K k −n−1 , σ k = σ k −1 + K k2 − K k2−n−1 , γ αk = k , σk
(3.19)
I k = xk − xk − N , K k = xk − N +1 − xk − N ,
ε k = Ik −αk Kk . An initialisation step is required and is given by n
γ k 0 = ∑ I k 0 −i K k 0 − i , i =0 n
σ k0 = ∑ K i =0
(3.20) 2 k 0 −i
.
The performance of the adaptive whitening filter, termed ‘ADAMO’ by Wiot (2004), in comparison with long and short Fourier filter is shown in Figure 3.4.
FIGURE 3.4: Adaptive Whitening Filter-based Transient Monitoring.
40
CHAPTER 3: ABRUPT CHANGE DETECTION
3.5
Abrupt Change Detection: The Linear Model based approach
In this section we describe the linear model-based approach for abrupt change detection in the power system fault and disturbance signals. The following methods are discussed here: •
Additive Spectral Changes
•
Autoregressive (AR) Modelling and Joint Segmentation
•
State-Space Modelling.
3.5.1 Additive Spectral Changes In this section we refer to more general and difficult cases in which changes occur in the spectral characteristics of the signal or system. In other words, these changes act as multiplicative changes in the transfer function. Here all the changes affect the dynamics of the system itself, making the problem more complex. For simplicity, we assume that the dynamics are modelled by the parameterised transfer function Tθ , and that the change in the dynamics is summarised as a change in the parameter θ from θ 0 to θ1 (Basseville & Nikoforov, 1993). Now we describe the spectral changes in an Autoregressive Moving Average (ARMA) model such as Yk = ∑i =1 AiYk −i + ∑ j =0 B jVk − j , p
q
(3.21)
where (Vk ) k is a white noise sequence with covariance matrix R and B0 = I . The spectral changes are changes in the shape of the spectrum of the observations. In other words, these changes are changes in the matrix coefficients Ai and B j :
CHAPTER 3: ABRUPT CHANGE DETECTION
Yk =
∑ ∑
p
i =1
Ai0Yk −i + ∑ j =0 B 0j Vk − j if k < t 0
41
q
A1Y + ∑ j =0 B1jVk − j if k ≥ t 0 i =1 i k −i p
q
.
(3.22)
The next step is to describe the spectral changes in a state-space model for which we consider the linear dynamic system described by the state-space representation of the observed signals (Yk ) k : X k +1 = F X k + GU k + Wk Yk = H X k + JU k + Vk
,
(3.23)
where the state X , the measurement Y , and the control U have dimensions n , r , and m , respectively, and where (Wk ) k and (Vk ) k are two independent Gaussian white noises, with covariance matrices Q and R , respectively.
The spectral changes can be modelled by changes in the state transition matrix F of the state-space model (3.23), in the following manner:
F=
F0 if k < t 0 F1 if k ≥ t o
.
(3.24)
These changes affect the denominator of the transfer function. A traditional approach to failure detection involves considering the fact that the design of the detection algorithm basically comprises two steps: (1) the generation of residuals, namely of artificial measurements that reflect the changes of interest; for example, these signals are ideally close to zero before the change and nonzero after
CHAPTER 3: ABRUPT CHANGE DETECTION
42
(2) the design of decision rules based upon these residuals. The key technique for detecting additive changes consists of achieving step 1 mentioned above in such a way that the change detection problem that results from the transformation of the observations into the residuals is exactly the problem of detecting changes in the mean value of a multidimensional Gaussian process. For the case of spectral changes, the key technique for achieving step 1 is neither the transformation to innovations nor the transformation to parity checks, which are no longer sufficient statistics in this case. The sufficient statistics here is the likelihood ratio (Basseville & Nikoforov, 1993).
3.5.2 Autoregressive (AR) Modelling and Joint Segmentation Autoregressive (AR) modelling is a well-known technique for stationary signal analysis. However, in practice, signals are generally nonstationary; a frequent class of them, referred to as quasi-stationary is characterised by abrupt changes between the stationary segments with different statistical properties. A complete solution to this problem involves finding the correct number of models (or stationary segments), their orders and parameters, and the transitions among them (Caballero, Prieto and Vidal, 1997). We consider that the observations are generated by switching between M different AR models of orders p1 ,...., p M , and coefficients a j = (a1, j ,...., a p j , j ) , that is:
M
x[n] = ∑ t j [n]x j [n] , j =1
where t j [n] selects the samples generated by model j ,
(3.25)
CHAPTER 3: ABRUPT CHANGE DETECTION
1, if x[n] is generated by model j
t j [ n] =
0, otherwise.
.
43
(3.26)
On the other hand, the output at time instant n for model j is given by
pj
x j [n] = −∑ ai , j x[n − i ] + e j [n]
j = 1,..., M ,
(3.27)
i =1
where e j [n] is a zero mean uncorrelated Gaussian noise with constant variance σ 2j . Equation (3.27) is generated using p j past samples of the observations: when a transition occurs these past values are used as the initial condition for the new model after change (Caballero, Prieto & Vidal, 1997). The method of generating the quasi-stationary signal is depicted in Figure 3.5, where A j ( z ) = 1 + ∑i =j1 ai , j z − i . p
FIGURE 3.5: Method and Model of generating a Quasi-stationary Signal.
The problem can be stated as: given the number of models M , their orders ( p1 ,...., p M ) , and the vector of observations x = ( x[0],...., x[ N - 1])T , first, determine the boundaries between segments; and second, find the best model for each segment, that is the, joint segmentation of the signal.
44
CHAPTER 3: ABRUPT CHANGE DETECTION
The objective is to find the likelihood function for the unknown parameters. It is needed to group the coefficients of the AR models, the variances of the noise sequences and the transition sequences in vectors a = (a1 ,...., a M )T , σ = (σ 12 ,....,σ M2 ) T and t = (t 1 ,...., t M ) , respectively (Caballero, Prieto and Vidal, 1995). Assuming independence among the stationary segments, we have the log likelihood function for the unknown parameters given by 2
pj M N −1 1 L(x; t, a,σ ) = − ∑ x [ n ] + t [ n ] a x [ n − i ] − ∑ ∑ ∑ plog) σ n . j i, j 2 n =max( n = max( p j ) 2σ n j =1 i =1 j N −1
(3.28)
A direct maximisation of (3.28) cannot be carried out in practice because of the discrete nature of the transition sequences t . To resolve this difficulty, a suboptimal procedure based on the Expectation-Maximisation (EM) algorithm (Caballero, Prieto and Vidal, 1995) is applied, which increases the likelihood of the obtained estimates, iteration by iteration, until a local maximum is reached. At each iteration the transitions are estimated as the posterior probabilities that a sample was generated by a given model (E-step); then, the new set of models is obtained by solving a least-squares problem (M-step). The prediction errors obtained for a given set of models are used to estimate the probabilities that a sample was generated by a particular model (E-step) and, using these estimates, a new set of optimal models is obtained (M-step). An adequate choice of the initial values obtained via competitive modelling achieves fast convergence and improved estimates. A practical result using a quasi-stationary process is shown in Figure 3.6. The quasistationary process is shown in plot (a); the posterior probabilities for model 2 at iteration 1 and 7 are shown in plot (b) and plot (c) respectively; and the evolution of the log likelihood function is shown in plot (d).
CHAPTER 3: ABRUPT CHANGE DETECTION
(a) A typical quasi-stationary process
45
(b) Posterior probabilities for model 2, iteration 1
(c)
(d)
Posterior probabilities for model 2, iteration 7
Evolution of the log likelihood function
FIGURE 3.6: Practical Application Result for Joint Segmentation using the EM Algorithm.
3.5.3 State-Space Modelling Mathematically, we can consider a sequence of observations depending only upon one scalar parameter θ . Before the unknown change time t 0 , the parameter θ is equal to θ 0 , and after the change it is equal to θ1 ≠ θ 0 . The problems are then to detect and estimate the changes in the parameter and the change time-instant t 0 (Basseville & Nikoforov, 1993).
CHAPTER 3: ABRUPT CHANGE DETECTION
46
Since power system has a deterministic system model, we can apply the model-based approach for abrupt change detection. Ljung and Söderström (1986) propose linear and general state-space modelling techniques. Linear state-space models can be solved by estimating the parameters using the well-known Kalman filters (Kalman, 1958). Important aspects of adaptive filtering using the Kalman filters are elucidated in the works of Jazwinski (1969), Franklin, Powell and Workman (1990) and Haykin (1996). To solve the linear state-space model, the Kalman filters are run in parallel, each of them corresponding to a particular assumption about when the system actually changed. The relative reliability of these assumed system behaviours is constantly judged, and unlikely hypotheses are replaced by new ones. General state-space models can be solved by the Extended Kalman filters (EKF) as studied by Ljung (1979). The relationship between the input u (t ) and the output y (t ) is described using the statespace model
x(t + 1) = Fx(t ) + Gu (t ), y (t ) = Hx(t ),
(3.29)
where u (t ) is an n-dimensional vector, and F, G, and H matrices of compatible dimensions. For linear state-space modelling, the parameter θ and the data are linearly related. Considering a general situation, (3.29) can be represented as
x(t + 1) = F (θ ) x(t ) + G (θ )u (t ) + w(t ), y (t ) = H (θ ) x(t ) + e(t ),
(3.30)
where {w(t )} and {e(t )} are sequences of independent random vectors with zero mean and
CHAPTER 3: ABRUPT CHANGE DETECTION
47
covariance matrices R1 (θ ) and R2 (θ ) respectively. The parameter θ can be estimated using the Kalman filters and the recursive least squares algorithm (Ljung & Söderström, 1986).
3.6
Abrupt Change Detection: The Model-free approach
Many abrupt change detection techniques depend on the knowledge of a signal statistical model. However, in some applications, because it may be difficult to design an accurate and tractable statistical model, model-free approaches need to be considered. In this section, a model-free, machine learning-based online algorithm for abrupt change detection in signals is described.
3.6.1 Support Vector Machines Support Vector Machine (SVM) is one of the relatively new and promising methods for learning separating functions in pattern recognition (classification) tasks, or for performing function estimation in regression problems. SVMs were originated from Statistical Learning Theory (SLT) by Vapnik (1995) at Bell Labs for “distribution free learning from data”. Desobry and Davy (2003) proposed a two-step algorithm for abrupt change detection in signals using SVM. First, informative descriptors (or vectors) localised in time, denoted by r xt , are extracted online from the signal. These can be cepstral coefficients computed on a
sliding window or short-time Fourier transforms (STFTs), and such like. Second, a kernelbased online stationarity index I (t ) is defined and computed in the descriptors space (or input space, denoted by X ), and geometrically defined in a feature space F induced by a Mercer kernel (Smola & Schölkopf, 2003) k (.,.) . Roughly, I (t ) is computed as follows: a
48
CHAPTER 3: ABRUPT CHANGE DETECTION
kernel k (.,.) is selected. At time t , a first ν -support vector (SV) novelty detector is
r r trained over the m1 last descriptors x1 = ( xt −m1 ,..., xt −1 ) , yielding a decision region R1 in r X . R1 is such that a vector x is considered to be similar to x1 iff x1 ∈ R1 . Next, a second
r r ν -SV novelty detector is trained over m2 future descriptors x2 = ( xt ,..., xt + m −1 ) , yielding 2
R2 . The regions R1 and R2 are representative of the probability density functions (pdfs)
which generated the sets x1 and x2 . Thus, comparing the geometries and locations of R1 and R2 is a robust way of comparing x1 and x2 . r r As described above, let us consider at time t , two subsets x1 = ( xt −m1 ,..., xt −1 ) and r r x2 = ( xt ,..., xt + m2 −1 ) of size m1 and m2 . Each of these subsets is used to train independently
ν -SV novelty detectors, yielding parameters w1 , ρ1 and w 2 , ρ 2 , and decision regions R1 and R2 . The idea underlying the abrupt changes detector is that a sudden change at time t r in the distribution of vectors xt may result in different locations/geometries of R1 and R2 .
An abrupt change corresponds to a large distance between the circle centres (w.r.t. their radii). In practice, at each time t , an index I (t ) is built which reflects the dissimilarity between x1 and x2 via a measure of the dissimilarity between R1 and R2 . The computation details of the updating of the parameters w 1 , ρ1 and w 2 , ρ 2 can be found in the work of Gretton and Desobry (2003). In the feature space, the shapes of the mapped regions R1 and R2 are simple (their boundaries are hypercircles C1 and C 2 ). A simple way to compare C1 and C 2 (i.e., R1 and R2 ), and to build I (t ) , involves considering an interclass (between R1 and R2 as separate
regions) / intraclass (within R1 or R2 ) ratio in F , such as
CHAPTER 3: ABRUPT CHANGE DETECTION
I (t ) =
distance between circle centers . radius of C1 + radius of C2
49
(3.31)
The algorithm is depicted in Figure 3.7. In Figure 3.7, in the feature space F , the training data are mapped on a hyperspace S with radius 1 and centre 0 . The ν -SV novelty detector related to x1 (resp. x2 ) yields a hypercircle C1 (resp. C 2 ) in the hyperplane W1 (resp. W2 ). C1 and C 2 are the (geometrical) centres of the regions occupied by the feature data,
and p1 and p 2 possess properties similar to mapped margin support vectors as they lie on the separating hyperplane.
FIGURE 3.7: A Geometrical depiction of the SVM-based Abrupt Change Detection Algorithm.
A practical application result on music signals using the SVM-based abrupt change detection is shown in Figure 3.8. In Figure 3.8, the music signal is shown in the top-plot and the corresponding SVM-based abrupt change detection index in bottom-plot. For a threshold
η equal to 0.85, all changes (dashed lines) are correctly detected, with only one false positive (circled).
CHAPTER 3: ABRUPT CHANGE DETECTION
50
FIGURE 3.8: An Application Result (Music Signals) for SVM-based Abrupt Change Detection.
3.7
Abrupt Change Detection: Nonparametric approach
In nonparametric approach, that is, when we do not assume any inherent signal model parameter that is changing, we consider the following methods: •
The Fourier Transform
•
The Wavelet Transform.
3.7.1 The Fourier Transform Fourier analysis is extremely useful for data analysis because it breaks down a signal into constituent sinusoids of different frequencies. For sampled vector data, Fourier analysis is performed using the Discrete Fourier transform (DFT) (Rabiner & Gold, 1975; Oppenheim & Schafer, 1989). The Fast Fourier transform (FFT) is an efficient algorithm for computing the DFT of a sequence; it is not a separate transform. It is particularly useful in areas
CHAPTER 3: ABRUPT CHANGE DETECTION
51
such as signal and image processing, where its uses range from filtering, convolution and frequency analysis to power spectrum estimation. For the input sequence x and its transformed version X (the discrete-time Fourier transform at equally spaced frequencies around the unit circle), the following two functions implement the relationships: N −1
X (k ) = ∑ x(n)W Nkn n =0
1 x ( n) = N
N −1
∑ X (k )W k =0
(3.32) − kn N
,
where
WN = e
2π − j N
.
(3.33)
The current fault analysis software in ESKOM uses a combination of superimposed current quantities and DFT (Stokes-Waller, 1998; Stokes-Waller & Keller, 1998). DFT is used to extract harmonically related voltage and current phasor magnitudes and angles from the recorded analogue signals (Stokes-Waller, 1998).
3.7.2 The Wavelet Transform The Wavelet transform (WT) is a mathematical tool, like Fourier transform for signal analysis. A wavelet is an oscillatory waveform of effectively limited duration that has an average value of zero. Fourier analysis consists of breaking up a signal into sine waves of various frequencies. Similarly, wavelet analysis is the breaking up of a signal into shifted and scaled versions of the original (or mother) wavelet. Figure 3.9 shows the basis func-
52
CHAPTER 3: ABRUPT CHANGE DETECTION
tions for Fourier transform (the sine wave) and WT (db10: Daubechies 10 [Daubechies, 1992] mother wavelet).
FIGURE 3.9: Basis Functions for Fourier Transform (sine wave) and Wavelet Transform (db10).
Fourier analysis does not provide good results for nonstationary signals, unlike the stationary signal, where the signal parameters change over the time, because in transforming the complete signal to the frequency domain, the time information is lost in Fourier analysis. This deficiency in Fourier analysis can be overcome to some extent by analysing a small section of the signal at a time - a technique called windowing the signal. This leads to an analysis technique called Short-Time Fourier Transform (STFT). But the drawback in STFT is that the size of the time-window is same for all frequencies. Wavelet analysis overcomes this deficiency by allowing a windowing technique with variable-sized regions, in other words, wavelet analysis allows the use of long time intervals where we want more precise low-frequency information, and shorter regions where we want high-frequency information. Figure 3.10 shows the time-based (Shannon), frequency-based (Fourier), STFT (Gabor) and wavelet views of signal analysis. Detailed mathematical descriptions of WT can be referred to in the seminal works of Daubechies (1992) and Mallat (1998).
CHAPTER 3: ABRUPT CHANGE DETECTION
53
FIGURE 3.10: Time, Frequency, STFT and Wavelet views of Signal Analysis.
In the past ten years, the application of the wavelet transform in the field of power systems has received particular attention. In power systems signal analysis, time-frequency resolution is needed, which states the reason for using the wavelet transform because it provides a local representation (both in time and frequency) of a given signal. It is different from the Fourier transform which provides a global representation of a signal. Examples of the use of the wavelet transform in power systems can be found in the works of Robertson, Camps and Mayer (1994), Robertson, Camps, Mayer et al. (1996), Wilkinson and Cox (1996). Santoso, Powers and Grady (1994) and Santoso, Powers, Grady et al. (1996) used the wavelet transform for power quality assessment. Ribeiro (1994) proposed the use of the wavelet transform for analysing the harmonics in power systems. Pillay and Bhattacharjee (1996) used the wavelet transform to model short-term power system disturbances.
54
3.8
CHAPTER 3: ABRUPT CHANGE DETECTION
Summary
We have discussed comparatively about various abrupt change detection techniques with reference to the topic of this thesis. We broadly classified the relevant techniques such as simple methods, the linear model-based approach, the model-free approach, nonparametric approach. We also discussed in details about these categories and the specific techniques included therein. This chapter provides the basic theoretical background for the rest of the thesis.
CHAPTER 4: THE RECURSIVE IDENTIFICATION METHOD
55
CHAPTER 4 THE RECURSIVE IDENTIFICATION METHOD
This chapter describes the application of the recursive identification technique to detect the abrupt changes in the disturbance recordings. The recursive identification technique uses M parallel Kalman filters. The main focus has been to estimate exactly the time-instants of the change in signal model parameters during the pre-fault condition and following the events such as the initiation of fault, the circuit-breaker opening, auto-reclosure of the circuit-breakers and such like. Ukil and Zivanovic (2005f) provide a detailed discussion of the topic.
4.1
Introduction
Many of these signals are quasi-stationary in the sense that they are composed of segments of stationary behaviour with abrupt changes in their characteristics in the transitions between different segments. It is important to find the time-instants at which the changes occur and to develop models for the different segments during which the system does not change. Some of the approaches depend on the signal statistical model namely, Generalised Likelihood Ratio (GLR) techniques and the Bayesian-theory type technique (Basseville & Nikoforov, 1993) where the signal under consideration can be assumed to be associated with a parametric system model. Typically the first processing step in a recognitionoriented signal processing can be automatic segmentation of the signal which can be done on-line or off-line (Basseville & Nikoforov, 1993).
56
CHAPTER 4: THE RECURSIVE IDENTIFICATION METHOD
Assuming a parametric system model (discussed in chapter 3.3; mentioned here for the sake of easy and continuous following of the current discussion), we consider a quasistationary sequence of k no. of observations x , with a d -dimensional parameter vector
θ which describes the properties of the observations. Before the unknown change timeinstant t 0 , the parameter θ is equal to θ 0 , and after the change it is equal to θ1 ≠ θ 0 . Two tasks are then: to detect the change time-instant t 0 and to estimate correspondingly the parameter vectors θ 0 and θ1 . With the primary focus on detecting the change time-instant t 0 , it is useful to consider t 0 a random unknown variable with unknown distribution (Basseville & Nikoforov, 1993). In the Bayesian approach the parameter vector θ is considered to be a random variable with a certain prior probability distribution i.e., before the observations have been made (Ljung, 1999). The output and input data, y t and u t will be correlated with the θ [the superscript indicates data record up to time t ]. At time t , we want to determine the posterior probability density function (PDF) for θ i.e., p(θ | yt ,ut ) . From the posterior PDF different estimates of θ , θˆ ( t ) can be estimated (Ljung & Söderström, 1986) like:
the mean of the posterior distribution i.e., the conditional expectation:
θˆ ( t ) = E (θ | y t , u t ) ,
(4.1)
the value for which the PDF attains its maximum i.e., the most likely value, which is known as maximum a posteriori (MAP) estimate (Ljung, 1999).
θˆMAP (t ) = arg max p (θ | y t , u t ) . θ
(4.2)
CHAPTER 4: THE RECURSIVE IDENTIFICATION METHOD
4.2
57
Recursive Identification
The recursive estimation algorithm builds segmented Autoregressive (AR), Autoregressive Moving Average (ARMA) type models from the observations. The polynomial form of the model is A(q) y (t ) = B (q)u (t ) + e(t ) ,
(4.3)
where,
u (t ) : Input, y (t ) : Output, e(t ) : Gaussian White Noise A(q) , B (q ) : Polynomial Transfer Functions q is the backward shift (or delay) operator
q −1 y (t ) = y (t − 1), A(q) = 1 + a1 q −1 + ...... + a n q − n ,
(4.4)
B(q) = b1 + b2 q −1 + ...... + bm q − m . The model in (4.3) and (4.4) describes the dynamic relationship between the input and output signals. It can be expressed in terms of the parameter vector (Ljung & Söderström, 1986), θ T = (a1 .....an ; b1 .....bm ) . We also introduce the vector of lagged input-output data (Ljung & Söderström, 1986),
ϕ T (t ) = (− y (t − 1)... − y (t − n); u (t − 1)...u (t − m) ) .
(4.5)
Then (4.3) can be rewritten as y (t ) = θ T ϕ (t ) + e(t ) .
(4.6)
58
CHAPTER 4: THE RECURSIVE IDENTIFICATION METHOD
A typical recursive identification algorithm is:
θˆ ( t ) = θˆ ( t − 1) + K ( t ) [ y ( t ) − yˆ ( t ) ], (4.7)
θˆ ( 0 ) = θ 0 , where θˆ(t) is the estimate of the parameter vector θ at time t , and y (t ) is the observed output at time t . yˆ (t ) is a prediction of the value y (t ) based on observations up to time t − 1 and also based on the current model (and possibly also earlier ones) at time t − 1 , that
is, the vector ϕ (t ) is a function of yt −1 , ut −1 . The gain K (t ) determines in what way the current prediction error y (t ) − yˆ (t ) affects the update of the parameter estimate (Ukil & Zivanovic, 2005f). e(t ) in (4.3) and (4.6) is a sequence of independent Gaussian variables with E{ e(t ) }= 0 and E{ e2 (t) }= R2 (t) . Then the posterior distribution p(θ | yt , ut ) is also Gaussian with mean θˆ ( t ) and covariance matrix P (t ) . θˆ ( t ) can be estimated using (4.7). P (t ) can be estimated according to (see Ljung and Söderström (1986) for the proof)
P (t ) = P(t − 1) −
P(t − 1)ϕ (t )ϕ T (t ) P(t − 1) , R2 (t ) + ϕ T (t ) P(t − 1)ϕ (t )
(4.8)
P (0) = P0 . An optimal choice of the gain K (t ) can then be computed using the Kalman filter (Ljung & Söderström, 1986)
K (t ) =
P(t − 1)ϕ (t ) . R2 (t ) + ϕ T (t ) P(t − 1)ϕ (t )
(4.9)
It is interesting to note that the model as in (4.9) can be seen as a linear state-space model (Ljung & Söderström, 1986)
CHAPTER 4: THE RECURSIVE IDENTIFICATION METHOD
θ (t + 1) = θ (t ) + w(t ),
59
(4.10)
y (t ) = ϕ T (t )θ (t ) + e(t ).
The segmentation algorithm used here to estimate the final segmented model is based on the model description as in (4.10), where the disturbance term w(t ) is zero most of the time, but occasionally it abruptly changes the system parameter vector θ(t) . If the disturbance w(t ) in (4.10) is Gaussian with covariance matrix R1 (t) , the best estimate of θ is given by the Kalman filter. If w(t ) is not Gaussian, the Kalman filter does not provide the optimal solution and the problem (4.10) becomes a nonlinear filtering problem. The recursive identification algorithm as used in the segmentation of the power system fault recording signals solves the nonlinear filtering problem using an approach based on finiteGaussian sum approximation (Andersson, 1985). We suppose that the posterior distribution of θ (t ) given y t −1 (i.e., all old y ’s up to and including y (t − 1) ) can be approximated with a sum of M Gaussian density functions (Andersson, 1985). For the n -dimensional Gaussian probability density we introduce the notation Gn ( x, m, P) = (2π ) − n / 2 (det P) −1 / 2 exp(− 12 ( x − m) T P −1 ( x − m)) .
(4.11)
The posterior density function for θ (t ) can then be written
M
p[θ (t ) | y t −1 ] = ∑ α i (t )Gn (θ (t ),θ i (t ), Pi (t )) ,
(4.12)
i =1
where α i is the posterior probability,
M
∑α (t ) and θ (t ) i =1
i
i
and Pi (t ) are the mean vectors
and covariance matrices respectively of the different Gaussian distributions at time t (Andersson, 1985).
60
CHAPTER 4: THE RECURSIVE IDENTIFICATION METHOD
This gives the recursive identification algorithm termed as adaptive forgetting through multiple models (AFMM) (Andersson, 1985) that can be summarised as: For i = 1,2,..., M
Pi (t − 1)ϕ (t )ϕ T (t ) Pi (t − 1) Pi (t ) = Pi (t − 1) − , R2 (t ) + ϕ T (t ) Pi (t − 1)ϕ (t )
(4.13)
ε i (t ) = y (t ) − ϕ T (t )θ i (t − 1) ,
(4.14)
θ i (t ) = θ i (t − 1) +
α i (t ) =
1 Pi (t )ϕ (t )ε i (t ) , R2 (t )
(4.15)
α i (t )
1 ε i2 (t ) , (4.16) exp − T ( R2 (t ) + ϕ T (t ) Pi (t − 1)ϕ (t )) 2 ( R2 (t ) + ϕ (t ) Pi (t − 1)ϕ (t ))
imin = arg min α i (t ) ,
(4.17)
imax = arg max α i (t ) ,
(4.18)
Pimin (t ) = R1 (t ) + Pimax (t ) ,
(4.19)
θ i (t ) = θ i (t ) .
(4.20)
i
i
min
max
We can assign the minimum posterior probability with a probability of jump denoted as q (Andersson, 1985),
αi
min
=q,
(4.21)
CHAPTER 4: THE RECURSIVE IDENTIFICATION METHOD
α i (t ) =
1 M
∑α k =1
k
α i (t ) .
61
(4.22)
(t )
The estimate θˆ(t ) , of θ (t ) then becomes
M
θˆ(t ) = ∑ α i (t )θ i (t ) .
(4.23)
i =1
It is noted that this algorithm consists of M recursive least squares algorithms (or Kalman filters) operating in parallel, equations (4.13) to (4.15), and computations of posterior probabilities α i , (4.16), (4.21) and (4.22), together with administration for cutting steps, (4.17) to (4.20). The administration for cutting steps is to allow only the most likely component (the one with the largest α i ) and also, to entirely cut off the least likely component (the one with the smallest α i ). The major computational burden lies in (4.13). Except for the initial values of θ i , Pi and α i , four choices have to be made in using this method. These design variables (Andersson, 1985) are as follows: •
R1 (t ) covariance matrix for the jumps,
•
q
probability of a jump,
•
M
number of Gaussian distributions used to approximate the distribution of the parameter vector,
•
R2 (t ) equation error noise variance.
Since R1 (t ) and q are only used to realise the algorithm at every step, the choice of these is not too critical. The only requirement of M should be that it is chosen large enough.
62
CHAPTER 4: THE RECURSIVE IDENTIFICATION METHOD
The choice of R2 (t ) is probably more critical, since the value of this is used in all the Kalman filters. Moreover, the relative reliability of a measurement (which is important when judging whether a jump has occurred) is directly related to R2 (t ) .
4.3
Implementation Steps
The implementation (Ukil & Zivanovic, 2005f) steps for this recursive identification algorithm are described below. •
M models [see (4.3) and (4.6)] of fixed order are estimated in parallel, from which the final one will be the estimated model.
•
Each of these models is based on a particular assumption about when the system actually changed.
•
Each of the M models is estimated recursively using the Kalman filter-type estimation algorithm.
•
At each time-step, the posterior probability for each model is calculated.
•
The time varying estimate is formed by weighting together the M different models with weights equal to their posterior probability.
•
The lowest posterior probability models are rejected at each time step.
•
A new model is started, assuming that the system parameters have jumped.
•
The surviving model with the highest posterior probability is considered once all the observation data have been examined.
CHAPTER 4: THE RECURSIVE IDENTIFICATION METHOD
•
4.4
63
This model is tracked back and the change time-instants marked.
Application Results
The MATLAB® (Mathworks Inc., 2002a) algorithm implemented for the segmentation of the power system fault analysis signals was developed according to the above discussed recursive parameter estimation with the heuristic smoothing filtering operations. The recursive identification algorithm detects the change time-instants and also models the corresponding segments. The change time-instants should be indicated as unit impulses. After normalising the original fault signal using its mean value, recursive parameter estimation is applied on it to determine the change time-instants and model the corresponding segments. Heuristic smoothing filter operations are then applied on this segmented model to perform sequentially the following smoothing operations: •
Since the original fault signals are embedded with noise, the change time-instants obtained using the recursive identification algorithm show multiple results. Hence a threshold checking is done using the universal threshold technique of Donoho and Johnstone (Donoho & Johnstone, 1994), that is, the multiple spikes are tested to determine which values actually exceed a given threshold. The threshold T, is computed according to T = σ 2 log e N ,
(4.24)
where σ is the median absolute deviation of the segment model, (Donoho & Johnstone, 1994) for the fault signal and N is the number of samples in the seg-
64
CHAPTER 4: THE RECURSIVE IDENTIFICATION METHOD
ment model. Instead of standard deviation, median absolute deviation is used because the median is hardly influenced by a small fraction of extreme values. •
The segments in the power system fault analysis signals occur during the pre-fault condition and subsequent events such as fault initiation, circuit-breaker opening and reclosing. These events as well as the respective segments are predefined. Hence any larger number of segments possibly indicates transients, power swings, and such like. Estimation of the number of segment(s) is also performed and checked.
•
Based on the modelling of the segments, analysis is done to estimate the eventcritical change instants, rejecting others.
Furthermore, the filtering operation removes confusing multiple close-spikes and combines them into single unit impulse for clear indication. Any unwanted glitches which could otherwise indicate a false alarm for segmentation are also removed. Figures 4.1 and 4.2 show the application results for the disturbance signals obtained from ESKOM DFRs. In Figure 4.1, the upper plot shows the original DFR recording for the voltage during a phase-to-ground fault in the RED-Phase, sampled at a sampling frequency of 2.5 kHz. The middle plot shows the modelling of the system using the recursive identification technique. The lower plot shows the time-instants of the changes in the signal characteristics, marked by the impulse indicators, indicating the different signal segments owing to different events during the fault. For example, segment A indicates the pre-fault section and the fault inception, segment B the fault, segment C the opening of the circuit-breaker, segment D auto-reclosing of the circuit-breaker and system restore.
CHAPTER 4: THE RECURSIVE IDENTIFICATION METHOD
65
FIGURE 4.1: Segmentation of RED-Phase Voltage Signal using the Recursive Identification Method.
In Figure 4.2, the upper plot shows the original DFR recording for the current during a phase-to-ground fault in the RED-Phase, sampled at a sampling frequency of 2.5 kHz. The middle plot shows the modelling of the system using the recursive identification technique. The lower plot shows the time-instants of the changes in the signal characteristics, marked by the impulse indicators, indicating the different signal segments owing to different events during the fault. For example, segment A indicates the pre-fault section and the fault inception, segment B the fault, segment C the opening of the circuit-breaker, segment D auto-reclosing of the circuit-breaker and system restore.
66
CHAPTER 4: THE RECURSIVE IDENTIFICATION METHOD
FIGURE 4.2: Segmentation of RED-Phase Current Signal using the Recursive Identification Method.
4.5
Comments on Application Results
Following the discussion of the recursive identification algorithm and the application results, the following comments are appropriate here (Ukil & Zivanovic, 2005f). • Since the intended application is not meant for real-time analysis, computation time
is not a critical factor. However, the proposed algorithm for the abrupt change detection and signal segmentation took an average computation time of 5.585 seconds. An Intel® Celeron® 1.9 GHz, 256 MB RAM computer was used for all the application tests using MATLAB® (Mathworks Inc., 2002a). It should be noted that the complete automatic disturbance recognition and analysis tasks have to be per-
CHAPTER 4: THE RECURSIVE IDENTIFICATION METHOD
67
formed within five minutes of acquiring the fault signals, abrupt change detection and segmentation being the first step. • Two-hundred disturbance records were tested with the proposed algorithm and
95% accuracy in correct segmentation was achieved. However, in certain kinds of disturbances, such as those showing gradual resistive decay, not indicating marked abrupt changes, some false positives occur. However, in those cases, the focus is on correctly determining the fault inception time and the effective use of the smoothing filter reduces the false positives. • The proposed algorithm using the recursive identification technique is considerably
faster and more robust compared with the traditional peak value detection and superimposed current quantities algorithms (Stokes-Waller & Keller, 1998). Also, the automatic segmentation based on abrupt change detection using the recursive identification technique facilitates further signal processing and analysis in the subsequent stages of automatic disturbance recognition and analysis. It focuses on the different segments and helping to determine quickly parameters such as duration of the fault, etc, directly from the abrupt change detection-based segmentation itself, which cannot be done using the traditional peak value detection and superimposed current quantities algorithms (Stokes-Waller & Keller, 1998).
4.6
Summary
We have presented in this chapter the recursive identification technique used for detecting the abrupt changes in the disturbance signals. Since the power system has a deterministic system model, we can apply the model based approach for abrupt changes detection. The
68
CHAPTER 4: THE RECURSIVE IDENTIFICATION METHOD
Kalman filter algorithm was used to perform the recursive parameter estimation. Several Kalman filters that estimate parameters were run in parallel, each of them corresponding to a particular assumption about when the system actually changed. The relative reliability of these assumed system behaviours was constantly judged and unlikely hypotheses replaced by new ones. The segmentation algorithm modelled the system as AR, ARMA type. This was followed by the heuristic smoothing filtering operation. We obtained fairly satisfactory results from the MATLAB® implementation. Hence the recursive identification method was reasonably effective in detecting the abrupt changes in the disturbance signals by modelling and segmenting the different piecewise constant segments from various events during the fault.
CHAPTER 5: THE WAVELET TRANSFORM METHOD
69
CHAPTER 5 THE WAVELET TRANSFORM METHOD
This chapter describes the application of the wavelet transform to detect the abrupt changes in the disturbance recordings. The key idea is to decompose the fault signals into effective detailed and smoothed version using the multiresolution signal decomposition technique based on the discrete wavelet transform. Then we apply the threshold method on the decomposed signals to estimate the change time-instants, segmenting the fault signals. Ukil and Zivanovic (2005b) provide a detailed discussion of the topic.
5.1
Introduction
The wavelet transform is particularly suitable for the power system disturbance and fault signals which may not be periodic and may contain both sinusoidal and impulse components. In addition, time-frequency resolution is needed for the power system fault analysis. This is another reason for using the wavelet transform because it provides a local representation (both in time and frequency) of a given signal, unlike the Fourier transform which provides a global representation of a signal. In particular, the ability of the wavelets to focus on short intervals for high-frequency components and long intervals for low-frequency components improves the decomposition of the fault signals into finer and detailed scales, facilitating further effective signal processing and analysis. In this proposed method, the wavelet transform is used to transform the original fault signal into finer wavelet scales, followed by a progressive search for the largest wavelet coef-
70
CHAPTER 5: THE WAVELET TRANSFORM METHOD
ficients on that scale (Craigmile & Percival, 2000). Large wavelet coefficients that are colocated in time across different scales provide estimates of the changes in the signal parameter. The change time-instants can be estimated by the time-instants when the wavelet coefficients exceed the first-order approximation of the universal threshold of Donoho and Johnstone (Donoho & Johnstone, 1994).
5.2
Wavelet Transform Analysis
Introduction about the wavelet transform (WT) has been discussed in chapter 3, section 3.7.2. While detailed mathematical descriptions of WT can be referred to in the seminal works of Daubechies (1992), Chui (1992), Meyer (1993), Strang and Nguyen (1996), Mallat (1998), a brief mathematical summary of WT needed for this thesis is provided in the following sections.
5.2.1 The Continuous Wavelet Transform The Continuous Wavelet Transform (CWT) is defined as the sum over all time of the signal multiplied by scaled and shifted versions of the wavelet function ψ . Let R be the set of real numbers. The CWT of a signal x(t) is defined as ∞
CWT (a, b) =
∫ x(t )ψ
* a ,b
(t )dt ,
(5.1)
−∞
where
ψ a ,b (t ) = a
ψ ((t − b) / a ) .
−1 / 2
(5.2)
CHAPTER 5: THE WAVELET TRANSFORM METHOD
71
ψ (t ) is the mother wavelet, the asterisk in (5.1) denotes a complex conjugate, and 0 ≠ a ∈ R, and b ∈ R respectively are the dilation (mathematicians use the term “dilation” to refer to both compression and expansion) and translation parameters respectively. The term | a | −1 / 2 in (5.2) is the normalisation value of ψ a ,b (t ) so that if ψ (t ) has a unit length, then its scaled version ψ a ,b (t ) also has a unit length.
5.2.2 The Discrete Wavelet Transform Instead of continuous scaling (“dilation”) and shifting (“translation”), the mother wavelet may be scaled and shifted discretely by choosing a = a 0m , b = na 0m b0 , t = kT in (5.1) and (5.2), where T = 1.0 and k, m, n ∈ Z , (Z is the set of positive integers). The Discrete Wavelet Transform (DWT) is thus given by DWT (m, n) = a 0− m / 2
(∑ x[k ]ψ
*
)
[(k − na 0m b0 ) / a 0m ] .
(5.3)
By careful selection of a0 and b0 , the family of scaled and shifted mother wavelets constitutes an orthonormal basis. An orthonormal basis comprises a set of vectors S such that u.v = 0 (dot product) for each distinct pair of u, v ∈ S . We can choose a 0 = 2 and b0 = 1
to constitute the orthonormal basis to have the WT referred to as a dyadic-orthonormal WT. The implications of the dyadic-orthonormal WT is that because of the orthonormal properties there will be no information redundancy in the decomposed signals. Also, with this choice of a0 and b0 , there exists a novel algorithm, known as the multiresolution signal decomposition (Mallat, 1989) technique, to decompose a signal into scales with different time and frequency resolution.
72
CHAPTER 5: THE WAVELET TRANSFORM METHOD
5.2.3 Multiresolution Signal Decomposition and Quadrature Mirror Filter The Multiresolution Signal Decomposition (MSD) (Mallat, 1989) technique decomposes a given signal into its detailed and smoothed versions. Let x[n] be a discrete-time signal. Then the MSD technique decomposes the signal in the form of WT coefficients at scale 1 into c1[n] and d1[n], where c1[n] is the smoothed version of the original signal, and d1[n] is the detailed version of the original signal x[n]. These are defined as c1 [n] = ∑ h[k − 2n]x[k ] ,
(5.4)
d 1 [n] = ∑ g[k − 2n]x[k ] ,
(5.5)
k
k
where h[n] and g[n] are the associated filter coefficients that decompose x[n] into c1[n] and d1[n] respectively. Downsampling is done in the process of decomposition so that the resulting c1[n] and d1[n] are each n/2 point signals. Thus, for the original n point signal x[n], after the decomposition we have n point signal together with c1[n] and d1[n], not 2n point. The next higher scale decomposition will be based on c1[n]. Thus, the decomposition process can be iterated, with successive approximations being decomposed in turn, so that the original signal is broken down into many lower resolution components. This is called the wavelet decomposition tree (Mallat, 1998). The MSD technique can be realised with the cascaded Quadrature Mirror Filter (QMF) (Mallat, 1998) banks. A QMF pair consists of two finite impulse response filters, one being a lowpass filter (LPF) and the other a highpass filter (HPF). The QMF pair divides the input signal into low-frequency and high-frequency components at the dividing point halfway between zero Hz and half the data sampling frequency. The output of the lowpass
CHAPTER 5: THE WAVELET TRANSFORM METHOD
73
filter is the smoothed version of the input signal and is used as the next QMF pair’s input. The output of the highpass filter is the detailed version of the original signal. Thus cascaded QMF pairs realise the MSD technique. Detail description about QMF can be found in Strang and Nguyen (1996). Figure 5.1 shows the MSD technique and QMF pair.
FIGURE 5.1: Multiresolution Signal Decomposition realised by Quadrature Mirror Filter banks.
5.3
Signal Decomposition
In this section, we describe the application of the multiresolution signal decomposition technique and quadrature mirror filter banks to decompose the fault signals from the DFRs into localised and detailed representation in the form of wavelet coefficients. Daubechies 1 and 4 wavelets are used as mother wavelets, that is, the filters h[n] and g[n] as in (5.4) and (5.5) are chosen with one and four coefficients respectively and calculated as proposed by Daubechies (1992). Daubechies 1 wavelet can also be referred as the Haar wavelet first introduced by Alfred Haar (Haar, 1910). Daubechies 1 and 4 wavelets were selected in stead of other choices of the mother wavelets, e.g., Meyer wavelet (Meyer, 1993), Coiflets, Gaussian wavelet, Mexican hat wavelet,
74
CHAPTER 5: THE WAVELET TRANSFORM METHOD
Morlet wavelet (Daubechies, 1992), because they are compactly supported wavelets (Daubechies, 1992) with an extremal phase and the highest number of vanishing moments for a given support width (Daubechies, 1992). Furthermore, the associated scaling filters are minimum-Phase filters (Daubechies, 1992). Hence from the point of view of fast implementation and varying patterns of the fault signals, Daubechies wavelets appear to be the optimal choice for the mother wavelet for this specific application. Daubechies 1 (Haar) wavelet can be described in the following mathematical terms: The scaling function φ ( x) is defined as
φ ( x) = 1, if x ∈ [0,1] , φ ( x) = 0, if x ∉ [0,1].
(5.6)
The wavelet function ψ (x) for this scaling function is defined as
ψ ( x) = 1, if x ∈ [0,0.5] , ψ ( x) = −1, if x ∈ [0.5,1] , ψ ( x) = 0, if x ∉ [0,1].
(5.7)
Figure 5.2 shows the scaling (a) and wavelet (b) function for the Daubechies 1 (Haar) mother wavelet.
FIGURE 5.2: Daubechies 1 (Haar) wavelet: (a) Scaling function and (b) Wavelet function.
CHAPTER 5: THE WAVELET TRANSFORM METHOD
75
For the Daubechies 4 wavelet, the scaling function φ ( x) has the form
φ ( x) = c0φ (2 x) + c1φ (2 x − 1) + c 2φ (2 x − 2) + c3φ (2 x − 3) ,
(5.8)
c0 = (1 + 3 ) / 4 ,
(5.9)
c1 = (3 + 3 ) / 4 ,
(5.10)
c2 = (3 − 3 ) / 4 ,
(5.11)
c4 = (1 − 3 ) / 4 .
(5.12)
where
It is generally not possible to solve directly for φ (x) ; the approach is to solve for φ (x ) iteratively until φ j (x) is very nearly equal to φ j −1 ( x) , where
φ j ( x) = c0φ j −1 (2 x) + c1φ j −1 (2 x − 1) + c2φ j −1 (2 x − 2) + c3φ j −1 (2 x − 3) .
(5.13)
The Daubechies 4 wavelet function ψ (x) for the four-coefficient scaling function is given by
ψ ( x) = −c3φ (2 x) + c 2φ (2 x − 1) − c1φ (2 x − 2) + c0φ (2 x − 3) .
(5.14)
The Daubechies 4 scaling and wavelet function are shown in Figure 5.3. After transforming the original fault signal using the mother wavelets described above, we obtain the smoothed and detailed versions, namely, c1[n] and d1[n] [see (5.4), (5.5)]. Signal d1[n] can be regarded as the difference between the original signal x[n] and c1[n], and
76
CHAPTER 5: THE WAVELET TRANSFORM METHOD
called the wavelet transform coefficient at scale one. We will use d1[n] for threshold checking to estimate the abrupt change time-instants described in the following section.
(a)
(b)
FIGURE 5.3: Daubechies 4 wavelet: (a) Scaling function and (b) Wavelet function.
5.4
Application of the Threshold Method
The threshold method on the wavelet transform coefficients of the original fault signal is used to detect the jumps and sharp cusps (Wang, 1998) in order to estimate the timeinstants of the abrupt changes. Mathematically, we say that signal x(t) has a sharp α -cusp at t if for 0 ≤ α < 1 , | x(t + ∆t ) − x(t ) | ≥ K | ∆t |α ,
(5.15)
as ∆t → 0 for some constant K ≥ 0 . We can consider it to be a jump if α = 0 . In practice, we can consider a cusp to be an abrupt change of the level of the trend over a small time period. As discussed in the previous section, after transforming the original fault signal using the wavelet transform, we will search progressively across the finer wavelet scales for the largest wavelet coefficients on that scale (Craigmile & Percival, 2000). Since wavelet
CHAPTER 5: THE WAVELET TRANSFORM METHOD
77
coefficients are the changes of the averages, a coefficient of large magnitude implies a large change in the original signal. Large wavelet coefficients that are co-located in time across different scales provide estimates of the cusp points (Craigmile & Percival, 2000) hence time-instants of the abrupt changes. The change time-instants can be estimated by the instants when the wavelet coefficients exceed the first order approximation to the universal threshold of Donoho and Johnstone (Donoho & Johnstone, 1994). The universal threshold T is given (as in 4.23) by T = σ 2 log e n ,
(5.16)
where σ can be the median absolute deviation of the wavelet coefficients, or standard deviation, and n is the number of samples of the wavelet coefficients. Median absolute deviation is a good choice because median is hardly influenced by a small fraction of extreme values (Wang, 1998). However, experimenting with the real disturbance signals, standard deviation is finally chosen to cover a wider range of disturbance signals than that covered by the median absolute deviation. After determining the time-instants when the wavelet coefficients of the fault signal exceed the threshold, we mark them using unit impulses, indicating the abrupt change timeinstants.
5.5
The Application Results
In this section, we present the practical application results of the power system fault analysis method developed according to the above-mentioned signal decomposition and representation using the WT and then applying the threshold method to the detailed version of
78
CHAPTER 5: THE WAVELET TRANSFORM METHOD
the fault signal, followed by the heuristic smoothing filtering operation (Ukil & Zivanovic, 2005b). MATLAB® with Wavelet toolbox (Mathworks Inc., 2002b) has been used for implementing the application. The whole procedure detects the change time-instants thus segments the fault signal. We are interested to indicate the change time-instants as unit impulses. After normalising the original fault signal using its mean value, it is transformed into the smoothed and detailed version using the WT, whereafter the threshold method is applied on the detailed version to determine the change time-instants (Ukil & Zivanovic, 2005b). Smoothing filter operations (Ukil & Zivanovic, 2005b) are then applied on this segmented model to perform the following smoothing operations sequentially: •
It removes confusing multiple close-spikes and combines them into single unit impulse.
•
It removes any unwanted glitches which can otherwise result in false positives for the abrupt changes.
•
The segments in the power system fault analysis signals are during the pre-fault condition and following events such as fault initiation, circuit-breaker opening and reclosing. These events are predefined and so are the number of segments. Thus any larger numbers of segmentation possibly indicates transients, power swings and the like. The number of segment(s) is also estimated and checked.
•
Based on the modelling of the segments, an analysis is conducted to estimate the event-critical change instants and reject the others.
CHAPTER 5: THE WAVELET TRANSFORM METHOD
79
Figure 5.4 shows the result for the fault signal, sampled at a sampling frequency of 2.5 kHz (Stokes-Waller, 1998), obtained from the ESKOM DFRs during a phase-to-ground fault. The WT uses the Daubechies 1 (Daubechies, 1992) mother wavelet. In Figure 5.4, the original DFR recording for the current during the phase-to-ground fault in the RED-Phase is shown in the top section, wavelet coefficients for this fault signal (in blue) and the universal threshold (in black, dashed) are shown in the middle section. The change time-instants computed using the threshold checking (middle section) followed by smoothing filtering are shown in the bottom section. It should be noted that only the highpass filter output of the QMF pair is shown, so the wavelet coefficients in the middle section indicate half of the total samples of the original signal. The time-instants of the
FIGURE 5.4: Segmentation of the RED-Phase Current Signal using the Wavelet Transform Method.
80
CHAPTER 5: THE WAVELET TRANSFORM METHOD
changes in the signal characteristics, in the lower plot in Figure 5.4, indicate the different signal segments owing to different events during the fault, for example, segment A indicates the pre-fault section and the fault inception, segment B the fault, segment C the opening of the circuit-breaker, segment D auto-reclosing of the circuit-breaker and system restore. Figure 5.5 shows another result for the RED-Phase voltage recording, sampled at a sampling frequency of 2.5 kHz (Stokes-Waller, 1998), from the ESKOM DFRs during a phase-to-ground fault. The WT uses Daubechies 4 (Daubechies, 1992) mother wavelet.
FIGURE 5.5: Segmentation of the RED-Phase Voltage Signal using the Wavelet Transform Method.
CHAPTER 5: THE WAVELET TRANSFORM METHOD
81
In Figure 5.5, the original DFR recording for the voltage during the fault in the RED-Phase is shown in the top section, wavelet coefficients for this fault signal (in blue) and the universal threshold (in black, dashed) are shown in the middle section and the change timeinstants computed using the threshold checking (middle section) followed by smoothing filtering is shown in the bottom section. It should be noted that only the highpass filter output of the QMF pair is shown, hence the wavelet coefficients in the middle section indicate half of the total samples of the original signal. The time-instants of the changes in the signal characteristics in the lower plot in Figure 5.5 indicate the different signal segments caused by different events during the fault. For example, segment A indicates the pre-fault section and the fault inception, segment B the fault, segment C the opening of the circuitbreaker, segment D the auto-reclosing of the circuit-breaker and system restore.
5.6
Comments on Application Results
Following the discussion of the applied algorithms and the application results, the following comments apply (Ukil & Zivanovic, 2005b). •
Since the intended application is not meant for real-time analysis, computation time is not a critical factor. However, the proposed algorithm for the abrupt change detection and signal segmentation took an average computation time of 0.431 seconds. An Intel® Celeron® 1.9 GHz, 256 MB RAM computer was used for all the application tests using MATLAB® (Mathworks Inc., 2002b). It should be noted that the complete automatic disturbance recognition and analysis tasks have to be performed within five minutes of acquiring the fault signals, abrupt change detection and segmentation being the first step.
82
CHAPTER 5: THE WAVELET TRANSFORM METHOD
•
Two-hundred disturbance records have been tested with the proposed algorithm and 99% accuracy in correct segmentation was achieved.
•
The proposed algorithm using the wavelet transform is considerably faster and more robust compared with the traditional peak value detection and superimposed current quantities algorithms (Stokes-Waller, 1998). Also, this algorithm which is based on the wavelet transform and threshold method facilitates further signal processing and analysis in the subsequent stages of automatic disturbance recognition and analysis, focusing on the different segments and helping to quickly determine parameters such as the duration of the fault directly from the abrupt change detection-based segmentation itself. This cannot be done using the traditional peak value detection and superimposed current quantities algorithms (Stokes-Waller, 1998).
•
Instead of the Fourier transform, the wavelet transform is particularly suitable for the power system disturbance and fault signals which may not be periodic and may contain both sinusoidal and impulse components.
•
Wavelet coefficients are greatly adaptive to the fault signal pattern variations.
•
The wavelet transform provides a local representation (both in time and frequency) of a given signal, thus achieving the necessary time-frequency resolution for the power system fault analysis. This is not possible with the traditional Fourier transform which provides a global representation of a signal.
CHAPTER 5: THE WAVELET TRANSFORM METHOD
5.7
83
Summary
We have presented in this chapter the wavelet transform used for detecting the abrupt changes in the disturbance signals. Power system disturbance and fault signals may not be periodic and may contain both sinusoidal and impulse components. Hence we propose the use of the wavelet transform, particularly the dyadic-orthonormal the wavelet transform to decompose the original fault signal into the smoothed and detailed version in terms of the wavelet coefficients using the multiresolution signal decomposition technique. Then we make a progressive search on that wavelet scale for the largest wavelet coefficients. The change time-instants can be estimated by the time-instants when the wavelet coefficients exceed the first order of approximation to the universal threshold of Donoho and Johnstone. This is followed by the smoothing operation. We have been mainly interested in estimating the change time-instants and obtained good results from the MATLAB® implementation. Hence the use of the dyadic-orthonormal the wavelet transform to transform the fault signals into the smoothed and detailed version, followed by threshold checking is extremely effective in detecting the abrupt changes in the signals originating from power system faults in order to segment them into the event-specific sections.
84
CHAPTER 6: ADAPTIVE ABRUPT CHANGE DETECTION
CHAPTER 6 ADAPTIVE ABRUPT CHANGE DETECTION
Abrupt change detection based on recursive identification and the wavelet transform technique as discussed in the previous chapters, is highly effective in detecting the abrupt changes and hence segmenting the signals recorded during disturbances in the electrical power network of South Africa for disturbance analysis. The proposed recursive identification and wavelet method estimates exactly the time-instants of the changes in the signal model parameters during the pre-fault condition, after the initiation of fault, after the circuit-breaker opening and auto-reclosure etc. These are effective for about 60% of the disturbance signals, having distinct abrupt changes in the signal model parameters. However, about 40% of the disturbance signals do not show distinct abrupt changes in the signal parameters. In those cases, we have to apply adaptive abrupt change detection. This is the focus of this chapter. For the adaptive abrupt change detection, we propose two methods, (1) the Adaptive Whitening Filter (Ukil & Zivanovic, 2005e) (2) the Adjusted Haar Wavelet (Ukil & Zivanovic, 2005g).
6.1
The Adaptive Whitening Filter
The Adaptive whitening filter based on the adjusted Fourier filter (Wiot, 2004) is used to pre-filter the original fault signal. The wavelet transform method is then used, as discussed in the chapter 5, to transform the pre-filtered fault signal into the finer wavelet scales, fol-
CHAPTER 6: ADAPTIVE ABRUPT CHANGE DETECTION
85
lowed by a progressive search for the largest wavelet coefficients on that scale (Ukil & Zivanovic, 2005b) to perform the abrupt change detection, as depicted in Figure 6.1.
FIGURE 6.1: Architecture of Abrupt Change Detection based on the Adaptive Whitening Filter.
The adaptive whitening filter along with the associated optimal Linear Prediction Error (LPE) filter, and the adjusted Fourier filter have been discussed in details in chapter 3.4.2 and 3.4.3, so we skip further discussion on them. Derivation of the complete recursive equations (see equations 3.2 to 3.20) of the adaptive whitening filter following the LMS method can be found in the discussions by Wiot (2004) and Philippot (1996). This filter has the advantage of perfectly extracting the main frequency component of the signal, strongly attenuating its harmonics and the dc component (Ukil & Zivanovic, 2005e). Characteristics of the adaptive whitening filter in terms of its frequency response
86
CHAPTER 6: ADAPTIVE ABRUPT CHANGE DETECTION
are shown in Figure 6.2.
FIGURE 6.2: Frequency Response Plots of the Adaptive Whitening Filter.
6.2
Application of the Adaptive Whitening Filter
The fault signals from the DFRs are first filtered with the adaptive whitening filter. A sampling frequency of 2.5 kHz, the same as that of the DFRs (Stokes-Waller & Keller, 1998), is used for the adaptive whitening filter with the primary focus on the 50 Hz component. Hence we set the network pulsation frequency at 51 Hz, at close proximity to 50 Hz. This choice of the network pulsation frequency is optimal for this application, following several tests with different fault signals (Ukil & Zivanovic, 2005e). The pre-filtering operation extracts the main frequency component to be used for the signal decomposition, strongly attenuating its harmonics and the dc component. This operation
CHAPTER 6: ADAPTIVE ABRUPT CHANGE DETECTION
87
increases the precision and sensitivity of the next operations, namely signal decomposition and abrupt change detection (Ukil & Zivanovic, 2005e). This pre-filtering operation also helps to distinguish a fault from a transient recovery, a short-term swing and the like, which are otherwise extremely difficult to estimate using only the wavelet transform (Ukil & Zivanovic, 2005e). This is followed by signal decomposition using the wavelet transform, threshold checking on the decomposed signal and smoothing filtering (Ukil & Zivanovic, 2005e) to perform the abrupt change detection. Figure 6.3 shows a result for the fault signal obtained from the DFRs in the electrical power network of ESKOM, South Africa, during a phase-to-Phase fault. The fault signal is sampled at a sampling frequency of 2.5 kHz.
FIGURE 6.3: Segmentation of RED-Phase Voltage Signal using the Adaptive Whitening Filter and the Wavelet Transform Method.
88
CHAPTER 6: ADAPTIVE ABRUPT CHANGE DETECTION
In Figure 6.3, the original DFR recording for the voltage in the RED-Phase during the fault involving the RED- and BLUE-Phases is shown in plot (i), pre-filtering of the fault signal using the adaptive whitening filter in plot (ii), wavelet coefficients for the filtered fault signal (in continuous line) and the universal threshold (in dashed line) in plot (iii), and the abrupt change time-instants computed using the threshold checking followed by smoothing filtering in plot (iv). Note that because only the highpass filter output of the QMF pair is shown for the wavelet coefficients, so the wavelet coefficients in plot (iii) indicate half of the total number of samples of the original signal (Ukil & Zivanovic, 2005e). Figure 6.3 shows an example in which the fault signal follows a resistive decay, which does not produce the sharp transitions among the different segments caused by the different events during the fault. Since we cannot determine the different sharp segments in these cases, we rather focus on correctly estimating the time-instant of the fault-inception for further signal processing and analysis based on it (Ukil & Zivanovic, 2005e). In Figure 6.3, this is shown in plot (iv), where the impulse indicator shows the fault-inception timeinstant, thus segmenting the fault signal into two segments, pre-fault and post-fault segment indicated by segments A and B respectively. The application of the adaptive whitening filter for pre-filtering the fault signal is of particular importance for these kinds of signals. As shown in plot (ii) of Figure 6.3, the filtered fault signal clearly suppresses the resistive decay part and highlights the fault-inception time-instant, which improves the accuracy of detection of the fault-inception time-instant using the subsequent wavelet decomposition and threshold operations (Ukil & Zivanovic, 2005e).
CHAPTER 6: ADAPTIVE ABRUPT CHANGE DETECTION
6.3
89
The Adjusted Haar Wavelet
The adaptive whitening filter technique does not work well for all the special disturbance signals without distinct abrupt changes in the signal parameters. It sometimes also degrades the accuracy of the wavelet transform method for normal disturbance signals that have distinct abrupt changes. Also, for these special disturbance signals not showing distinct abrupt changes in the signal parameters, standard mother wavelets such as, Haar (Daubechies 1), Daubechies 4 (Daubechies, 1992) fail to achieve correct event-specific signal segmentations. Hence Ukil and Zivanovic (2005g) propose a new adjustment technique to the standard Haar wavelet by introducing 2n adjusting zeroes in the Haar wavelet scaling filter. This technique is fairly effective in segmenting these fault signals into preand post-fault segments, and is an improvement on the standard mother wavelets for this application.
6.3.1 Overview of the Haar Wavelet Haar wavelet was first introduced by Alfred Haar in 1910 (Haar, 1910). It has the following mathematical description (discussed in chapter 5.3, mentioned here briefly for the sake of easy and continuous following of the current discussion). The scaling function φ ( x) is defined as
φ ( x) = 1, if x ∈ [0,1] , φ ( x) = 0, if x ∉ [0,1]. The wavelet function ψ (x) for this scaling function is defined as
(6.1)
90
CHAPTER 6: ADAPTIVE ABRUPT CHANGE DETECTION
ψ ( x) = 1, if x ∈ [0,0.5] , ψ ( x) = −1, if x ∈ [0.5,1] , ψ ( x) = 0, if x ∉ [0,1].
(6.2)
Refer to Figure 5.2 (chapter 5.3) for the graphical displays of the scaling function and wavelet function for the Haar wavelet.
6.3.2 Adjustment to the Haar Wavelet In this section, we discuss the proposed adjusted Haar wavelet in terms of the key properties of the wavelets. In general, the FIR scaling filter for the Haar wavelet is h = 0.5 [1 1] , where 0.5 is the normalisation factor. As an adjustment to and hence an improvement of the characteristics of the Haar wavelet, we propose to introduce 2n zeroes (n is a positive integer) in the Haar wavelet scaling filter, keeping the first and last coefficients 1 as proposed by Ukil and Zivanovic (2005g). Following the orthogonality property for the scaling filter, the filter length has to be even (Strang &Nguyen, 1996). It is thus necessary to introduce 2n adjusting zeroes, n being the adjustment parameter (Ukil & Zivanovic, 2005g). The Haar wavelet corresponds to n = 0 . The introduced additional zeroes in the filter kernel have zero coefficients. The scaling filter kernel for the adjustment parameter is shown below. h = 0.5[1 1] h = 0.5[1 0 0 1]
for n = 0 for n = 1
h = 0.5[1 0 0 0 0 1] for n = 2
. . .
(6.3)
CHAPTER 6: ADAPTIVE ABRUPT CHANGE DETECTION
91
It can be shown that the adjusting zeroes improve the Haar wavelet characteristics, especially for application in this research. Also, it can be shown mathematically that the introduction of the adjusting zeroes does not violate the key wavelet properties like compact support, orthogonality and perfect reconstruction. We prove these results in section 6.3.3.
6.3.3 Compact Support LEMMA 6.1: The adjusted Haar wavelet with 2n of adjusting zeroes has compact support. Proof: The scaling filter of the Haar wavelet comes from an FIR filter, with finite length.
The original Haar wavelet has the closed and bounded interval [0,1] as its support and so the support is compact. Compact support corresponds to FIR. The adjusted Haar wavelet scaling filter with the additional adjusting zeroes also forms an FIR filter kernel. Also, because the adjusting zeroes have zero coefficients, it also has compact support: [0,1] (Ukil & Zivanovic, 2005g). ■
6.3.4 Orthogonality Real vectors are orthogonal (perpendicular) when x.y = 0 (dot or inner product). Real functions are said to be orthogonal when
∫ X (ω )Y (ω )dω = 0 . If the vectors or the func-
tions are complex valued, it is necessary to consider complex conjugates of one vector or one function. The discrete analogue of an orthonormal transform is a square matrix with orthonormal columns (Strang &Nguyen, 1996). This is an “orthogonal” matrix if real, a “unitary” matrix if complex. The orthogonal filter bank comes from an orthogonal matrix
92
CHAPTER 6: ADAPTIVE ABRUPT CHANGE DETECTION
AtT At = I
and At AtT = I .
(6.4)
LEMMA 6.2: The adjusted Haar scaling filter with 2n adjusting zeroes is a symmetric,
orthogonal FIR filter. Proof: An FIR filter H ( z ) is symmetric when z N H ( z ) = H ( z −1 ) . For orthogonality, N is
odd and the filter length must be even (Strang &Nguyen, 1996). The adjusted Haar scaling filter with two nonzero coefficients and 2n adjusting zeroes with zero coefficients, form a filter kernel with N = 2n + 1 . If this is a symmetric filter, it has a length of 2n + 2 and has the form (h(0), h(1), h(2),..., h(n), h(n),..., h(2), h(1), h(0) ) , and by convention, h(0) is the first nonzero coefficient. This vector must be orthogonal to all its double shifts (Strang & Nguyen, 1996). The inner product with its shift by 1 must be 2h(0)h(1) = 0 ; so h(1) = 0 . Then the inner product with its shift by 2 gives 2h(0)h(2) = 0 ; so h(2) = 0 . Continuing in this vein, the inner product with its shift by n gives 2h(0)h(n) = 0 ; so h(n) = 0 . Hence the only nonzero coefficient for the symmetric, orthogonal filter is the h(0) at both ends of the filter. Therefore the adjusted Haar scaling filter with two nonzero coefficients (equal to 1) at both ends and 2n adjusting zeroes with zero coefficients embedded in between form a symmetric, orthogonal FIR filter kernel (Ukil & Zivanovic, 2005g). ■
6.3.5 Perfect Reconstruction
The perfect reconstruction condition for a lowpass filter P0 ( z ) is P0 ( z ) − P0 (− z ) = 2 z −l .
(6.5)
93
CHAPTER 6: ADAPTIVE ABRUPT CHANGE DETECTION
The equation (6.5) can be simplified as discussed by Strang and Nguyen (1996), so that the perfect reconstruction condition states that the filter P( z ) must be a “halfband filter” (Strang & Nguyen, 1996), so that
P( z ) + P(− z ) = 2 .
(6.6)
LEMMA 6.3: The introduction of the 2n adjusting zeroes to the Haar wavelet scaling filter
satisfies the perfect reconstruction condition. Proof: The original Haar wavelet scaling filter has the form h = [1 1] , that is,
(
H 0 (ω ) = 1 + e − jω
gives
the
(
)
(
)
and H 0 ( z ) = 1 + z −1 . The introduction of the 2n adjusting zeroes
scaling
filter
as
h = [1 0 0 ... 0 0 1] ,
i.e.,
(
H 0 (ω ) = 1 + e − j ( 2 n +1)ω
)
and
)
H 0 ( z ) = 1 + z − ( 2 n +1) , where n is a positive integer. The original Haar wavelet corresponds
to n = 0 . Thus, as per the perfect reconstruction condition shown in (6.6), for the adjusted Haar wavelet scaling filter we obtain,
(
) (
)
H 0 ( z ) + H 0 (− z ) = 1 + z − ( 2 n +1) + 1 − z − ( 2 n +1) = 2 .
(6.7)
This completes the proof (Ukil & Zivanovic, 2005g). ■
6.3.6 The Adjusted Scaling Function Following LEMMA 6.1, 6.2 and 6.3, we have established that the adjusted Haar wavelet scaling filter, with 2n adjusting zeroes, satisfies the key wavelet properties such as compact support, orthogonality and perfect reconstruction (Ukil & Zivanovic, 2005g).
94
CHAPTER 6: ADAPTIVE ABRUPT CHANGE DETECTION
In the frequency domain, the adjusted Haar wavelet scaling filter with 2n adjusting zeroes is given as
H 0 (ω ) =
(
)
1 1 + e − j ( 2 n +1)ω . 2
(6.8)
Obviously, the function H (ω ) in (6.8) is a lowpass filter with H 0 (0) = 1 and H 0 (π ) = 0 . The formula for the Fourier transform of the scaling function based on the lowpass filter is given below as discussed by Qian (2002). ∞ ω Φ (ω ) = H 0 (ω2 )Φ (ω2 ) = H 0 (ω2 )H 0 (ω4 )Φ (ω4 ) = ∏ H 0 k Φ(0) . 2 k =1
(6.9)
Using (6.8) and (6.9), the Fourier transform of the adjusted Haar wavelet scaling function is ω − j ( 2 n +1) k ω ∞ 1 2 Φ(ω ) = ∏ H 0 k = ∏ 1 + e 2 k =1 2 k =1 ∞
.
(6.10)
So,
∞
Φ(ω ) = ∏ e
− j ( 2 n +1)
ω 2 k +1
k =1
∞
= ∏e k =1
− j ( 2 n +1)
ω 2 k +1
ω
ω
− j ( 2 n +1) k +1 1 j ( 2 n +1) 2 k +1 2 e +e 2
(2n + 1)ω cos k +1 2
∞ 1 ∞ (2n + 1)ω = exp− j (2n + 1)ω ∑ k +1 ∏ cos k +1 2 k =1 2 k =1
(6.11)
CHAPTER 6: ADAPTIVE ABRUPT CHANGE DETECTION
=e
− j ( 2 n +1)
ω ∞ 2
95
(2n + 1)ω . 2 k +1
∏ cos k =1
As per formula 1.439, on page 38 of Gradshteyn and Ryshik (1965), ∞
ω sin (ω 2 ) , = k +1 ω 2
∏ cos 2 k =1
(6.12)
we obtain the adjusted scaling function as,
Φ(ω ) = e
− j ( 2 n +1)
ω 2
sin ((2n + 1) ω 2) . (2n + 1)ω 2
(6.13)
Figures 6.4, 6.5, and 6.6 show the pole-zero plots of the adjusted Haar wavelet scaling filters for n = 0,1,2 respectively. It should to be noted that the original Haar wavelet scaling filter corresponds to n = 0 , and additional complex conjugate pairs of zeroes for each n > 0 are introduced (Ukil & Zivanovic, 2005g).
FIGURE 6.4: Pole-Zero plot of the Adjusted Haar Wavelet Scaling Filter, for n=0, which corresponds to the Original Haar Wavelet.
96
CHAPTER 6: ADAPTIVE ABRUPT CHANGE DETECTION
FIGURE 6.5: Pole-Zero plot of the Adjusted Haar Wavelet Scaling Filter with Adjusting Zeroes, for n=1, which introduces one pair of Complex Conjugate Zeroes.
FIGURE 6.6: Pole-Zero plot of the Adjusted Haar Wavelet Scaling Filter with Adjusting Zeroes, for n=2, which introduces two pairs of Complex Conjugate Zeroes.
6.3.7 The Adjusted Wavelet Function To compute orthogonal mother wavelets from the lowpass filter H 0 (ω ) , another function H 1 (ω ) is required such that
CHAPTER 6: ADAPTIVE ABRUPT CHANGE DETECTION
H 0 (ω ) H 1* (ω ) + H 0 (ω + π ) H 1* (ω + π ) = 0 .
97
(6.14)
This is condition for the quadrature mirror filter as shown by Strang and Nguyen (1996) and Qian (2002). One solution to (6.14) is H 1 (ω ) = −e − jω H 0* (ω + π ) .
(6.15)
Substituting H 0 (0) = 1 and H 0 (π ) = 0 into (6.15) yields H 1 (0) = 0 and H 1 (π ) = 1 , respectively. This means that H 1 (ω ) in (6.15) is a highpass filter. Thus, for the adjusted Haar wavelet scaling (lowpass) filter H 0 (ω ) [see (6.8)], the highpass filter H 1 (ω ) is given by
H 1 (ω ) = −e − jω H 0* (ω + π ) = −e − jω =
(
(
1 1 − e j ( 2 n +1)ω 2
)
(6.16)
)
1 j 2 nω e − e − jω . 2
Obviously, H 0 (ω ) and H 1 (ω ) constitute quadrature mirror filters, specified by (6.14). We can compute the Fourier transform of the wavelet function as discussed by Strang and Nguyen (1996), by using ω ω Ψ (ω ) = H 1 Φ . 2 2
(6.17)
Now, we establish the main result, which is as follows.
THEOREM 6.4: The introduction of the 2n adjusting zeroes to the Haar wavelet scaling filter improves the frequency characteristics of the adjusted wavelet function by an order of 2n+1.
98
CHAPTER 6: ADAPTIVE ABRUPT CHANGE DETECTION
Proof: For any n ≥ 0 , according to (6.17), the Fourier transform of the adjusted wavelet
function is using (6.13) and (6.16),
(
)
− j ( 2 n +1) 1 sin ((2n + 1) ω 4) − jω 4 Ψn (ω ) = e jnω − e 2 e 2 (2n + 1)ω 4
=
ω
(6.18)
)
(
1 j ( 2 n−1) ω4 − j ( 2 n +3) ω4 sin (( 2n + 1) ω 4 ) e −e , 2 (2n + 1)ω 4
with magnitude
2 sin ((2n + 1) ω 4) Ψn (ω ) = 2 {1 − cos((4n + 2) ω 4)} (2n + 1)ω =
=
{2 sin ((2n + 1) ω 4)}
2
{sin ((2n + 1) ω 4)}2 (2n + 1) ω 4
<
2 sin ((2n + 1) ω 4) (2n + 1)ω
2
2
4 . (2n + 1)ω
(6.19)
The factor 2n+1 in the denominator of (6.19) improves the frequency characteristics of the adjusted Haar wavelet function, by decreasing the ripples (as n > 0 ). This completes the proof (Ukil & Zivanovic, 2005g). ■ The following figures (6.7 to 6.9) illustrate the proof above. The ripple is measured by the difference of the heights of the main (first) lobe and the secondary (second onwards) lobes in the magnitude plot. It is to be noted in Figure 6.7 to 6.9, that the aforesaid differences of heights increase gradually for a specific frequency band. This indicates improvement of the ripple as the ripples die away gradually more quickly from Figure 6.7 to 6.9.
CHAPTER 6: ADAPTIVE ABRUPT CHANGE DETECTION
99
FIGURE 6.7: The Fourier Spectrum of the Adjusted Haar Wavelet, for n=0, which corresponds to the Original Haar Wavelet.
FIGURE 6.8: The Fourier Spectrum of the Adjusted Haar Wavelet, for n=1, which decreases the strong ripples of the Original Haar Wavelet.
100
CHAPTER 6: ADAPTIVE ABRUPT CHANGE DETECTION
FIGURE 6.9: The Fourier Spectrum of the Adjusted Haar Wavelet, for n=2, which further decreases the strong ripples of the Original Haar Wavelet.
6.4
Application of the Adjusted Haar Wavelet
After normalising the original fault signal using its mean value, it is transformed into the smoothed and detailed version using the wavelet transform, with the adjusted Haar wavelet as the mother wavelet. Adjustment parameter n=2 is applied, that is, four adjusting zeroes are included in the adjusting Haar wavelet scaling filter (Ukil & Zivanovic, 2005g). The threshold method (Ukil & Zivanovic, 2005b) is then applied to the detailed version to determine the abrupt change time-instants. This is followed by smoothing filter operations (Ukil & Zivanovic, 2005b) to indicate the change time-instants as unit impulses. MATLAB® with the Wavelet toolbox (Mathworks Inc., 2002b) was used to implement the application.
CHAPTER 6: ADAPTIVE ABRUPT CHANGE DETECTION
101
Figures 6.10 (a, b) to 6.13 (a, b) show the comparative results of the application of the original Haar wavelet (plot a) and the adjusted Haar wavelet (plot b) on the fault signals, sampled at a sampling frequency of 2.5 kHz (Stokes-Waller, 1998), obtained from the DFRS of ESKOM, South Africa, during various disturbances. For these special disturbance signals, the original Haar wavelet fails to achieve correct segmentation whereas the adjusted Haar wavelet correctly segments the fault signals into pre- and post-fault segments shown as A and B respectively in the plot (b) of Figure 6.10 to 6.13 (Ukil & Zivanovic, 2005g). Figures 6.10 and 6.11 show segmentation of the voltage waveform during phase-to-ground faults. The voltage waveform is shown in continuous lines and segmentations as vertical dashed unit impulses. Plot (a) shows the incorrect segmentation using the Haar wavelet as the mother wavelet, and plot (b) shows the correct segmentation using the adjusted Haar wavelet. In plot (b), the fault signal is segmented into pre- and post-fault segments, indicated as A and B respectively, based on the fault inception timing.
(a)
(b)
FIGURE 6.10: Segmentation of the Voltage Waveform during a Phase-to-Ground Fault. Plot (a) shows the Result using the Haar Wavelet, while Plot (b) shows the Result using the Adjusted Haar Wavelet.
102
CHAPTER 6: ADAPTIVE ABRUPT CHANGE DETECTION
(a)
(b)
FIGURE 6.11: Segmentation of the Voltage Waveform during a Phase-to-Ground Fault. Plot (a) shows the Result using the Haar Wavelet, while Plot (b) shows the Result using the Adjusted Haar Wavelet.
Figures 6.12 and 6.13 show segmentation of the voltage waveform during phase-to-phase faults. The voltage waveform is shown in continuous lines, and segmentations as vertical dashed unit impulses. Plot (a) shows the incorrect segmentation using the Haar wavelet as the mother wavelet, whereas plot (b) shows the correct segmentation using the adjusted Haar wavelet. In plot (b), the fault signal is segmented into pre- and post-fault segments, indicated as A and B respectively, based on the fault inception timing.
(a)
(b)
FIGURE 6.12: Segmentation of the Voltage Waveform during a Phase-to-Phase Fault. Plot (a) shows the Result using the Haar Wavelet, while Plot (b) shows the Result using the Adjusted Haar Wavelet.
CHAPTER 6: ADAPTIVE ABRUPT CHANGE DETECTION
(a)
103
(b)
FIGURE 6.13: Segmentation of the Voltage Waveform during a Phase-to-Phase Fault. Plot (a) shows the Result using the Haar Wavelet, while Plot (b) shows the Result using the Adjusted Haar Wavelet.
6.4.1 Comments on the Application Results The following comments apply to the application results of the adjusted Haar wavelet method (Ukil & Zivanovic, 2005g). •
The proposed algorithm using the adjusted Haar wavelet for the abrupt change detection and signal segmentation took an average computation time of 0.501 seconds, which is the same as that using the original Haar wavelet as discussed in Ukil and Zivanovic (2005b). Hence the adjustment algorithm has no impact on the speed of operation. An Intel® Celeron® 1.9 GHz, 256 MB RAM computer was used for all the application tests using MATLAB® (Mathworks Inc., 2002b).
•
One-hundred and fifty critical disturbance records, which were unsuccessfully segmented using the Haar wavelet, were tested with the proposed adjusted Haar wavelet. The accuracy of correct segmentation, in terms of pre- and post-fault segments, was 98%.
104
•
CHAPTER 6: ADAPTIVE ABRUPT CHANGE DETECTION
Since signals did not show distinct abrupt changes in the signal parameters, segmentation was done on the basis of the fault-inception time. Further processing of the signals concentrated on the post-fault segment. The proposed adjusted Haar wavelet is extremely accurate in estimating the fault-inception time for the disturbance signals without distinct abrupt changes in the signal parameters.
6.5
Summary
In this chapter we have discussed about the adaptive abrupt change detection for the disturbance signals not showing distinct changes in the signal model parameters. We have proposed two methods, namely the adaptive whitening filter and the adjusted Haar wavelet. The first method is an improvement technique based on the adaptive whitening filter. We first filter the fault signals using the adaptive whitening filter, focusing on the main frequency component to be used for signal decomposition, strongly attenuating its harmonics and the dc component. Then we propose the use of the wavelet transform, as discussed in Chapter 5, to decompose the filtered fault signal into a smoothed and detailed version in terms of wavelet coefficients, finally making a progressive search on that wavelet scale for the largest wavelet coefficients to accomplish the abrupt change detection. In the second method, we propose a general adjusting technique to the standard Haar wavelet. We propose to introduce 2n adjusting zeroes in the Haar wavelet scaling filter. It was also mathematically established that the adjusted Haar wavelet scaling filter, with 2n adjusting zeroes, and the resulting adjusted wavelet function satisfy the key wavelet properties such as compact support, orthogonality and perfect reconstruction. This generates a
CHAPTER 6: ADAPTIVE ABRUPT CHANGE DETECTION
105
series of new adjusted Haar wavelets which can be successfully used as the basis (mother wavelet) for the wavelet transform algorithm for accomplishing the adaptive abrupt change detection task.
106
CHAPTER 7: THE COMPLETE ALGORITHM
CHAPTER 7 THE COMPLETE ALGORITHM
This chapter describes the overall optimised algorithm for abrupt change detection-based segmentation. The aim is to take into consideration all the individual algorithms developed and described in the previous chapters and combine them in an optimised fashion to form the overall algorithm. The overall algorithm is the best optimised solution for segmentation of all kinds of disturbance signals. We also present the software implementation details.
7.1
Introduction
Figure 7.1 shows the top-down architecture of the complete segmentation algorithm.
FIGURE 7.1: Top-down Architecture of the Complete Segmentation Algorithm.
CHAPTER 7: THE COMPLETE ALGORITHM
107
According to Figure 7.1, the input into the segmentation algorithm is the disturbance signal from the digital fault recorders (DFRs), while the output is the segmented disturbance signal. The complete algorithm has been implemented in MATLAB® (MATHWORKS INC., 2002a), (MATHWORKS INC., 2002b). The sequential blocks are as follows: (1) the Disturbance Signal Read Module (2) the Signal Representation Algorithms (3) the Threshold Checking Algorithm (4) the Heuristic Smoothing Filtering (5) the Decision-Making Algorithm. We will explain them in details in the following sections.
7.2
The Disturbance Signal Read Module
This module reads the COMTRADE data (*.DAT) (IEEE STANDARD C37.111-1991, 1991) file. A detailed discussion of the COMTRADE format is provided in Appendix A. The read module reads the data columns from the data file for the analogue and binary values as needed. It also performs the normalisation of the read value using the mean of the read data. In the implementation step, there are two possibilities: 1. to read a single data column specified by the user 2. to read data column automatically in a sequential manner. Details of the implementation of the read module can be found in any of the MATLAB® scripts and functions on segmentation in Appendix B. The read module is generally found at the beginning of the scripts and functions.
108
7.3
CHAPTER 7: THE COMPLETE ALGORITHM
The Signal Representation Algorithms
After the signal values have been read and normalised, signal representation algorithms are applied to represent the signal effectively for the purpose of event-specific segmentation. The individual algorithms developed and described in the previous chapters serve as the signal representation algorithms. These are as follows: (1) the Recursive Identification method (2) the Wavelet Transform method (3) the Adaptive Whitening Filter method (4) the Adjusted Haar Wavelet method.
7.3.1 The Recursive Identification method Details of the recursive identification method can be found in Chapter 4. The recursive identification algorithm performs the segmentation assuming a parametric system model, and considering the disturbance signal as a quasi-stationary sequence. The recursive identification technique uses M parallel Kalman filters (Ukil & Zivanovic, 2005f). The software implementation of the recursive identification method can be found in Appendix B in the MATLAB® scripts and functions: ‘abrupt_recursive.m’, ‘func_ recursive.m’ and ‘test_recursive.m’.
7.3.2 Wavelet Transform method Details of the wavelet transform method are provided in Chapter 5. The wavelet transform is used to transform the original fault signal into finer wavelet scales, followed by a progressive search for the largest wavelet coefficients on that scale. Large wavelet coefficients
CHAPTER 7: THE COMPLETE ALGORITHM
109
that are co-located in time across different scales provide estimates of the changes in the signal parameter (Ukil & Zivanovic, 2005b). This method is faster than the recursive identification method. The Daubechies 4 (Daubechies, 1992) mother wavelet is mainly used. The software implementation of the wavelet transform method can be found in Appendix B in the MATLAB® scripts and functions: ‘abrupt_wavelet.m’, ‘func_wavelet.m’ and ‘test_wavelet.m’.
7.3.3 The Adaptive Whitening Filter method This method is required for special kinds of disturbance signals which do not show distinct abrupt changes in the signal parameters. Details of the adaptive whitening filter method are provided in Chapter 6 under the topic of adaptive abrupt change detection. The adaptive whitening filter based on the adjusted Fourier filter (Wiot, 2004) is used to pre-filter the original fault signal, followed by the wavelet transform method. The adaptive whitening filter method acts as an improvement step, extracting the main frequency component to be used for the signal decomposition, strongly attenuating its harmonics and the dc component (Ukil & Zivanovic, 2005e). The software implementation of the adaptive whitening filter and the associated segmentation technique is provided in Appendix B in the MATLAB® scripts and functions: ‘abrupt_adamo.m’, ‘func_adamo.m’ and ‘test_adamo.m’.
7.3.4 The Adjusted Haar Wavelet method This method is required for special kinds of disturbance signals which do not show distinct abrupt changes in the signal parameters. Details of the adjusted Haar wavelet method are provided in Chapter 6 under the topic of adaptive abrupt change detection. A new adjust-
110
CHAPTER 7: THE COMPLETE ALGORITHM
ment technique applied to the standard Haar wavelet is proposed by Ukil and Zivanovic (2005g), by introducing 2n adjusting zeroes in the Haar wavelet scaling filter. This general adjustment creates a new series of adjusted Haar mother wavelets. This technique is reasonably effective in segmenting those fault signals into pre- and post-fault segments, and is an improvement on the standard mother wavelets for this application. The software implementation of the general adjusted Haar wavelet is provided in Appendix B in the MATLAB® scripts: ‘abhicreate.m’ and ‘abhiwavelet.m’. Application of the new series of adjusted mother wavelets for the segmentation is provided in Appendix-B in the MATLAB® scripts and functions: ‘abrupt_mywavelet.m’, ‘func_mywavelet.m’ and ‘test_mywavelet.m’.
7.4
The Threshold Checking Algorithm
We use the threshold method on the signal representation to detect the jumps and sharp cusps (Wang, 1998) in order to estimate the time-instants of the abrupt changes. Details of the threshold checking algorithm are provided in Chapter 5 (section 5.4). Details of the implementation of the threshold checking algorithm can be found in any of the MATLAB® scripts and functions on segmentation in Appendix B. The threshold checking algorithm is mainly used with the wavelet transform method, as well as in the smoothing filter part of the recursive identification method.
7.5
Heuristic Smoothing Filtering
Heuristic smoothing filtering (Ukil & Zivanovic, 2005b) is applied to the segmented model to perform the following smoothing operations sequentially:
CHAPTER 7: THE COMPLETE ALGORITHM
111
(1) It removes confusing multiple close spikes and combines them into single unit impulse. (2) It removes any unwanted glitches which could otherwise result in false positives for the abrupt changes. (3) The segments in the power system fault analysis signals occur during the pre-fault condition and after events such as fault initiation, circuit-breaker opening and reclosing. These events are predefined, as are the number of segments. Thus any larger number of segmentation possibly indicates transients, power swings, and the like. An estimation of the number of segment(s) is also performed and checked. (4) Based on the modelling of the segments, analysis is done to estimate the eventcritical abrupt change instants and reject others. Of these, operation 1 is the main objective of the smoothing filter. This is implemented by taking the group or band of close signal points which pass a given threshold, and then compressing them into one effective unit impulse indicating the abrupt change. Each point in the band of signal is checked against its contiguous data points for consistency in the abrupt changes, and is either selected or rejected against its neighbours. The median is selected for a group of similar data points. The software implementation of this operation is provided in Appendix B in the MATLAB® function: ‘func_compressband.m’. Details of the implementation of the operation (2),(3) and (4) listed above can be found in any of the MATLAB® scripts and functions on segmentation given in Appendix B. These smoothing filtering operations can be generally found at the end of the scripts, functions, after operation (1), that is, after calling the function ‘func_compressband.m’.
112
7.6
CHAPTER 7: THE COMPLETE ALGORITHM
The Decision-Making Algorithm
The decision making algorithm optimises the selection of signal representation algorithms for the best possible result considering the specific disturbance signal to be processed. It takes into account the number of segments as the test parameter and decides whether or not the segmentation is correct. If the segmentation is incorrect and needs to be refined, the decision making algorithm chooses the proper signal representation algorithm and performs the segmentation again. This is repeated according to a specific optimisation algorithm as described below, until the final proper segmentation is achieved. The optimisation algorithm is depicted by the flowchart in Figure 7.2. The decision making algorithm only needs to handle signals with fewer than four segments, because a larger number of segmentations possibly indicates power swings, and is handled in the preceding heuristic smoothing filtering step (see Figure 7.1). The initial signal representation algorithm is the adjusted Haar wavelet with adjustment parameter = 2. The decision making algorithm subsequently uses the following: 1. the Recursive Identification method 2. the Wavelet Transform method with the Daubechies 4 mother wavelet 3. the Adaptive Whitening Filter method. The following considerations are made in the flowchart depicted in Figure 7.2. •
Output refers to the final segmented signal output.
•
The Check (algorithm) step comprises three procedures,
CHAPTER 7: THE COMPLETE ALGORITHM
113
(1) Signal Representation (2) Threshold Checking (3) Heuristic Smoothing Filtering. •
The signal representation algorithms in the Check (algorithm) are represented like: (1) the Recursive Identification method: ‘Recursive’ (2) the Wavelet Transform method, with the Daubechies 4 mother wavelet: ‘db4’ (3) the Adaptive Whitening Filter method: ‘ADAMO’ (4) the Adjusted Haar Wavelet method: ‘abhi’.
114
CHAPTER 7: THE COMPLETE ALGORITHM
FIGURE 7.2: The Decision-Making Algorithm Flowchart.
CHAPTER 7: THE COMPLETE ALGORITHM
115
The decision making algorithm and the flowchart shown in Figure 7.2 are explained by the pseudo code as follows.
Start with the initial segmented signal based on ‘abhi’ algorithm Is No. of Segment < 3? No, Output the segmented signal. Yes, Operation (Switch) based on No. of Segment Case: No. of Segment = 0 Check ‘db4’ Output the segmented signal. End Case: No. of Segment = 0
Case: No. of Segment = 1 Check ‘db4’ Operation (Switch) based on No. of Segment Case: No. of Segment > 2 Check ‘ADAMO’ Output the segmented signal. End Case Case: No. of Segment < 2 Output the segmented signal. End Case
116
CHAPTER 7: THE COMPLETE ALGORITHM
Case: No. of Segment = 2 Check ‘Recursive’ Is No. of Segment ≥ 2? No, Output the segmented signal. Yes, Check ‘db4’ Output the segmented signal. End Case End Operation End Case: No. of Segment = 1 Case: No. of Segment = 2 Check ‘db4’ Operation (Switch) based on No. of Segment Case: No. of Segment = 0 Output the segmented signal. End Case Case: No. of Segment = 1 Check ‘Recursive’ Is No. of Segment ≤ 2? Yes, Output the segmented signal. No,
CHAPTER 7: THE COMPLETE ALGORITHM
Check ‘abhi’ Output the segmented signal. End Case Case: No. of Segment = 2 Check ‘abhi’ Output the segmented signal. End Case Case: No. of Segment = 3 Check ‘ADAMO’ Output the segmented signal. End Case Case: No. of Segment > 3 Check ‘Recursive’ Is No. of Segment ≤ 3? Yes, Output the segmented signal. No, Check ‘abhi’ Output the segmented signal. End Case
End Operation
End Case: No. of Segment = 2
End Operation
117
118
CHAPTER 7: THE COMPLETE ALGORITHM
The final output of the decision-making algorithm is the properly segmented signal for further processing. The software implementation of the decision making algorithm is provided in Appendix B in the MATLAB® scripts and functions: ‘abrupt_segment.m’, ‘test_ segment.m’ and ‘func_segment.m’. The algorithm took an average computation time of 2.809 seconds, and was tested on 550 signals with 99% accuracy (Ukil & Zivanovic, 2006).
7.7
Summary
We have discussed in this chapter the complete segmentation algorithm. We have also discussed in detail the architectural structure of the software implementation of the overall algorithm. The complete segmentation algorithm has five major modules: the Disturbance Signal Read Module, the Signal Representation Algorithms, the Threshold Checking Algorithm, Heuristic Smoothing Filtering and the Decision Making Algorithm. After the disturbance signal is read, it is represented using the signal representation algorithms which comprise all the developed algorithms for the segmentations described in the previous chapters. The starting algorithm is the adjusted Haar wavelet method. Segmentation is then performed using the threshold-checking algorithm. The segmentation is refined using the heuristic smoothing filtering. The final decision-making algorithm checks whether the segmentation is correct for the specific disturbance signal. If it is incorrect and needs to be refined, decision-making algorithm iteratively chooses the proper signal representation algorithm and performs the segmentation until it is correct. Final output is the properly segmented disturbance signal. The complete algorithm takes into consideration all segmentation algorithms and is optimised for all kinds of disturbance signals based on testing on real signals.
CHAPTER 8: APPLICATIONS
119
CHAPTER 8 APPLICATIONS
The various abrupt change detection algorithms described in the previous chapters are extremely effective in detecting the abrupt changes in the signals recorded during disturbances in the electrical power network for disturbance analysis. The aim is to estimate exactly the time-instants of the changes in the signal model parameters during the pre-fault condition, after initiation of fault, after circuit-breaker opening and auto-reclosure etc. After the event-specific segmentations, synchronisation of the different digital fault recorder (DFR) recordings is necessary on the basis of the fault inception timings. The synchronised signals are segmented again. This synchronised segmentation is the first step towards automatic disturbance recognition, facilitating further complex feature vector analysis and pattern recognition. Moreover, the synchronised, segmented recordings can be directly used to monitor the relay performance and analyse certain kinds of disturbances. This chapter describes the various applications of the abrupt change detection in the power systems domain.
8.1
Abrupt Change Detection-based Segmentation
Figure 8.1 shows the abrupt change detection result using the wavelet transform and threshold method as discussed by Ukil and Zivanovic (2005b). In Figure 8.1, the original DFR recording from the ESKOM transmission network for the voltage during a phase-toground fault in the BLUE-Phase, sampled at a frequency of 2.5 kHz (Stokes-Waller, 1998), is shown in the top section; wavelet coefficients for this fault signal and the univer-
120
CHAPTER 8: APPLICATIONS
sal threshold (dashed) are shown in the middle section and the change time-instants as unit impulses computed using threshold checking (middle section); and the smoothing filtering (Ukil & Zivanovic, 2005b) is shown in the bottom section. Note that only the highpass filter output of the QMF pair is shown, hence the wavelet coefficients in the middle section indicate half of the total samples of the original signal (Ukil & Zivanovic, 2005b). The time-instants of the changes in the signal characteristics in the lower plot in Figure 8.1 indicate the different signal segments caused by different events during the fault. For example, segment A indicates the pre-fault section and the fault inception, segment B the fault, segment C the opening of the circuit-breaker, segment D the auto-reclosing of the circuit-breaker and system restore.
FIGURE 8.1: Abrupt Change Detection-based Segmentation of the BLUE-Phase Voltage Recording during a Phase-to-Ground Fault.
CHAPTER 8: APPLICATIONS
8.2
121
Synchronisation
Usually many DFRs, employed for different substations, trigger for any abnormal condition in the power network. All these simultaneous recordings only differ by some timedelays. First, all these simultaneous recordings are segmented using the abrupt change detection algorithms described in the previous chapters. However, before any further global analysis can be done, an important step is required, namely ‘Synchronisation’ (Ukil & Zivanovic, 2005c). This is essential because the time-bases of the DFRs triggering for the same disturbance are not perfectly synchronised (Chantler, Pogliano, Aldea et al., 2000) which can lead to erroneous conclusions especially when analysing the performances of the protective relays etc. Figure 8.2 shows the voltage recordings of three different DFRs triggering for the same phase-to-ground fault involving the BLUE-Phase. First, all the recordings are segmented using the abrupt change detection algorithms as described in the previous chapters. This is reflected in all of the three plots in Figure 8.2 by the dashed vertical lines, segmenting the voltage signals into segments like A, B, C and D. For all of them, segment A indicates the pre-fault section and the fault inception, segment B the fault, segment C the opening of the circuit-breaker, segment D the auto-reclosing of the circuit-breaker and system restore. It is clear from Figure 8.2 that all these segments are not synchronised for the three recordings, although they essentially represent the same event (Ukil & Zivanovic, 2005c).
122
CHAPTER 8: APPLICATIONS
FIGURE 8.2: Three Segmented but Unsynchronised DFR Voltage Recordings for the same Phaseto-Ground Fault.
To synchronise the recordings for further analysis, we will use the fault inception timing, that is, the borderline between segments A and B in all of the segmented recordings. It is important to note that in Figure 8.2, the X-axis indicates the number of samples. It is necessary to divide it by the sampling frequency of 2.5 kHz (Stokes-Waller, 1998) to obtain the fault inception time in milliseconds (Ukil & Zivanovic, 2005c). Table 8.1 lists the fault inception timings of the three DFR recordings shown in Figure 8.2.
CHAPTER 8: APPLICATIONS
123
TABLE 8.1: Fault Inception Time of the DFR Recordings shown in Figure 8.2
The complete synchronisation algorithm (Ukil & Zivanovic, 2005c) is described below.
First, we segment the different DFR recordings for the same disturbance using abrupt change detection.
Then we estimate the individual fault inception timings of the segmented but unsynchronised recordings.
We choose the recording with the minimum fault inception time as the reference one. (Here it is the DFR-3 recording as evident from Table 8.1.)
We synchronise the rest of the recordings with the reference recording by equating their fault inception times with the reference fault inception time; that is, we left-shift the rest of the recordings, their fault inception times equated with the reference one.
Then we again perform the abrupt change detection-based segmentation on these synchronised recordings to have the synchronised, segmented recordings for further analysis.
Application of the synchronisation algorithm on the unsynchronised recordings shown in Figure 8.2 results in the synchronised, segmented recordings as depicted in Figure 8.3. DFR-3 recording with the minimum fault inception time is chosen as the reference
124
CHAPTER 8: APPLICATIONS
recording and the other two recordings are synchronised accordingly (Ukil & Zivanovic, 2005c).
FIGURE 8.3: Three Segmented and Synchronised DFR Voltage Recordings for the same Phaseto-Ground Fault.
After synchronising the analogue signals (voltage signals in this case) of the different DFRs, the respective binaries are also synchronised and matched against the synchronised, segmented analogue signals (Ukil & Zivanovic, 2005c). One such example for the DFR-3 recording, analogue voltage and binaries for the fault duration and circuit-breaker auto-
CHAPTER 8: APPLICATIONS
125
reclosure is shown in Figure 8.4.
FIGURE 8.4: Synchronised Analogue Voltage Recording Plot and Binary Plots for the Fault Duration and Auto-Reclosing of the Circuit-Breaker.
After this, using the matched binary plots as a cross-check against any possible discrepancies, we estimate the change time-instants of the synchronised, segmented analogue recordings of the different DFRs (Ukil & Zivanovic, 2005c). Table 8.2 lists the change timeinstants of the recordings shown in Figure 8.3.
126
CHAPTER 8: APPLICATIONS
TABLE 8.2: Change Time-instants of the Synchronised, Segmented Analogue DFR Recordings shown in Figure 8.3
8.3
Relay Performance Monitoring
The synchronised, segmented analogue signals with the matched synchronised binaries can be used to monitor different performance characteristics of the protective relays, for example, fastest relay operating time, auto-reclosing of the circuit-breakers, and the main-1 and main-2 relay operation (Ukil & Zivanovic, 2005c).
8.3.1 Fastest Relay Operating Time To determine the fastest relay operating time, first we need to determine the fault duration. This can be done by estimating the duration of the segment B (fault) in the synchronised analogue plots shown in Figure 8.3, with the help of Table 8.2. Table 8.3 lists the fault duration times of the three DFR recordings shown in Figure 8.3 in terms of number of samples and milliseconds.
TABLE 8.3: Fault Duration in the Synchronised, Segmented Analogue DFR Recordings shown in Figure 8.3
CHAPTER 8: APPLICATIONS
127
To determine the fastest relay operation we need to select the minimum fault duration from the different synchronised DFR recordings (Ukil & Zivanovic, 2005c). In this case, this is the DFR-3 fault duration (70.4 milliseconds) from Table 8.3. The formula for calculating the fastest relay operating time is given below (Ukil & Zivanovic, 2005c). Fastest Relay Operating Time = Fault Duration – Trip Time.
(8.1)
Most of the circuit-breakers in the ESKOM transmission system are two-cycle breakers (Stokes-Waller, 1998), that is, the expected tripping time is in the region of 40 milliseconds (50 Hz system). Using this information and (8.1), we can compute the fastest relay operating time during the disturbance, which in this example case is 70.4 – 40 = 30.4 milliseconds. The fastest relay operating time gives a good idea whether the relays are operating correctly or whether they require maintenance (Ukil & Zivanovic, 2005c).
8.3.2 Auto-reclosing of the Circuit-Breakers From the synchronised, segmented analogue signals and their matched binaries, it is possible to analyse auto-reclosing of the circuit-breakers. By comparing the signal parameter values of segments A and D in the synchronised analogue plots as shown in Figure 8.3, it is possible to determine whether or not the autoreclosing is successful, following the relay operation (Ukil & Zivanovic, 2005c). In this case, matching the segment A and D signal parameters values of the synchronised analogue plots in Figure 8.3 shows that auto-reclosing of the circuit-breakers was successful. The length of the auto-reclosing can be determined by estimating the duration of the segment C in the synchronised analogue plots shown in Figure 8.3 (Ukil & Zivanovic, 2005c). Using the synchronised analogue plots of Figure 8.3 and Table 8.2 values, we can
128
CHAPTER 8: APPLICATIONS
compute the length of auto-reclosing (segment C) for all the three DFR recordings, as listed in Table 8.4.
TABLE 8.4: Circuit-Breaker Auto-Reclosing time of the Synchronised, Segmented Analogue DFR Recordings shown in Figure 8.3
8.3.3 Main-1 and Main-2 Relay Operation Every feeder in ESKOM transmission system is equipped with two identical relays, main-1 and main-2, with the same settings for distance protection (Stokes-Waller, 1998). The idea behind it is to avoid the possibility of a fault on the protection side when the incident occurs. Information on main-1 and main-2 relay operations can be obtained from the synchronised, matched binaries depending on the distance protection scheme. Analysis of main-1 and main-2 relay operation data gives some idea of whether or not the relays employed for the distance protection are working properly (Ukil & Zivanovic, 2005c). The relay operating times can be used to determine the relationship between the fault location and the tripping speed. Using this relationship and the system impedance ratio, it is possible to draw impedance plots for the relays, which can be used to analyse the distance protection (Ukil & Zivanovic, 2005c). Distance protection at ESKOM has three zones of protection, each with a certain impedance reach with respect to the impedance of the line: Zone 1 - 80%, Zone 2 - 120% and Zone 3 - 150% of the line impedance (Stokes-Waller, 1998).
CHAPTER 8: APPLICATIONS
8.4
129
Disturbance Analysis Examples
Abrupt change detection-based segmentation is the first step towards automatic fault recognition and disturbance analysis, followed by feature vector construction and pattern matching. In this section, we would discuss about the analysis of certain kinds of disturbances directly from the segmented recordings before conforming to any further significant and complex feature vector analysis.
8.4.1 Cleared Single Phase Fault Figure 8.5 shows the original DFR recording from the ESKOM transmission network for the voltage during a phase-to-ground fault in the RED-Phase. Segmentations of the voltage signal are indicated by the vertical dashed lines in Figure 8.5.
FIGURE 8.5: Segmented Analogue Voltage Recording during a Phase-to-Ground Fault Successfully Cleared.
130
CHAPTER 8: APPLICATIONS
In Figure 8.5, segment A indicates the pre-fault section and the fault inception, while segment D indicates auto-reclosing of the circuit-breaker and system restore. Comparison of the signal parameter values of segments A and D shows that the fault was successfully cleared (Ukil & Zivanovic, 2005d).
8.4.2 Uncleared Single Phase Fault By evaluating the induced voltage on the open phase, it is possible to see whether or not the secondary arc is extinguished. Very little voltage is required to sustain the secondary arc. Should the breaker close before the secondary arc is extinguished the fault will reappear and a three-Phase trip will be issued (Keller, 2004). Figure 8.6 shows an application example of how the uncleared fault can be evaluated using the abrupt change detection-based segmentation.
FIGURE 8.6: Segmented Analogue Current Recording during a Phase-to-Ground Fault not Successfully Cleared.
CHAPTER 8: APPLICATIONS
131
In Figure 8.6, the current during the phase-to-ground fault in the BLUE-Phase is shown along with the segmentations. Reappearance of the fault not cleared is evident from the comparison of segments B and D, representing the fault, and segment C representing the breaker opening and reclosing (Ukil & Zivanovic, 2005d).
8.4.3 Circuit-Breaker Restrike The restriking of a circuit breaker normally happens when the breaker has to interrupt capacitive current. A line breaker will interrupt capacitive current when it is the second breaker to open, which means that the remote end breaker is already open. During this condition, only charging current will flow, which is capacitive. During the switching of capacitive current, an extremely high voltage can be seen across the breaker contacts which results in the reappearance of the arc. The reappearance is called a restrike if it occurs 5 ms (¼ of a cycle for 50 Hz system) after breaker contact separation. The end that is restriking can be determined by examining the current waves. The voltage will cross through zero instantaneously and this will be associated with a spike in the current (Keller, 2004). Figure 8.7 shows the voltage and current recordings during a circuit-breaker restrike in the upper and lower plot respectively. Using the adjusted Haar wavelet method (Ukil & Zivanovic, 2005g) for abrupt change detection, both the voltage and current recordings are segmented into two segments, namely, A and B (i.e. before and after the restrike of the circuit-breaker). It is to be noted that the segments of the voltage and the current are matched. A closer look at segment B of the segmented current plot also shows the spike, indicating the restriking end (Ukil & Zivanovic, 2005d).
132
CHAPTER 8: APPLICATIONS
FIGURE 8.7: Voltage (Upper Plot) and Current (Lower Plot) Recordings during a Circuit-Breaker Restrike.
8.4.4 Reactor Ring Down When a line reactor is connected to a line and the circuit-breakers are opened at both ends, the voltage does not disappear. Instead, an oscillating voltage waveform can be found which slowly reduces in magnitude. This phenomenon is called reactor ring down (Keller, 2004). It is a result of the interaction between the reactor and the capacitance of the line. This forms an oscillatory circuit as depicted in Figure 8.8.
CHAPTER 8: APPLICATIONS
133
FIGURE 8.8: Schematic Oscillating Circuit Diagram of the Reactor Ring Down phenomenon.
Figure 8.9 shows the segmented voltage recording during a rector ring down incident. The segmentations are shown as vertical dashed lines.
FIGURE 8.9: Voltage Recording during a Reactor Ring Down.
Usually, many segments for the oscillating signals are recorded during the reactor ring down phenomenon as indicated in Figure 8.9. Comparing the signal parameter values of
134
CHAPTER 8: APPLICATIONS
segments B, C, D, E and F in Figure 8.9 indicates that the signal is oscillating and slowly decreasing in magnitude (Ukil & Zivanovic, 2005d).
8.4.5 Capacitive Voltage Transformer Transient Behaviour Capacitive Voltage Transformer (CVT) transient behaviour is caused by the discharge of energy stored in the capacitive and inductive elements of the CVT when there is a sudden change in the primary voltage. The transient behaviour is seen as oscillations in the secondary voltage. The transients are influenced by the burden on the CVT. For a resistive burden the transients are normally minute and die down rapidly. For zero or small burdens the transients are extremely prominent (Keller, 2004).
FIGURE 8.10: Segmented Voltage Recording for the Transient Behaviour of the Capacitive Voltage Transformer (CVT).
Figure 8.10 shows a segmented voltage recording for the transient behaviour of the CVT. From Figure 8.10, it should be noted that the transient behaviour is reflected in segment B.
CHAPTER 8: APPLICATIONS
135
For zero or small burdens the transients are extremely prominent as can be seen at the start of segments B and C (Ukil & Zivanovic, 2005d).
8.4.6 Energising of a Transformer Energising of a transformer often goes hand in hand with high magnetising inrush currents. Transformer protection must be set so that the transformer does not trip for this inrush current. This high current is the result of the remnant flux in the transformer core when it was switched out, and depends where on the sine wave the transformer is switched back in.
FIGURE 8.11: Segmented Current Recording reflecting the Energising of a Transformer, which is associated with High Magnetising Inrush Currents.
The segmented current recording for the energising of a transformer is shown in Figure 8.11. In Figure 8.11, the current signal is segmented into two segments, A and B, that is, before and after the transformer is energised. A closer look at segment B reveals the fact of high magnetising current inrush (Ukil & Zivanovic, 2005d).
136
CHAPTER 8: APPLICATIONS
8.5
Analysis of Power Signals from the Mexican Power Network
Ruiz-Vega, Messina and Enriquez-Harper (2005) discussed the use of nonlinear, nonstationary analysis techniques to characterise the forced inter-area oscillations problem in power systems, recorded in the Mexican interconnected system. In order to analyse the active power flow oscillations, they used nonlinear spectral representation of the data in the form of the wavelet transform and the Hilbert-Huang transform (Ruiz-Vega, Messina & Enriquez-Harper, 2005). They also compared the results with the linear spectral representation of the data in form of Fourier spectral analysis and Prony analysis (RuizVega, Messina & Enriquez-Harper, 2005). Figure 8.12 shows the active power flow oscillations recorded in the Mexican interconnected system (MZD-DGD) (Ruiz-Vega, Messina & Enriquez-Harper, 2005).
FIGURE 8.12: Active Power Flow Oscillations in the Transmission Line of the Mexican Interconnected System (MZD-DGD).
Based on their analysis, the recorded signal was divided into four main observation (time) windows. Each time window was then segmented into subintervals to investigate specific
CHAPTER 8: APPLICATIONS
137
characteristics of interests as shown in Figure 8.13.
FIGURE 8.13: Selected Time windows for Linear Spectral Analysis.
In collaboration with the Mexican team, we also tested our segmentation algorithms on the power oscillation signals obtained from the Mexican interconnected system (MZD-DGD). We tested our algorithms on the same signal shown in Figure 8.12. We tested the adjusted Haar wavelet algorithm and the wavelet algorithm (with the db4 and the db1 mother wavelet). Figure 8.14 shows the result of the test using the adjusted Haar wavelet method.
FIGURE 8.14: Segmentation of the Power Oscillation Signal from the Mexican Power Network, using the Adjusted Haar Wavelet Method.
138
CHAPTER 8: APPLICATIONS
Figures 8.15 and 8.16 show the results of the segmentation using the wavelet method based on the db4 and db1 mother wavelets respectively.
FIGURE 8.15: Segmentation of the Power Oscillation Signal from the Mexican Power Network, using the Wavelet Method (db4 mother wavelet).
FIGURE 8.16: Segmentation of the Power Oscillation Signal from the Mexican Power Network, using the Wavelet Method (db4 mother wavelet).
CHAPTER 8: APPLICATIONS
139
Comparison of the original time-window selection (Figure 8.13) and our signal segmentations (Figure 8.14 to 8.16), reveals that the adjusted Haar wavelet based approach (Figure 8.14) is reasonably accurate in proper signal segmentation for analysing the active power flow oscillations.
8.6
Summary
In this chapter we have discussed the various applications of the abrupt changes detectionbased segmentation on real disturbance signals. The primary application is event-specific segmentation for further signal processing towards automatic disturbance recognition. We proposed a novel synchronisation algorithm for the different DFRs triggering for the same disturbance. The synchronisation algorithm works on the basis of the fault-inception time estimated by the abrupt change detection algorithms. It selects the reference signal with the minimum fault-inception time and synchronises others with it by left-shifting them, their fault-inception times being equal. Using the synchronised and segmented analogue signals and their respective matched binaries, further effective analysis can be done for monitoring the performances of the protective relays, for example, the fastest relay operating time, auto-reclosing of the circuit-breakers, main-1 and main-2 relay operation. The analysis of certain kinds of disturbances directly from the synchronised, segmented recordings are also discussed in this chapter. This included cleared and uncleared single phase fault, circuit-breaker restrike, reactor ring down, capacitive voltage transformer transient behaviour, energising of a transformer. We also looked at signal segmentation on power oscillation signals from the Mexican interconnected system.
140
CHAPTER 9: COMMERCIAL IMPLEMENTATION
CHAPTER 9 COMMERCIAL IMPLEMENTATION
This thesis and work is part of the overall project entitled Automatic Disturbance Recognition in Electrical Power Systems. The whole project was initiated because of an industrial need, as most of the power utilities worldwide nowadays focus more and more on analysis of disturbances to prevent them. Synonymous with this is the commercial potential of the whole project. We looked at effective commercial implementation possibilities of the whole project utilising cutting-edge technologies. This chapter discusses the idea of the commercialisation, possible outcomes and the key technologies associated with it. Ukil and Zivanovic (2005h) provide a detailed discussion of the associated technologies.
9.1
Introduction
The analysis of faults and disturbances was in the past, and will continue to be in the future, a fundamental foundation for a secure and reliable electrical power supply. The introduction of digital recording technology opened up new dimensions in quantity and quality to fault and disturbance data acquisition. Nowadays, the challenge is to automatically convert data to knowledge, which frees the human resources to implement preventive action. If implemented, the automatic analysis system proposed and developed in the project entitled “Automatic Disturbance Recognition in Electrical Power Systems” will perform the mathematical analysis of all (presently known) patterns of incorrect behaviour (disturbances and faults) automatically and only short message would be issued for the control-
CHAPTER 9: COMMERCIAL IMPLEMENTATION
141
lers (operating personnel). This message would summarise the required information and knowledge regarding the faults and the corrective measures. The automatic analysis system would be available for use via the application service provider technology over the Internet as part of the project. The commercial implementation will be entitled as “Application Service Provider of Automatic Disturbance Analysis in
Power Systems”. This will have two modes of operations: Indirect and Direct. In the indirect mode of use, users will need to upload their fault signals via the Internet into the analysis system and receive the analysis report via the Web and e-mail. In the direct mode of use, users will be able to use the analysis system remotely (hosted in a central server computer) via the Internet after logging into the system using their login ID and password. All the application implementation, maintenance, upgrading and customisation of the analysis system will be done on the central host computer, which means that the users will not have to install any software in their computers. Instead, they will simply use it over the Internet at any time without any computational overheads.
9.2
Outcomes of Commercial Implementation
The commercial implementation focuses on developing a fully automatic disturbance analysis system for power systems and usage via the application service provider technology over the Internet. These two objective outcomes are explored in detail below.
9.2.1 First Objective Outcome: Automatic disturbance analysis System The automatic disturbance analysis system will be developed as the main research aim. The background and related technical information of the whole project are discussed in
142
CHAPTER 9: COMMERCIAL IMPLEMENTATION
detail in Chapter 2. In line with the first objective outcome of the commercial implementation, the automatic analysis system will be enhanced to make it more robust and customised for different users from a commercial point of view. The analysis system will accept the fault and disturbance signals from the users in the IEEE COMTRADE format, standardised and specified by IEEE STANDARD C37.111-1991. (1991), analyse the signals and return a report describing the fault and necessary corrective measures in accordance with the specific needs of the users.
9.2.2 Second Objective Outcome: Application Service Provider According to the second objective of the project, the automatic analysis system will be made available for use via the application service provider technology over the Internet. This will have two modes of operation: direct and indirect. In the direct mode of use, users will be able to use the analysis system remotely (hosted in a central server computer) via the Internet after logging into the system using their login ID and password as in any other Website. In the indirect mode of use, users will need to upload their fault signals via the Internet into the analysis system and will receive the analysis report via the Web and Email. Hence the definite outcomes of the project will be the automatic disturbance analysis system as an application service over the Internet. The idea of an application service provider and the key technologies associated with it will be described in the following sections.
CHAPTER 9: COMMERCIAL IMPLEMENTATION
9.3
143
The Application Service Provider (ASP)
The Application Service Provider (ASP) solution offers individuals or enterprises access to software applications and related support services over the Internet. Obtaining these applications, services, solutions from an outside supplier is a cost-effective solution to the demands of systems ownership: up-front capital expenses, implementation challenges, and a continuing need for maintenance, upgrading and customisation. Thus ASP is a promising technology ideal for a cost-effective, competitive information technology marketplace (Ukil & Zivanovic, 2005h). ASP is a third-party entity that manages and distributes software-based applications, services and solutions to customers across a wide area network (WAN) such as the, Internet from a central data server. Basically, ASPs are Independent Software Vendors (ISVs) or Internet Service Providers (ISPs) that are now using the Internet as the delivery vehicle to make software applications available (Ukil & Zivanovic, 2005h). Delivering access to applications in this manner allows small to medium enterprises to eliminate the time and costs associated with installing, managing and supporting new applications.
9.4
Overview of ASP
ASP provides a contractual service offering to deploy, host, manage and rent access to an application from a centrally managed facility. In the following sections, we discuss various aspects of ASP.
9.4.1 ASP Categories Generally ASPs are broken down into five categories as specified by ASPstreet.com:
144
CHAPTER 9: COMMERCIAL IMPLEMENTATION
•
Enterprise ASPs which deliver high-end business applications
•
Local/Regional ASPs which supply wide variety of application services for smaller businesses in a local area
•
Specialist ASPs which provide applications for a specific need, such as Web site services or human resources
•
Vertical Market ASPs which provide support to a specific industry, such as healthcare
•
Volume Business ASPs which supply general small/medium-sized businesses with pre-packaged application services in volume.
9.4.2 Benefits There are many benefits to renting applications over the Internet rather than maintaining the software and hardware on site. Some of them are mentioned below as discussed by Ukil and Zivanovic (2005h): •
The problem of the shortage of information technology (IT) employees is resolved.
•
Use of the ASP model means that, there is no need to make large up-front payments for software licences and hardware. Instead, a monthly subscription fee is paid and the IT costs can be spread over time.
•
Access to high-end applications is possible which otherwise might be unaffordable.
•
There is guaranteed performance and uptime.
CHAPTER 9: COMMERCIAL IMPLEMENTATION
•
Implementation time is faster and the product cycle shorter.
•
Security is enhanced.
145
9.4.3 Key Features Currently, many vendors are trying to ASP-enable their applications so that they can be delivered through the ASP model. Enabling them in this way is sometimes a difficult task and this process must conform to certain features of ASPs for successful ASP operation. Some of these are mentioned below on the basis of Ukil and Zivanovic’s work (2005h).
Type of Customer
Whether the user is a large, medium or small company, or even a home-user, an ASP arrangement saves time and money by eliminating massive investments in software, deployment time and IT personnel.
Reliability
Reliability is a key issue and customers who turn to ASPs should find better application reliability and availability compared with their experience with their internal IT organisations.
Cost Saving
This is one of the key factors for successful ASPs. Studies by ASPstreet.com have indicated that by leasing an application from an ASP, customers save between 33 and 53% compared with purchasing and managing the hardware and software for the application themselves.
146
CHAPTER 9: COMMERCIAL IMPLEMENTATION
Flexible Payment scheme and Agreement
For economic use of the service, several payment models should be offered, such as, payment per usage/CPU usage: computing power/per time/per result etc. This should be done according to structured service-level agreement (SLA) (Ukil & Zivanovic, 2005h) which is a contractual obligation between an ASP and its clients. The SLA should spell out the details of service, such as allowable downtime, connection speed, application access, security issues etc.
Kinds of Application
Applications which are difficult to procure and maintain are generally well suited to ASPs because companies prefer to rent rather than procure them. The applications should be Web-enabled, profitable and of interest to clients.
Availability
Another important issue related to ASP enabling of applications is guaranteed availability. The application and its infrastructure should always be available around the clock, seven days a week under unusually high loads and other issues.
Security
Because ASPs depend on the Internet, proper security is vital. A secure access method should be offered for usage, maintenance as well as for data transfer with strong data encryption technology (secure protocol such as, SSH, IPSec and Virtual Private Network [VPN] and the like).
Load Balancing
During high peak time, the same performance should be guaranteed at the same quality level.
CHAPTER 9: COMMERCIAL IMPLEMENTATION
147
Multi-user capability
The solution should be able to handle multiple remote user connectivity simultaneously.
Adequate Bandwidth
The bandwidth should be adequate to support the required number of remote operations.
Universal operation
There should be no dependency on the operating system unless the aim of the ASP is for platform-specific application(s).
Strong Customer Service
Because ASP is Web-based and replaces in-house software applications, customer service is a critical factor. Customer service should be technically strong, flexible, available round the clock and in various modes (Internet based, telephonic, fax, etc).
High-quality Management
Top quality management is required for both technical and administrative issues to ensure smooth and successful operation.
9.4.4 Key Technologies The key technologies for ASPs considered for this commercial implementation are: •
Web Services
•
Thin Client Computing: Tarantella®, Windows® Terminal Services and Citrix®
•
WinConnect® Server XP™ technology.
148
CHAPTER 9: COMMERCIAL IMPLEMENTATION
These are discussed in detail in the following sections.
9.5
Web Services
This section describes the Web services in detail.
9.5.1 Introduction A Web Service provides service at some level, and is primarily used as a means for businesses to communicate with each other and with clients. Web services allow organisations to communicate data without intimate knowledge of each other's IT systems behind the firewall and describe a standardised way of integrating Web-based applications using the XML, SOAP, WSDL and UDDI open standards (SUN MICROSYSTEMS INC.) over an Internet protocol backbone. The primary advantage of a Web service is the low cost in investment compared with explicit licensing of software. Also, it eases the information management and staff involvement because the Web service provider is responsible for software updates, patches, version control, security concerns, back-up policy and so forth. A Web service provides an open module concept, that is, the possibility for further extension of the service as required (Ukil & Zivanovic, 2005h). Hence the main concerns regarding Web service are data control, security of operation, sharing of responsibility, and so forth. These can be managed by handling access via strong authentication, access control lists, shadowed password database etc (SUN MICROSYSTEMS INC.).
CHAPTER 9: COMMERCIAL IMPLEMENTATION
149
9.5.2 Web Service Technology Web services, unlike the traditional client-server models, such as a Web server-Web page system, do not provide the user with a graphical user interface (GUI). Instead, Web services instead share business logic, data and processes through a programmatic interface across a network. Web services allow different applications from different sources to communicate with each other without time-consuming custom coding, and because all communication is in extensible markup language (XML), Web services are not tied to any one operating system or programming language. For example, Java™ can talk with Perl™, Windows® applications can talk with UNIX® applications. Also, it is possible to use secure protocols for data transmission e.g., SSH2, HTTPS etc for Web service connectivity. Tunnelling of the corporate network is possible as well as terminal emulation is also possible for Web and X-Windows® based applications (SUN MICROSYSTEMS INC., 2001).
9.5.3 Web Service Architecture Web services communicate with clients by exchanging XML documents. The use of the XML standard facilitates the clients to retrieve data from a Web service without having knowledge of the technology underlying the data source. Clients connect to a Web servicepublished data source, execute a request, and receive response data formatted as an XML document, but are unaware of underlying protocols, component models, application programming interfaces (APIs), and operating systems as shown in Figure 9.1 (SUN MICROSYSTEMS INC., 2001).
150
CHAPTER 9: COMMERCIAL IMPLEMENTATION
FIGURE 9.1: Web-Service Architecture: XML-based Interface to Back-end Services.
The presentation server is a Web application that transforms XML documents into documents that are appropriate for the client type, for example, hypertext markup language (HTML) or wireless markup language (WML) documents. Internally, the Web service architecture depends on the Request-Response mechanism (SUN MICROSYSTEMS INC., 2001) as described below. When a Web service receives a client request, it forwards the request in the form of an XML document (the XML input document) to the appropriate XML operation. The XML operation calls one or more methods on business components. The XML operation transforms the return values of these method calls into an XML document (the XML output document) and returns it (Ukil & Zivanovic, 2005h). When an XML operation is executed, the Web service does the following: (1) It parses the XML input document, mapping the document’s elements to the parameters of the methods that the XML operation is defined to call. (2) It calls the methods defined in the XML operation, in their specified order.
CHAPTER 9: COMMERCIAL IMPLEMENTATION
151
(3) It formats return values of the methods into an XML output document according to the definition of the XML operation. (4) It returns the XML output document. Request-Response mechanism is depicted in Figure 9.2.
FIGURE 9.2: Simplified Service-Client Request-Response Mechanism.
9.5.4 Web Service Building Blocks According to standards specified in the Web service arena, Web service basic building blocks are as follows: •
Describing Web Services: Web Service Description Language (WSDL)
WSDL (WORLD WIDE WEB CONSORTIUM (W3C), 2001) describes not only the possible operations (including messages, parameters, return values and complex types) but also how to access them (binding information, ports, port Types). WSDL documents closely resemble remote procedure call (RPC), components or distributed objects, interface description language (IDL), but are defined in XML.
152
•
CHAPTER 9: COMMERCIAL IMPLEMENTATION
Publishing Web Service Metadata: Universal Description, Discovery and Integration (UDDI)
The Universal Description, Discovery and Integration (UDDI) standard (ORGANISATION FOR THE ADVANCEMENT OF STRUCTURED INFORMATION STANDARDS [OASIS] CONSORTIUM) defines how to publish business information and associated technical service descriptions relating to Web services. •
Transporting Information: Simple Object Access Protocol (SOAP)
The Simple Object Access Protocol (SOAP) (WORLD WIDE WEB CONSORTIUM (W3C), 2000) is generally regarded as the underlying transport for Web services. It carries messages and their parameters encoded as XML and can send messages synchronously (over hypertext transfer protocol [HTTP]) or asynchronously (over email, file transfer protocol [FTP] or mobile multicast protocol [MoM]). The SOAP specification defines the XML tags used to delimit the message structure and header elements. •
Extensible Markup Language (XML), Electronic Business XML (ebXML)
XML is a set of rules for defining semantic tags that break a document into parts and identify the different parts of the document. It is a meta-markup language that defines a syntax used to define other domain-specific, semantic, structured markup languages (Harold, 1999).
9.5.5 Web Service Development Forte™ for Java™ (SUN MICROSYSTEMS INC., 2001) integrated development environment (IDE) module provides powerful technology that makes possible the publication of back-end business services over the World Wide Web as Web services. On the server side, the user needs to create a language binding (SUN MICROSYSTEMS INC., 2001), imple-
CHAPTER 9: COMMERCIAL IMPLEMENTATION
153
ment the functionality and then register or create registry information. On the client side, the user needs to create a language proxy (SUN MICROSYSTEMS INC., 2001) and then call methods or send messages.
9.6
Thin Client Computing
Because of the ongoing support costs for Client-Server technology, many companies have expressed a wish to return to economies of the ‘expensive’ centralised computing environment they had previously abandoned. The concept of deploying Java™ programs to simple network attached computers running software from a central server was developed. Java™ would provide a modern, colourful, windowing user interface at the workstations, and the workstations would have their own computing power for good performance. All the software would live back on the central server where it could be easily and inexpensively maintained. This is the concept of modern Thin Client Computing (Ukil & Zivanovic, 2005h). Three different thin client computing solutions are discussed in the following sections.
9.6.1 Tarantella® Tarantella® (TARANTELLA INC., 2000) software provides fast, secure remote access to all server-based applications, including those hosted on Microsoft® Windows® NT Terminal Services, Microsoft® Windows® 2000 servers, all UNIX® and Linux servers, Mainframes and AS/400® servers. Using standard graphical and character protocols, such as Microsoft® RDP, X11, VT, TN3270 and TN5250, Tarantella® software ensures full compatibility with all application servers, now and in the future.
154
CHAPTER 9: COMMERCIAL IMPLEMENTATION
It allows users to quickly access all required applications without resorting to development work. It allows employees, customers or partners to have fast application access, locally or remotely, all through the simplicity of a Web browser. When a user logs on to the Tarantella® server (just like accessing any other Web page on the Internet), Tarantella® looks to see if user has the client Java™ application. If not, it is automatically downloaded to the workstation (only about 280KB). The company's Web applications menu then loads into the browser screen. The user can start any applications offered by the company by simply pointing to the icon and clicking (TARANTELLA INC., 2000). All administrative updates are managed centrally on Tarantella® servers, by designated Tarantella® administrators. Without visiting either the application servers or the client devices, administrators can publish applications to users and make them immediately available for use. Administration uses the directory service model, with objects representing people and resources within an organisation, arranged in a hierarchy that mirrors the structure of the organisation. Each user has his/her own Webtop, namely, the application launchpad (TARANTELLA INC., 2000). The Webtop travels with the user that is, it does not remain with the client device. Users can access their own personalised Webtop from any client device with a network connection to the Tarantella® server. Applications run on servers, not on client devices, and can be made resumable. This means that users can log out of their Webtop without exiting the application. Later, they can resume the same application from any client device, and pick up where they left off. The working principle of Tarantella® is depicted in Figure 9.3 (TARANTELLA INC., 2000).
CHAPTER 9: COMMERCIAL IMPLEMENTATION
155
9.6.2 Windows® Terminal Services Windows® Terminal Services (WTS) (MICROSOFT CORPORATION, 1999) is Microsoft's® thin client strategy. One popular use of WTS is running, say a centralised Windows® accounting program, and using WTS to serve up satisfactory performance at remote offices. This avoids having to start over with a new UNIX® based or Client-Server accounting program when a business opens new branches.
FIGURE 9.3: Tarantella® Working Principle.
9.6.3 Citrix® Citrix® Metaframe® XP™ (CITRIX SYSTEMS INC.) presentation server centrally manages enterprise applications and accesses them from anywhere. Citrix® Metaframe® XP™ is the easiest way to manage enterprise applications from a central location and access
156
CHAPTER 9: COMMERCIAL IMPLEMENTATION
them from anywhere. The foundation of the Citrix® Metaframe® Access Suite, Metaframe® XP™ presentation server is the world’s most widely deployed presentation server for centrally managing Windows®, Web and legacy applications and delivering their functionality as a service to workers, wherever they may be. Metaframe® XP™ presentation server is certified to run on Windows® 2000 and 2003 servers, and supports virtually any custom or commercially packaged Windows® applications.
9.7
WinConnect® Server XP™
WinConnect® Server XP™ (THINSOFT PTE LTD, 2002) enables a Windows® XP™ computer (Host PC) to allow up to 21 remote desktop sessions. It allows Remote Desktop Protocol (RDP) (MICROSOFT CORPORATION, 2000) 4.0, 5.0 and 5.1-enabled Thin Client devices (such as Terminals, Internet/Information Appliances, Tablet PCs and Personal Digital Assistants [PDAs]) to connect to a Host PC to run Windows® applications simultaneously and independently. WinConnect® takes advantage of Microsoft's® Remote Desktop Connection (RDC) (MICROSOFT CORPORATION, 2000) technology, which allows a user on a remote PC to log on to a host Windows® XP™ system. Such a remote user can use applications, printers and Internet connections running on the host, and store data on the host's hard drive exactly as if the user were sitting at the host system. The remote system does not have to be a purpose-built thin client device. It can be any PC running RDC client software - including Linux and even DOS® systems. The host can be any system running Windows® XP™, Home or Professional - although there are some restrictions with Home. The means of connection can be a local area network (LAN) or the
CHAPTER 9: COMMERCIAL IMPLEMENTATION
157
Internet. A broadband Internet connection over a digital subscriber line (DSL, ADSL) or Integrated Services Digital Network (ISDN) connection improves performance, but is not necessary. This is because Remote Desktop transfers only the minimal data (such as display data and keyboard data) to remotely operate with the Host computer. Therefore, even low-bandwidth Internet connections or a wired or wireless TCP/IP connection, LAN, WAN, dial-up (Internet), VPN connection will suffice (THINSOFT PTE LTD, 2002). Depending on the type and licensing of software to be used with WinConnect®, all users can use and work with the same software at the same time. Files can be stored on the central WinConnect® Server XP™ Host or on the local storage of the Remote Desktop Client (if the Remote Desktop Client is a RDP 5.1-enabled device only and the WinConnect® Server XP™ Host is installed with Windows® XP™ Professional Edition) (Ukil & Zivanovic, 2005h). Data can be stored in the central WinConnect® Server XP™ Host instead of the local PC hard drive. Normally, the server is located in a secure place and is only accessible by IT personnel who have the necessary rights. The data transmitted from the Remote Desktop Client, usually, takes a path along a secured communication/network, to the WinConnect® Server XP™ Host. The working principle of WinConnect® is depicted in Figure 9.4 (THINSOFT PTE LTD, 2002).
158
CHAPTER 9: COMMERCIAL IMPLEMENTATION
FIGURE 9.4: WinConnect® Server XP™ Working Principle.
9.8
Application in Commercial Implementation
In this section, we will discuss the utilisation of the various ASP technologies in a comparative manner for the commercial implementation of the ASP of Automatic disturbance analysis in Power Systems. The automatic analysis system would be available for use via the ASP technology over the Internet in two modes. In the indirect mode of use, users will need to upload their fault signals via the Internet into the analysis system and will receive the analysis report via the Web and E-mail. In the direct mode of use, users will be able to use the analysis system remotely (hosted in a central server computer) via the Internet after logging into the system using their login ID and password. It will be necessary to Web-enable and integrate various systems in the indirect mode, for example, the fault analysis engine, the Web interface, feedback system etc, which are
CHAPTER 9: COMMERCIAL IMPLEMENTATION
159
cross-platform from the point of view of the operating system and programming languages. Hence the Web service would be the perfect technology to integrate and implement the indirect mode of the ASP of the automatic disturbance analysis system (Ukil & Zivanovic, 2005h). The direct mode requires extensive and robust remote connectivity and operation. Hence for cross-platform remote operation, we propose the use of Thin Client Computing technology, namely, Tarantella® tool for realising the direct mode of the ASP. By contrast, for the predominant Windows® based remote operation, WinConnect® Server XP™ would be a better suitable tool to realise the direct mode of the ASP, from the point of view of operation and implementation complexity and cost-effectiveness (Ukil & Zivanovic, 2005h). Application of the ASP technologies in the commercial implementation of the project, in the direct and indirect mode is depicted in the Figure 9.5.
FIGURE 9.5: Commercial Implementation: ASP of Automatic Disturbance Analysis in Electrical Power Systems.
160
9.9
CHAPTER 9: COMMERCIAL IMPLEMENTATION
Highlights and Advantages of Commercial Implementation
The key highlights and advantages of the commercial implementation of the project and its outcomes are outlined below. •
Innovation and Technology Transfer
The research and development of the automatic disturbance analysis system is a direct innovative breakthrough in the power industry, to make possible secure, fast, reliable and fully automatic fault recognition and correction. The commercial implementation of the project would initiate a direct technology transfer to the industry. For example, the developed system can be directly applied to power utilities, such as, ESKOM, South Africa, augmenting its present semi-automatic fault and disturbance analysis system. •
Prototype Product and Business Development
Because of the direct industrial application and commercial availability through cuttingedge ASP technology, the commercial implementation of the project could be successfully further extended to real business development with the proposed outcomes, market and competitor analysis, financial, development and operational planning, being implemented in parallel with the technological aspects. •
Job Creation
Since this project has direct industrial application and true business potential, its commercial implementation would certainly create numerous sustainable job opportunities. •
Economic Development
The industry-oriented, advanced technological outcomes of the commercial implementati-
CHAPTER 9: COMMERCIAL IMPLEMENTATION
161
on would certainly contribute towards an all-round economic development, from job creation to better life, being a fundamental research foundation of the electrical power industry which is a key factor in today’s successful economy and development. •
Intellectual Capital Development and Strengthening of Institutional Resources
The commercial implementation, its proposed outcomes, execution and technology would be a world-class research and development endeavour with off-the shelf industry application. This would signify high-profile intellectual capital for the institution which could go a long way in strengthening the institutional resources for further innovative and industryspecific research and development.
9.10 Summary We have discussed about the potential commercial implementation of the whole project as an application service provider (ASP) solution in this chapter. Impregnated with real industrial need for the power utilities world-wide, the Automatic Disturbance Recognition in Electrical Power Systems project could initiate successful technology transfer to the power industry. We have investigated and presented the idea of the commercialisation, possible outcomes, effective benefits, key technologies for the ASP implementation and utilisation of the ASP technologies for the commercial implementation of our project. The automatic disturbance analysis system was developed as the main research aim, the first objective outcome of the commercial implementation of the project. The second objective of the project, the automatic analysis system would be made available for use via the application service provider technology over the Internet in direct and indirect modes.
162
CHAPTER 10: CONCLUSION
CHAPTER 10 CONCLUSION
In this thesis a new approach to event-specific signal segmentation for the purpose of automatic disturbance recognition in electrical power systems is provided. Different novel segmentation algorithms were proposed and developed. All the algorithms were tested on real disturbance recordings obtained from digital fault recorders of the South African power network. It was shown that the event-specific segmentation facilitates and augments further signal processing and disturbance recognition and analysis tasks. The method of synchronisation of the segmented disturbance signals was also developed. Various disturbance analysis tasks, relay performance monitoring directly from synchronised, segmented disturbance signals were performed. Also, we looked into the commercialisation of the project as Application Service Provider solution.
10.1 Automatic Disturbance Recognition in Electrical Power Systems The analysis of faults and disturbances in power systems is a fundamental foundation for a secure and reliable electrical power supply. In the project, Automatic Disturbance Recognition in Electrical Power Systems, we focus on automated disturbance recognition and analysis for the power transmission network of South Africa, in collaboration with South African power utility, ESKOM. The purpose of this study is to augment the existing fault analysis and recognition system at National Control, ESKOM, with more robust and accurate algorithms and automated techniques in an effort to make it fully automated. This could streamline the complex task of converting massive amounts of disturbance data into
CHAPTER 10: CONCLUSION
163
the knowledge which is extremely complex, time-consuming and slow when done manually. So, in this thesis we have highlighted the urgent need for the automatic disturbance recognition in the power industry. The problem of lack of fast, fully automated processing tools for the massive amount of disturbance records available from advanced instruments such as the digital fault recorders was identified. Following the industrial motivation for the whole project and hence proper specification, we described the idea of an efficient, robust and advanced automated analysis system in this thesis. Then we looked into the detail architecture of the proposed automated analysis system and sub-divided the overall project of Automatic Disturbance Recognition in Electrical Power Systems into four major submodules. Typically the first processing step in a recognition-oriented signal processing can be automatic segmentation of the signal. A similar approach was proposed in this thesis by performing the event-specific signal segmentation based on detection of abrupt changes in the disturbance signal model parameters. This is the first step towards the proposed automatic disturbance analysis and primary focus of this thesis, followed by feature vector construction for the specific segments and applying pattern-matching algorithms to accomplish the automatic disturbance recognition and analysis tasks. Judging the real commercial potential of the whole project, we discussed in this thesis a commercial implementation of the fully automatic fault analysis system for power systems and usage of that via the application service provider technology over the Internet.
164
CHAPTER 10: CONCLUSION
10.2 Abrupt Change Detection based Signal Segmentation The focus of this thesis, abrupt change detection is used as a segmentation algorithm for disturbance signals. Detection of abrupt changes in signal characteristics is a much studied subject with many different approaches. We first investigated and categorised the different approaches of abrupt change detection-based segmentation for this research. Following those approaches, we developed the specific algorithms of abrupt change detection based signal segmentation and tested them on the real disturbance signals.
10.2.1
Abrupt Change Detection Algorithms
In this thesis, we categorised the abrupt change detection techniques into four classes: simple methods, the linear model based approach, the model free approach and the nonparametric approach. From a practical point of view, we also classified the disturbance signals into two broad categories: signals with and without distinct abrupt changes in the signal parameters. In this thesis, we developed two main segmentation algorithms. One is based on linear state-space modelling and utilises recursive parameter identification technique. The other is based on the nonparametric approach, using the wavelet transform and universal thresholding method. Of these two algorithms, the wavelet transform-based algorithm is more robust for various kinds of disturbance signals. These two algorithms are well suited to different kinds of disturbance signals, especially for signals showing distinct abrupt changes in the signal parameters. For the signals not showing distinct changes in signal parameters, we developed a new adaptive abrupt change detection technique. This consists of two methods, namely the ada-
CHAPTER 10: CONCLUSION
165
ptive whitening filter method and the adjusted Haar wavelet method. The adaptive whitening filter method acts as an improvement pre-filter to be used in conjunction with the wavelet method. In the second method, we developed a general adjusting technique to the standard Haar wavelet. This generates a series of new adjusted Haar wavelets which we successfully used as the basis for the wavelet transform algorithm for accomplishing the adaptive abrupt change detection task. Finally, we combined all the developed algorithms into an optimised overall algorithm which is well-suited for the event-specific segmentation task for all kinds of disturbance signals. Practical application results show about 99% accuracy which is fairly satisfactory. Long, medium and short lines do not have significant impacts on the segmentation algorithms. Increases in the fault levels at the substation buses due to commissioning of new lines are tackled by the adaptive algorithm, hence do not influence the overall segmentation algorithm. VAR equipment switching has no effect on the algorithms either.
10.2.2
Applications of Abrupt Change Detection Algorithms
Development of the abrupt change detection-based segmentation algorithms was properly utilised in the scope of this work and thesis for many different applications. Of primary significance is the event-specific segmentation of the disturbance signals for further signal processing towards automatic disturbance recognition and analysis. This helps us to identify automatically the specific, finite duration signal segments, such as, fault duration, circuit-breaker opening, auto-reclosing, etc. From the automated recognition-oriented signal processing point of view, this step improves the accuracy, speed and robustness of the whole automatic recognition task. This step would allow us to concen-
166
CHAPTER 10: CONCLUSION
trate on and process the specific segments indicated by the segmentation, instead of the whole signal. In this thesis, we also developed a new synchronisation algorithm for the different simultaneous unsynchronised disturbance recordings. This is essential because the time-bases of the DFRs triggering for the same disturbance are not perfectly synchronised. This could lead to erroneous conclusions, especially when analysing the performances of the protective relays and so forth. We also looked at different applications directly from the synchronised, segmentation without conforming to any further complex feature vector analysis. This included monitoring the relay performances: estimating the fastest relay operating time, analysing autoreclosing of the circuit-breakers and analysing the performance of the main-1 and main-2 relays. We also looked at analysis of certain kinds of disturbances, which included cleared, uncleared single phase fault, circuit-breaker restrike, reactor ring down phenomenon, transient behaviour of the capacitive voltage transformer and energising of a transformer. Beside the applications that were focused on the power network in South Africa, we also applied our abrupt change detection based algorithm for the segmentation of the power oscillation signals obtained from the Mexican interconnected system and achieved satisfactory result.
10.3 Future Work In this section, we review some future research directions as discussed by Ukil and Zivanovic (2006) extending the work presented in this thesis. These are mentioned below.
CHAPTER 10: CONCLUSION
10.3.1
167
On-line Abrupt Change Detection in Electrical Power Systems
The abrupt change detection algorithms presented in this work are intended for off-line operation. However, the accuracy, speed of operation of the individual and the complete segmentation algorithm(s) is/are deemed to be fairly favourable. Hence a possible future endeavour could be to apply them for on-line signal segmentation in the electrical power systems. This would also be the first step towards on-line disturbance recognition and analysis.
10.3.2
Early Disturbance Prediction and Prevention
Abrupt change detection-based segmentation is critical for successful and reliable automated disturbance recognition and analysis, as shown by the various applications in this thesis. However, the whole analysis is done only after the disturbance or fault has actually occurred. This can also be extended along with on-line detection procedures, to detect any anomaly early in the power network signals, which could lead to certain system instability, system overloading, disturbances and the like. Prediction of the possibility of system instability, deterministic disturbances beforehand could be critical in preventing blackouts. It is to be noted that Thottan and Ji (2003) used abrupt change detection, statistical analysis and pattern matching together for anomaly detection in IP networks in order to predict and prevent anomalies before they occur.
10.3.3 Application in other Domains The abrupt change detection and automatic segmentation technique, can be successfully applied to other domains, for example, segmentation of the music signals for automatic music cognition, blind source separation and so on.
168
BIBLIOGRAPHY
BIBLIOGRAPHY
ANDERSSON, P. 1985. Adaptive forgetting in recursive identification through multiple models. International Journal of Control, 42(5): 1175 – 1193.
ASPstreet.com. White papers on Application Service Provider. Available: http://www.aspstreet.com
BARTH, P., LUDWIG, A. and SCHEGNER, P. 2003. Development of the Fault and Disturbance Data Handling in the German HV-network. CIGRE SC-34 Colloquium, Sydney.
BARTYLAK, A. 2002. Application of Disturbance Recorders as Near Real Time Information Support for National Control in ESKOM. In Proceedings of DPSP Conference, Amsterdam, Netherlands.
BASSEVILLE, M. 1988. Detection change in signals and systems – A survey. Automatica, 24: 309 – 326.
BASSEVILLE, M., ESPIAU, B. and GASNIER, J. 1981. Edge detection using sequential methods for change in lever – Part I: A sequential edge detection algorithm. IEEE Transactions on Acoustics, Speech, Signal Processing, ASSP-29(1): 24 – 31.
BASSEVILLE, M. and NIKOFOROV, I.V. 1993. Detection of Abrupt Changes – Theory and Applications. Englewood Cliffs, NJ : Prentice-Hall.
BIBLIOGRAPHY
169
CABALLERO, I.S., PRIETO, C.P. and VIDAL, A.R.F. 1997. Joint Segmentation and AR Modeling of Quasistationary Signals using the EM Algorithm. In Proceedings of IEEE Workshop on Nonlinear Signal and Image Processing (NSIP’97), Mackinac Island, Michigan.
CHANTLER, M., POGLIANO, P., ALDEA, A., TORNIELLI, G., WYATT, T. and JOLLEY, A. 2000. The Use of Fault-Recorder Data for Diagonising Timing and Other Related Faults in Electricity Transmission Networks. IEEE Transactions on Power Systems, 15(4): 1388 – 1393.
CHUI, C. K. 1992. Wavelets: a tutorial in theory and applications. New York: Academic Press.
CITRIX SYSTEMS INC. Citrix® Metaframe® white papers and guide. Ft. Lauderdale, Florida. Available: http://www.citrix.com
CRAIGMILE, P. F. and PERCIVAL, D. B. 2000. Wavelet-Based Trend Detection and Estimation. Department of Statistics, Applied Physics Laboratory, University of Washington, Seattle, WA.
DAUBECHIES, I. 1992. Ten Lectures on Wavelets. Philadelphia: Society for Industrial and Applied Mathematics.
DESOBRY, F. and DAVY, M. 2003. Support Vector-Based Online Detection of Abrupt Changes. In Proceedings of IEEE ICASSP Conference, Hong-Kong, China.
170
BIBLIOGRAPHY
DONOHO, D.L. and JOHNSTONE, I.M. 1994. Ideal Spatial Adaptation by Wavelet Shrinkage. Biometrika, 81(3): 425 – 455.
FRANKLIN, G.F., POWELL, J.D. and WORKMAN, M.L. 1990. Digital Control of Dynamic Systems. 2nd Edition. Reading, MA: Addison-Wesley.
GELEB, A. and VELDE, W.E.V. 1968. Multiple-input Describing Functions and Nonlinear System Design. New York: McGraw-Hill.
GIRSHICK, M.A. and RUBIN, H. 1952. A Bayes approach to a quality control model. Annals Mathematical Statistics, 23: 114 – 125.
GRADSHTEYN, I. S. and RYSHIK, I. M. 1965. Table of Integrals, Series, and Products, New York: Academic Press.
GRETTON, A. and DESOBRY, F. 2003. On-line one-class ν–support vector machines. An application to signal segmentation. In Proceedings of IEEE ICASSP Conference, Hong-Kong, China.
HAAR, A. 1910. Zur Theorie der orthogonalen Funktionen-Systeme. Math. Ann., 69: 331 – 371.
HANNAN, E.J. 1980. Recursive estimation based on ARMA models. Annals of Statistics, 8: 762 – 777.
BIBLIOGRAPHY
171
HAROLD, E.R. 1999. XML™ Bible. Foster City, California: IDG Books Worldwide.
HASTIE, T., TIBSHIRANI, R. and FRIEDMAN, J. 2001. The Elements of Statistical Learning. New York: Springer.
HAYKIN, S. 1996. Adaptive Filter Theory. Englewood Cliffs, NJ: Prentice-Hall.
HENZE, C. 2005. Automatic Fault Classification using Support Vector Machines. Practical Training Report. Eskom Transmission, Germiston, South Africa (Unpublished).
IEEE STANDARD C37.111-1991. 1991. IEEE Standard Common Format for Transient Data Exchange. Version 1.8.
JAZWINSKI, A.H. 1969. Adaptive filtering. Automatica, 5: 475 – 485.
JORDAAN, J.A. and ZIVANOVIC, R. 2002. Nonparametric Time-Varying Phasor Estimation in Power Systems by making use of Non-quadratic Criterion. In Proceedings of IEEE AFRICON’02, South Africa.
KALMAN, R.E. 1958. Design of a self-optimizing control system. Transaction ASME, 80: 468 – 478.
KELLER, P. 2004. Analysis of Digital Records. Technical Report. Eskom Transmission, Germiston, South Africa (Unpublished).
172
BIBLIOGRAPHY
KELLER, P., HENZE, C. and ZIVANOVIC, R. 2005. Automated Analysis of Digital Fault Records: A Utility Perspective. In Proceedings of Southern African Universities Power Engineering Conference (SAUPEC’05), Johannesburg, South Africa.
KUNG, S.Y., ARUN, K.S. and RAO, D.V.B. 1983. State-Space and Singular Value Decomposition-Based Approximation Methods for the Harmonic Retrieval Problem. J. Opt. Soc. Am., 73(12): 1799 – 1811.
LAΪ, T.L. 1974. Control charts based on weighted sums. Annals Statistics, 2(1): 134 – 147.
LEMMERLING, P. 1999. Structured Total Least Squares: Analysis, Algorithms and Applications. Doctorate Thesis, Katholieke Universiteit Leuven, Belgium.
LJUNG, L. 1979. Asymptotic behaviour of the extended Kalman filer as a parameter estimator for linear systems. IEEE Transactions on Automatic Control, AC-24: 36 – 50.
LJUNG, L. 1999. System Identification – Theory for the User. 2nd edition. Upper Saddle River, NJ: Prentice Hall.
LJUNG, L. and SÖDERSTRÖM, T. 1986. Theory and Practice of Recursive Identification. Cambridge, MA: MIT Press.
LORDEN, G. 1971. Procedures for reacting to a change in distribution. Annals Mathematical Statistics, 42: 1897 – 1908.
BIBLIOGRAPHY
173
MALLAT, S. 1989. A theory for multiresolution signal decomposition: the wavelet representation. IEEE Transactions on Pattern Analysis and Machine Intelligence, 11(7): 674 – 693.
MALLAT, S. 1998. A wavelet tour of signal processing. New York: Academic Press.
MATHWORKS INC. 2002a. MATLAB® Documentation- System Identification Toolbox. Version 6.5.0.180913a (R13). Nattick, MA.
MATHWORKS INC. 2002b. MATLAB® Documentation- Wavelet Toolbox. Version 6.5.0.180913a (R13). Nattick, MA.
MEYER, Y. 1993. Les ondelettes. Algorithmes et applications. Paris: Colin Ed., 2nd edition. (English translation: Wavelets: algorithms and applications, SIAM).
MICROSOFT CORPORATION. 1999. Windows® Terminal Services guide. Redmond, WA. Available: http://www.microsoft.com/windows2000/technologies/terminal/default.asp
MICROSOFT CORPORATION. 2000. Remote Desktop Protocol (RDP) Features and Performances. Redmond, WA. Available: http://www.microsoft.com/windows2000/techinfo/howitworks/rdpfandp.asp
OPPENHEIM, A.V. and SCHAFER, R.W. 1989. Discrete-Time Signal Processing. Englewood Cliffs, NJ: Prentice Hall.
174
BIBLIOGRAPHY
ORGANISATION FOR THE ADVANCEMENT OF STRUCTURED INFORMATION STANDARDS (OASIS) CONSORTIUM. Universal Description, Discovery and Integration (UDDI) resources. Billerica, MA. Available: http://www.uddi.org
PAGE, E.S. 1954. Continuous inspection schemes. Biometrika, 41: 100 – 115.
PILLAY, P. and BHATTACHARJEE, A. 1996. Application of Wavelets to model shortterm power system disturbances. IEEE Transactions on Power Systems, 11(4): 2031 – 2037.
PHILIPPOT, L. 1996. Parameter Estimation and Error Estimation for Line Fault Location and Distance Protection in Power Transmission Systems. PhD dissertation. Université Libre de Bruxelles, Belgium.
QIAN, S. 2002. Introduction to Time-Frequency and Wavelet Transforms, Upper Saddle River, NJ: Prentice Hall.
RABINER, L.R., and Gold, B. 1975. Theory and Application of Digital Signal Processing. Englewood Cliffs, NJ: Prentice Hall.
RIBEIRO, P.F. 1994. Wavelet transform: an advanced tool for analysing non-stationary harmonic distortions in power systems. In Proceedings of IEEE International Conference on Harmonics in Power Systems, Bologna, Italy.
ROBERTS, S.W. 1959. Control charts based on geometric moving averages. Technomet-
BIBLIOGRAPHY
175
rics, 1: 239 – 250.
ROBERTSON, D., CAMPS, O. I. and MAYER, J.S. 1994. Wavelets and power system transients: Feature detection and classification. In Proceedings of SPIE International Symposium on Optical Engineering in Aerospace Sensing, pp. 474 – 487.
ROBERTSON, D., CAMPS, O.I., MAYER, J.S. and GISH, W.B. 1996. Wavelets and electromagnetic power system transients. IEEE Transactions on Power Delivery, 11(2): 1050 – 1056.
RUIZ-VEGA, D., MESSINA, A.R. and ENRIQUEZ-HARPER, G. 2005. Analysis of Inter-area Oscillations via Non-linear Time Series Analysis Techniques. In Proceedings of Power System Computation Conference (PSCC’05), Belgium.
SANTOSO, S., POWERS, E.J. and GRADY, W.M. 1994. Electric power quality disturbance detection using wavelet transform analysis. In Proceedings of IEEE-SP International Symposium on Time-Frequency and Time-Scale Analysis, Philadelphia. pp. 166 – 169.
SANTOSO, S., POWERS, E.J., GRADY, W.M. and HOFFMAN, P. 1996. Power Quality Assessment via Wavelet Transform Analysis. IEEE Transactions on Power Delivery, 11(2): 924 – 930.
SHEWHART, W.A. 1931. Economic Control of Quality Manufactured Product. Princeton, NJ: D. Van Nostrand Reinhold.
176
BIBLIOGRAPHY
SIEMENS. 1999. System Program OSCOP P V6.1. Program manual.
SMOLA, A. and SCHÖLKOPF, B. 2002. Learning with Kernels. Cambridge, MA: MIT Press.
SOLO, V. 1980. Some aspects of recursive parameter estimation. International Journal of Control, 32: 395 – 410.
STOKES-WALLER, E. 1998. Automated Digital Fault Recording Analysis on the Eskom Transmission System. In Proceedings of Southern African Conference on Power System Protection, South Africa.
STOKES-WALLER, E. and KELLER, P. 1998. Power Network and Protection Performance Analysis on the Eskom Transmission System based on Digital Fault Records. In Proceedings of Southern African Conference on Power System Protection, South Africa.
STRANG, G. and NGUYEN, T. 1996. Wavelets and filter banks. Wellesley, MA: Wellesley-Cambridge Press.
SUN MICROSYSTEMS INC. Java™ Web Services white papers. Palo Alto, CA. Available: http://java.sun.com/Webservices
SUN MICROSYSTEMS INC. 2001. Building Web Services Forte™ for Java™ Programming Series. Palo Alto, CA.
BIBLIOGRAPHY
177
TARANTELLA INC. 2000. Tarantella® white papers and guide. Santa Cruz, CA. Available: http://www.tarantella.com
THINSOFT PTE LTD. 2002. WinConnect® Server XP™ manual and guide. Singapore. Available: http://www.thinsoftinc.com
THOTTAN, M. and JI, C. 2003. Anomaly Detection in IP Networks. IEEE Transactions on Signal Processing, 51(8): 2191 – 2204.
UKIL, A. and ZIVANOVIC, R. 2005a. Detection of Abrupt Changes in Power System Fault Analysis: A Comparative Study. In Proceedings of Southern African Universities Power Engineering Conference (SAUPEC’05), Johannesburg, South Africa.
UKIL, A. and ZIVANOVIC, R. 2005b. Abrupt Change Detection in Power System Fault Analysis using Wavelet Transform. In Proceedings of International Conference on Power Systems Transients (IPST’05), Montreal, Canada.
UKIL, A. and ZIVANOVIC, R. 2005c. Application of Abrupt Change Detection in Relay Performance Monitoring. In Proceedings of International Universities Power Engineering Conference (UPEC’05), Cork, Ireland.
UKIL, A. and ZIVANOVIC, R. 2005d. Application of Abrupt Change Detection in Power Systems Disturbance Analysis and Relay Performance Monitoring. IEEE Transactions on Power Delivery (Under review).
178
BIBLIOGRAPHY
UKIL, A. and ZIVANOVIC, R. 2005e. Abrupt Change Detection in Power System Fault Analysis using Adaptive Whitening Filter and Wavelet Transform. Electric Power Systems Research, Elsevier (In press).
UKIL, A. and ZIVANOVIC, R. 2005f. The Detection of Abrupt Changes using Recursive Identification for Power System Fault Analysis. Electric Power Systems Research, Elsevier (Under review).
UKIL, A. and ZIVANOVIC, R. 2005g. Adjusted Haar Wavelet for application in the Power Systems Disturbance Analysis. Digital Signal Processing, Elsevier (Under review).
UKIL, A. and ZIVANOVIC, R. 2005h. Practical Approaches in Application Service Provider: A Comparative Study. In Proceedings of Second World Congress on Lateral-Computing (WCLC 2005), Bangalore, India.
UKIL, A. and ZIVANOVIC, R. 2006. Automatic Signal Segmentation based on Abrupt Change Detection for Power Systems Applications. In Proceedings of 2006 IEEE Power India Conference, New Delhi, India (Accepted).
VAPNIK, V. 1995. The Nature of Statistical Learning Theory. New York: Springer.
WANG, Y. 1998. Change curve estimation via wavelets. Journal of American Statistical Association, 93: 163 – 172.
BIBLIOGRAPHY
179
WILKINSON, W. A. and COX, M. D. 1996. Discrete Wavelet Analysis of Power System Transients. IEEE Transactions on Power Systems, 11(4): 2038 – 2044.
WIOT, D. 2004. A New Adaptive Transient Monitoring Scheme for Detection of Power System Events. IEEE Transactions on Power Delivery, 19(1): 42 – 48.
WORLD WIDE WEB CONSORTIUM (W3C). 2000. Simple Object Access Protocol (SOAP) 1.1. Available: http://www.w3.org/TR/2000/NOTE-SOAP-20000508
WORLD WIDE WEB CONSORTIUM (W3C). 2001. Web Services Description Language (WSDL) 1.1. Available: http://www.w3.org/TR/wsdl
ZIVANOVIC, R. 1999. Frequency Estimation Algorithm based on Local Polynomial Approximation. In Proceedings of International Universities Power Engineering Conference (UPEC’99), Leicester, UK.
ZIVANOVIC, R. 2001. Nonparametric Frequency Estimation for Power System Applications. In Proceedings of IEEE Porto Power Tech’2001, Porto, Portugal.
180
APPENDIX A:COMTRADE
APPENDIX A
COMTRADE
The variety of sources of transient data, such as Digital Fault Recorders, Digital Protective Relays and Transient Simulation Programmes from different manufacturers using proprietary or different standard formats, made it necessary to introduce the IEEE Standard for
Common Format for Transient Data Exchange (Comtrade) for Power Systems as specified by the IEEE STANDARD C37.111-1991 (1991).
The standard files should be ASCII files. Each event should have three types of files associated with it. Each of the three types carries a different class of information: header (*.HDR), configuration (*.CFG) and data (*.DAT).
The intent of the header (*.HDR) is to provide supplementary information in a narrative form for the user to better understand the conditions of the transient record. The header file is not intended to be manipulated by an applications program. The following elements should be included in the header file:
•
Description of the power system prior to disturbance
•
Name of the station
•
Identity of the line, transformer, reactor, capacitor or circuit breaker that experienced the transient
APPENDIX A:COMTRADE
•
Length of the faulted line
•
Positive and zero sequence resistances and reactances
•
Capacitances
•
Mutual coupling between parallel lines
•
Locations and ratings of shunt reactors and series capacitors
•
Nominal voltage ratings of transformer windings
•
Transformer power ratings and windings connections
•
Positive and zero sequence impedance of the source
•
Description of how the data was obtained
•
Description of the anti-aliasing filters used
•
Description of analog mimic circuitry
•
Number of discs on which the case data is stored
•
The format in which the data is recorded
•
The headings of the columns of the data table.
181
The intent of the configuration file (*.CFG) is to provide the information necessary for a computer program to read and interpret the data values in the associated data files. Since configuration file is in a predefined, fixed format, a computer program does not have to be customised for each configuration file. The content consists of the following elements:
•
Station name and identification
182
APPENDIX A:COMTRADE
•
Number and type of channels
•
Channel names, units and conversion factors
•
Line frequency
•
Sample rate and number of samples at this rate
•
Date and time of first data value
•
Date and time of trigger point
•
File type.
The data file (*.DAT) contains the data values in rows and columns where each row consists of a set of data values preceded by a sequence number and the time for that set of data values. No other information is contained in the data file. The first column contains the sample number of the data set in that row. The second column gives the time of the data in microseconds from the beginning of the record. The third and remaining columns contain the data values that represent voltages, currents and status information.
APPENDIX B: MATLAB® Implementation Codes
183
APPENDIX B
MATLAB® Implementation Codes
In this section we provide the Matlab® implementation codes for the different algorithms. The codes are presented in the alphabetical order of the Matlab® filenames.
****************************** abhicreate.m ****************************** % Creates Abhisek Wavelet and adds to Wavelet toolbox % Written by Abhisek Ukil %filename contains wavelet filter name (m-files): abhiwavf.m, abhiwavelet.m filename='abhiwavelet'; %filename='abhifilter'; % Delete any previously registered Abhisek Wavelet wavemngr('del','Abhisek'); % Add any previously registered Abhisek Wavelet wavemngr('add','Abhisek','abhi',1,'',filename);
*************************************************************************
184
APPENDIX B: MATLAB® Implementation Codes
****************************** abhiwavelet.m ***************************** % Creates Abhisek Mother-Wavelet using Quadrature Mirror Filter % Based on FIR filter % To see wavelet functions: wavefun('abhi',0,'plot'); % Written by Abhisek Ukil function F=abhiwavelet(wname); %F=[0.5 0.5]; %haar %F=[0.5 0 0 0.5]; %Adjustment parameter=1 F=[0.5 0 0 0 0 0.5]; %Adjustment parameter=2
*************************************************************************
APPENDIX B: MATLAB® Implementation Codes
185
******************************** abrupt.m *******************************
% Example Abrupt Change Detection % Written by Abhisek Ukil t=0:pi/50:2*pi; xa=5*cos(50.*t); xb=2*cos(50.*t); a=ones(101,1); a(20:40)=0; a(60:80)=0; a=a'; b=zeros(101,1); b(20:40)=1; b(60:80)=1; b=b'; x=xa.*a + xb.*b ; x=x'; %seg = segment([x,ones(size(x))],[0 1 0],0.5); %Using System %Identification Toolbox seg = segment([x,ones(size(x))],[1 0 1],0.1); subplot(2,1,1) plot(x) title('Signal with Abrupt changes') subplot(2,1,2) %plot([seg x]) plot(seg) %title('Time-Instants of Abrupt changes (Up or Down Spike shows sudden Increase or Decrease of Value)') title('Time-Instants of Abrupt changes')
*************************************************************************
186
APPENDIX B: MATLAB® Implementation Codes
***************************** abrupt_adamo.m **************************** % Abrupt Change Detection Script using Adaptive Whitening Filter and % The wavelet transform % Written by Abhisek Ukil
% Sampling Freqency=2.5 KHz, Fundamental Frequency=50 Hz Fs=2500; Ts=1/Fs; Ffund=50; C=round(Fs/Ffund); % Pulsation Frequency Fpulse=51; W0=2*pi*Fpulse; % Alpha and Beta coefficients alpha=sin(W0*Ts*C)/sin(W0*Ts); beta=cos(W0*Ts*C)-alpha*cos(W0*Ts); % Creating the FIR filter % B(z)=1-alpha*z^(-C+1)-beta*z^(-C) b=zeros(1,C); b(1)=1; b(C-1)=-alpha; b(C)=-beta; a=[1]; % for FIR filter no a % Frequency Response of the FIR filter [db,mag,pha,grd,w] = freqz_m(b,a); % Plotting Frequency Responses figure(1) subplot(2,2,1); plot(db);grid; title('Log-Magnitude Response'), xlabel('Frequency in Hz'), ylabel('LogMagnitude, Decibels'), subplot(2,2,2); plot(mag); grid; title('Magnitude Response'),xlabel('Frequency in Hz'),ylabel('Magnitude') subplot(2,2,3); plot(pha);grid; title('Phase Response'),xlabel('Frequency in Hz'),ylabel('Phase') subplot(2,2,4); plot(grd),title('Group Delay'),xlabel('Frequency in Hz'),ylabel('Group Delay'), %Signal Filtering
APPENDIX B: MATLAB® Implementation Codes
187
[name,path,flag]=uigetfile({'*.dat','Comtrade DATA-files *.dat)'},'Select Input Comtrade DATA File'); % flag=1 means DAT File is Loaded, If flag=0 then Cancel is clicked i.e. %Null file if flag ~=0 filename=strcat(path,name); % Make Full Path String for Input file file_a=load(filename); datacolumn=input('Data Column : '); in=file_a(:,datacolumn); % Normalise signal norm_in=in./mean(abs(in)); % Filter the Input with ADAMO output=filter(b,a,norm_in); % Abrupt Change Detection using Wavelet Analysis ti=func_wavelet(output); %ti=func_mywavelet(output,'abhi'); %Plotting output figure(2) subplot(3,1,1); plot(in) title('Original Fault signal') subplot(3,1,2); plot(output) title('ADAMO filtered signal') subplot(3,1,3); plot(ti) title('Time-instants of Abrupt Changes') end
*************************************************************************
188
APPENDIX B: MATLAB® Implementation Codes
*************************** abrupt_mywavelet.m *************************** % Abrupt Change Detection Script using Adjusted Wavelet (abhi) method % Written by Abhisek Ukil warning off MATLAB:divideByZero % Dialog box to input file [name,path,flag]=uigetfile({'*.dat','Comtrade DATA-files (*.dat)'},'Select Input Comtrade DATA File'); % flag=1 means DAT File is Loaded, If flag=0 then Cancel is clicked i.e. %Null file if flag ~=0 filename=strcat(path,name); % Make Full Path String for Input file a=load(filename); datacolumn=input('Data Column : '); in=a(:,datacolumn); % Normalise signal norm_in=in./mean(abs(in)); % Wavelet Analysis [c1,c2]=dwt(norm_in,'abhi'); % Filtering Response %Change the negatives into positives c2(:)=abs(c2); for i=1:size(c2) if c2(i)<0 c2(i)=-c2(i); end end % If whole signal is same i.e. unchanged then some low value of c2 % overall for i=1:size(c2) if c2(i)<0.05 c2(i)=0; end end
APPENDIX B: MATLAB® Implementation Codes
189
%Calculating Threshold sd=std(c2); % Standard Deviation of c2 %sd=func_DJmd(c2); %Median Absolute deviation threshold1=sd*sqrt(2*log(length(c2)));
% Length after wavelet filter
a=zeros(size(c2)); a(1:size(a))=threshold1; %Filtering x=c2; for i=1:size(x) if x(i)
% +++++++ Refinement based on number of segments +++++++++ timeinstant=zeros(size(x)); segment=4; if length(passed2)>segment threshold2=std(m)*sqrt(2*log(length(m))); b=zeros(size(x));
190
APPENDIX B: MATLAB® Implementation Codes
b(1:size(b))=threshold2; timeinstant=x>threshold2; %refinegroup=func_compressband(timeinstant,x,20,20); refinegroup=func_compressband(timeinstant,x,40,20); timeinstant(:)=0; if length(refinegroup)<=segment timeinstant(refinegroup)=1; end
else %timeinstant(:)=0; timeinstant(passed2)=1; end % ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ else timeinstant=x; end %Removing Extreme End and Front glitches end_distance=500; for i=1:length(timeinstant) if i>(length(timeinstant)-end_distance) | i<(end_distance/5) if timeinstant(i)==1 timeinstant(i)=0; end end end %Adjust the Half-sampled timeinstant [ti_i,ti_j]=find(timeinstant); timeinstant=zeros(size(in)); if ti_i for k=1:length(ti_i) timeinstant(ti_i(k)*2)=1; end end
APPENDIX B: MATLAB® Implementation Codes
191
% Plotting figure(1) subplot(3,1,1) plot(in) title('Fault Signal with Abrupt Changes') subplot(3,1,2) %plot([c2,a]) plot(c2,'b'); hold; plot(a,'k:') title('Wavelet Coefficients (blue), Threshold (dashed, black)') subplot(3,1,3) plot(timeinstant) title('Change Time-instants') end
*************************************************************************
192
APPENDIX B: MATLAB® Implementation Codes
**************************** abrupt_recursive.m *************************** % Abrupt Change Detection Script using Recursive Identification method % Written by Abhisek Ukil % Dialog box to input file [name,path,flag]=uigetfile({'*.dat','Comtrade DATA-files (*.dat)'},'Select Input Comtrade DATA File'); % flag=1 means DAT File is Loaded, If flag=0 then Cancel is clicked i.e. %Null file if flag ~=0 filename=strcat(path,name); % Make Full Path String for Input file a=load(filename); datacolumn=input('Data Column : '); in=a(:,datacolumn); % Normalise signal norm_in=in./mean(abs(in)); % Segmentation [seg,v,thm,r2e] = segment([norm_in,ones(size(norm_in))],[4 4 1],0.001); %Using System Identification Toolbox segm=seg(:,8); %segm=thm(:,8); %Changing -ve into +ve segm(:)=abs(segm); %Filtering threshold=std(segm)*sqrt(2*log(length(segm))); timeinstant=segm>threshold; a=zeros(size(segm)); a(1:size(a))=threshold; group=func_compressband(timeinstant,segm,60,48); %group=func_compressband(timeinstant,segm,40,20); timeinstant(:)=0; timeinstant(group)=1; %Removing End glitches
APPENDIX B: MATLAB® Implementation Codes
193
end_distance=1500; for i=1:length(timeinstant) if i>(length(timeinstant)-end_distance) if timeinstant(i)==1 timeinstant(i)=0; end end end %Removing if number of segments more than a specified number segmentlimit=8; if length(find(timeinstant))>segmentlimit timeinstant(:)=0; end % Plotting figure(1) subplot(3,1,1) plot(in) title('Signal with Abrupt Changes') subplot(3,1,2) plot(seg(:,2)) title('System Modelling'),ylabel('Scaled value') %axis([0,7000,-3,3]) subplot(3,1,3) plot(timeinstant) title('Time-Instants of Abrupt Changes'),ylabel('Impulse indicator')
end
*************************************************************************
194
APPENDIX B: MATLAB® Implementation Codes
**************************** abrupt_segment.m *************************** % Abrupt Change Detection Script using Complete Algorithm % Written by Abhisek Ukil %clear warning off MATLAB:divideByZero % Dialog box to input file [name,path,flag]=uigetfile({'*.dat','Comtrade DATA-files (*.dat)'},'Select Input Comtrade DATA File'); if flag ~=0 % flag=1 means DAT File is Loaded, If flag=0 then Cancel is clicked i.e. Null file filename=strcat(path,name); % Make Full Path String for Input file a=load(filename); datacolumn=input('Data Column : '); in=a(:,datacolumn); % Abrupt Change Detection using Wavelet Analysis timeinstant=func_segment(in);
figure(1) plot(in), title(['Original Fault Signal, File: ' name ', Data Column: ' num2str(datacolumn) ', Segmentation (dashed)']), hold on, plot(timeinstant.*(1.5*max(abs(in))),'k:'), hold on, plot(timeinstant.*(1.5*max(abs(in))),'k:')
end
*************************************************************************
APPENDIX B: MATLAB® Implementation Codes
195
**************************** abrupt_wavelet.m *************************** % Abrupt Change Detection Script using Wavelet method % Written by Abhisek Ukil warning off MATLAB:divideByZero % Dialog box to input file [name,path,flag]=uigetfile({'*.dat','Comtrade DATA-files (*.dat)'},'Select Input Comtrade DATA File'); % flag=1 means DAT File is Loaded, If flag=0 then Cancel is clicked i.e. %Null file if flag ~=0 filename=strcat(path,name); % Make Full Path String for Input file a=load(filename); datacolumn=input('Data Column : '); in=a(:,datacolumn); % Normalise signal norm_in=in./mean(abs(in)); % Wavelet Analysis % wavelets: db2, db3, db4, db5, db6,.. sym1, sym2, sym3.., % coif1,coif2,...,coif5 % Best: sym1, sym2, db1, db4 [c1,c2]=dwt(norm_in,'db4'); % Filtering Response %Change the negatives into positioves c2(:)=abs(c2); for i=1:size(c2) if c2(i)<0 c2(i)=-c2(i); end end % If whole signal is same i.e. unchanged then some low value of c2 % overall for i=1:size(c2)
196
APPENDIX B: MATLAB® Implementation Codes
if c2(i)<0.05 c2(i)=0; end end %Calculating Threshold sd=std(c2); % Standard Deviation of c2 %sd=func_DJmd(c2); %Median Absolute deviation threshold1=sd*sqrt(2*log(length(c2)));
% Length after wavelet filter
a=zeros(size(c2)); a(1:size(a))=threshold1; %Filtering x=c2; for i=1:size(x) if x(i)
% ++++++++++ Refinement based on number of segments ++++++++++++ timeinstant=zeros(size(x)); segment=4;
APPENDIX B: MATLAB® Implementation Codes
197
if length(passed2)>segment threshold2=std(m)*sqrt(2*log(length(m))); b=zeros(size(x)); b(1:size(b))=threshold2; timeinstant=x>threshold2; %refinegroup=func_compressband(timeinstant,x,20,20); refinegroup=func_compressband(timeinstant,x,40,20); timeinstant(:)=0; if length(refinegroup)<=segment timeinstant(refinegroup)=1; end
else %timeinstant(:)=0; timeinstant(passed2)=1; end %++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ else timeinstant=x; end %Removing Extreme End and Front glitches end_distance=500; for i=1:length(timeinstant) if i>(length(timeinstant)-end_distance) | i<(end_distance/5) if timeinstant(i)==1 timeinstant(i)=0; end end end %Adjust the Half-sampled timeinstant [ti_i,ti_j]=find(timeinstant); timeinstant=zeros(size(in)); if ti_i for k=1:length(ti_i) timeinstant(ti_i(k)*2)=1; end
198
APPENDIX B: MATLAB® Implementation Codes
end
% Plotting figure(1) subplot(3,1,1) plot(in) title('Fault Signal with Abrupt Changes') subplot(3,1,2) %plot([c2,a]) plot(c2,'b'); hold; plot(a,'k:') title('Wavelet Coefficients (blue), Threshold (dashed, black)') subplot(3,1,3) plot(timeinstant) title('Change Time-instants') end
*************************************************************************
APPENDIX B: MATLAB® Implementation Codes
199
***************************** func_adamo.m ****************************** % Developes the ADAMO Filter and returns the FIR filter coefficients % Written by Abhisek Ukil function [b,a]=func_adamo() % Sampling Freqency=2.5 KHz, Fundamental Frequency=50 Hz Fs=2500; Ts=1/Fs; Ffund=50; C=round(Fs/Ffund); % Pulsation Frequency Fpulse=50.1; W0=2*pi*Fpulse; % Alpha and Beta coefficients alpha=sin(W0*Ts*C)/sin(W0*Ts); beta=cos(W0*Ts*C)-alpha*cos(W0*Ts); % Creating the FIR filter % B(z)=1-alpha*z^(-C+1)-beta*z^(-C) b=zeros(1,C); b(1)=1; b(C-1)=-alpha; b(C)=-beta; a=[1]; % for FIR filter no a % Frequency Response of the FIR filter [db,mag,pha,grd,w] = freqz_m(b,a); % Plotting Frequency Responses of the ADAMO filter figure(1) subplot(2,2,1); plot(db);grid; title('Log-Magnitude Response'), xlabel('Frequency in Hz'), ylabel('LogMagnitude, Decibels'), subplot(2,2,2); plot(mag); grid; title('Magnitude Response'),xlabel('Frequency in Hz'),ylabel ('Magnitude') subplot(2,2,3); plot(pha);grid; title('Phase Response'),xlabel('Frequency in Hz'),ylabel('Phase') subplot(2,2,4); plot(grd),title('Group Delay'),xlabel('Frequency in Hz'),ylabel('Group Delay')
*************************************************************************
200
APPENDIX B: MATLAB® Implementation Codes
*************************** func_compressband.m ************************* % Smoothing Filtering: compresses multiple lines band into one line. % Written by Abhisek Ukil function group= func_compressband(test,master,bandlength,grouplength) %bandlength=20; passed=find(test); count=1; buffer(count)=passed(1); bufferval=master(buffer); groupcount=0; for i=2:length(passed) if (passed(i)-passed(i-1))
APPENDIX B: MATLAB® Implementation Codes
201
end % This section is for end segment %buffer=buffer groupcount=groupcount+1; bufferval=master(buffer); if count
*************************************************************************
202
APPENDIX B: MATLAB® Implementation Codes
******************************* func_DJmd.m ****************************
% Calculates Median Absolute Deviation % Written by Abhisek Ukil function md=func_DJmd(x) n=length(x); sq=x.^2; sumsq=sum(sq); med=median(x); %md=(sqrt(sumsq/n-med^2))/0.6725; md=sum(abs(x-med)); md=sqrt(md);
*************************************************************************
APPENDIX B: MATLAB® Implementation Codes
203
*************************** func_mywavelet.m **************************** % Performs the signal decomposition of the given signal using Wavelet % Transform, followed by Threshold checking based on 'Universal % Threshold' technique of Donoho and Johnstone. % Inputs: signal to be processed, mother wavelet name. % Returns the abrupt change time-instant(s) vector. % Written by Abhisek Ukil function timeinstant=func_mywavelet(in,wavename) % Wavelet Analysis % Best: abhi [c1,c2]=dwt(in,wavename);
%abhi
% Filtering Response %Change the negatives into positioves c2(:)=abs(c2); % If whole signal is same then some low value of c2 overall for i=1:size(c2) if c2(i)<0.05 c2(i)=0; end end %Calculating Threshold sd=std(c2); % Standard Deviation of c2 %sd=func_DJmd(c2); % Donoho-Johnstone MAD threshold1=sd*sqrt(2*log(length(c2))); a=zeros(size(c2)); a(1:size(a))=threshold1; %Filtering x=c2; for i=1:size(x) if x(i)
% Length after wavelet filter
204
APPENDIX B: MATLAB® Implementation Codes
end %passed=find(x); if (find(x)) %group=func_compressband(x,x,20,20); group=func_compressband(x,x,40,20); y=zeros(size(x)); for i=1:length(group) y(group(i))=x(group(i)); end passed2=find(y); for i=1:size(passed2) m(i)=x(passed2(i)); end
% +++++++++++ Refinement based on number of segments ++++++++++++ timeinstant=zeros(size(x)); segment=4; if length(passed2)>segment threshold2=std(m)*sqrt(2*log(length(m))); b=zeros(size(x)); b(1:size(b))=threshold2; timeinstant=x>threshold2; %refinegroup=func_compressband(timeinstant,x,20,20); refinegroup=func_compressband(timeinstant,x,40,20); timeinstant(:)=0; if length(refinegroup)<=segment timeinstant(refinegroup)=1; end else %timeinstant(:)=0; timeinstant(passed2)=1; end % ++++++++++++++++++++++++++++++++++++++++++++++++++++++
APPENDIX B: MATLAB® Implementation Codes
205
else timeinstant=x; end % Removing Extreme End and Front glitches end_distance=500; for i=1:length(timeinstant) if i>(length(timeinstant)-end_distance) | i<(end_distance/5) if timeinstant(i)==1 timeinstant(i)=0; end end end
%Adjust the Half-sampled timeinstant [ti_i,ti_j]=find(timeinstant); timeinstant=zeros(size(in)); if ti_i for k=1:length(ti_i) timeinstant(ti_i(k)*2)=1; end end
%
% Plotting
%
figure(1)
%
subplot(3,1,1)
%
plot(in)
%
title('Fault Signal with Abrupt Changes')
%
subplot(3,1,2)
%
%plot([c2,a])
%
plot(c2,'b'); hold; plot(a,'k:')
%
title('Wavelet Coefficients (blue), Threshold (dashed, black)')
%
subplot(3,1,3)
%
plot(timeinstant)
%
title('Change Time-instants')
*************************************************************************
206
APPENDIX B: MATLAB® Implementation Codes
***************************** func_recursive.m **************************** % Performs the abrupt change detection based on recursive identification. % Input: signal to be processed. % Returns the abrupt change time-instant(s) vector. % Written by Abhisek Ukil function timeinstant=func_recursive(norm_in) % Segmentation [seg,v,thm,r2e] = segment([norm_in,ones(size(norm_in))],[4 4 1],0.001); %Using System Identification Toolbox segm=seg(:,8); %Changing -ve into +ve segm(:)=abs(segm); %Filtering threshold=std(segm)*sqrt(2*log(length(segm))); timeinstant=segm>threshold; a=zeros(size(segm)); a(1:size(a))=threshold; group=func_compressband(timeinstant,segm,40,20); timeinstant(:)=0; timeinstant(group)=1; %Removing End glitches end_distance=1500; for i=1:length(timeinstant) if i>(length(timeinstant)-end_distance) if timeinstant(i)==1 timeinstant(i)=0; end end end %Removing if number of segments more than a specified number segmentlimit=6; if length(find(timeinstant))>segmentlimit timeinstant(:)=0;
APPENDIX B: MATLAB® Implementation Codes
207
end %
figure(1)
%
subplot(2,1,1)
%
plot(in)
%
title('Signal with Abrupt Changes')
%
% subplot(3,1,2)
%
% plot([segm,a])
%
subplot(2,1,2)
%
plot(timeinstant)
%
title('Time-Instants of Abrupt Changes')
*************************************************************************
208
APPENDIX B: MATLAB® Implementation Codes
***************************** func_segment.m **************************** % Performs Automatic Segmentation based on Abrupt Change Detection. % % Takes as input the signal to be segmented as vector, % Returns the change time-instants indicated as unit impulses as vector, % if there is no change, returns null. % % Format: timeinstant=func_segment(in) % % The return vector is of same length as the input signal, with all elements being zero, % except at those time-instants where the input signal has abrupt changes, if any. % These positions in the return vector are marked with unit impulses i.e., 1. % So, the index of the non-zero elements of the return vector indicate the % change time-instants of the original signal. If there is no abrupt change % in the input signal, all the elements of the return vector will be zero. % % Written by Abhisek Ukil function timeinstant=func_segment(in) % Perform Normalisation of the Input Signal norm_in=in./mean(abs(in)); % Starting Algorithm 'abhi' timeinstant=func_mywavelet(norm_in,'abhi'); if length(find(timeinstant))<3 switch length(find(timeinstant)) case 0 timeinstant=func_wavelet(norm_in); case 1 timeinstant=func_wavelet(norm_in); len=length(find(timeinstant)); if len>2 % Load the ADAMO filter
APPENDIX B: MATLAB® Implementation Codes
[ADAMO_b,ADAMO_a]=func_adamo; % Filter the Input with ADAMO output=filter(ADAMO_b,ADAMO_a,norm_in); % Abrupt Change Detection using Wavelet Analysis timeinstant=func_wavelet(output); elseif len<2 %timeinstant=timeinstant; timeinstant=func_mywavelet(norm_in,'abhi'); else len==2 timeinstant=func_recursive(norm_in); if length(find(timeinstant))>=2 timeinstant=func_wavelet(norm_in); end end case 2 timeinstant=func_wavelet(norm_in); switch length(find(timeinstant)) case 0 timeinstant=timeinstant; case 1 timeinstant=func_recursive(norm_in); if length(find(timeinstant))>2 timeinstant=func_mywavelet(norm_in,'abhi'); end case 2 timeinstant=func_mywavelet(norm_in,'abhi'); case 3 % Load the ADAMO filter [ADAMO_b,ADAMO_a]=func_adamo; % Filter the Input with ADAMO output=filter(ADAMO_b,ADAMO_a,norm_in); % Abrupt Change Detection using Wavelet Analysis timeinstant=func_wavelet(output); otherwise
209
210
APPENDIX B: MATLAB® Implementation Codes
timeinstant=func_recursive(norm_in); if length(find(timeinstant))>3 timeinstant=func_mywavelet(norm_in,'abhi'); end end end end
*************************************************************************
APPENDIX B: MATLAB® Implementation Codes
211
***************************** func_wavelet.m ***************************** % Performs the signal decomposition of the given signal using Wavelet % Transform, followed by Threshold checking based on 'Universal % Threshold' technique of Donoho and Johnstone. % Inputs: signal to be processed, mother wavelet name. % Returns the abrupt change time-instant(s) vector. % Written by Abhisek Ukil function timeinstant=func_wavelet(in,wavename) % wavelets: db2, db3, db4, db5, db6,.. sym1, sym2, sym3.., % coif1,coif2,...,coif5 % Best: sym1, sym2, db1, db4, abhi [c1,c2]=dwt(in,'db4');
%db4
% Filtering Response %Change the negatives into positioves c2(:)=abs(c2); % If whole signal is same then some low value of c2 overall for i=1:size(c2) if c2(i)<0.05 c2(i)=0; end end %Calculating Threshold sd=std(c2); % Standard Deviation of c2 %sd=func_DJmd(c2); % Donoho-Johnstone MAD threshold1=sd*sqrt(2*log(length(c2))); a=zeros(size(c2)); a(1:size(a))=threshold1; %Filtering x=c2; for i=1:size(x) if x(i)
% Length after wavelet filter
212
APPENDIX B: MATLAB® Implementation Codes
end %passed=find(x); if (find(x)) %group=func_compressband(x,x,20,20); group=func_compressband(x,x,40,20); y=zeros(size(x)); for i=1:length(group) y(group(i))=x(group(i)); end passed2=find(y); for i=1:size(passed2) m(i)=x(passed2(i)); end
% +++++++++++ Refinement based on number of segments ++++++++++++ timeinstant=zeros(size(x)); segment=4; if length(passed2)>segment threshold2=std(m)*sqrt(2*log(length(m))); b=zeros(size(x)); b(1:size(b))=threshold2; timeinstant=x>threshold2; %refinegroup=func_compressband(timeinstant,x,20,20); refinegroup=func_compressband(timeinstant,x,40,20); timeinstant(:)=0; if length(refinegroup)<=segment timeinstant(refinegroup)=1; end else %timeinstant(:)=0; timeinstant(passed2)=1; end % +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
APPENDIX B: MATLAB® Implementation Codes
213
else timeinstant=x; end % Removing Extreme End and Front glitches end_distance=500; for i=1:length(timeinstant) if i>(length(timeinstant)-end_distance) | i<(end_distance/5) if timeinstant(i)==1 timeinstant(i)=0; end end end
%Adjust the Half-sampled timeinstant [ti_i,ti_j]=find(timeinstant); timeinstant=zeros(size(in)); if ti_i for k=1:length(ti_i) timeinstant(ti_i(k)*2)=1; end end
%
% Plotting
%
figure(1)
%
subplot(3,1,1)
%
plot(in)
%
title('Fault Signal with Abrupt Changes')
%
subplot(3,1,2)
%
%plot([c2,a])
%
plot(c2,'b'); hold; plot(a,'k:')
%
title('Wavelet Coefficients (blue), Threshold (dashed, black)')
%
subplot(3,1,3)
%
plot(timeinstant)
%
title('Change Time-instants')
*************************************************************************
214
APPENDIX B: MATLAB® Implementation Codes
****************************** test_adamo.m ***************************** % Script to test ADAMO based Algorithm quickly % Written By: Abhisek Ukil clear % Load the ADAMO filter [ADAMO_b,ADAMO_a]=func_adamo; % Specify the start string startstr='s'; while lower(startstr)=='s' % Load the Signal [name,path,flag]=uigetfile({'*.dat','Comtrade DATA-files (*.dat)'},'Select Input Comtrade DATA File'); if flag ~=0 % flag=1 means DAT File is Loaded, If flag=0 then Cancel is clicked i.e. Null file
filename=strcat(path,name); % Make Full Path String for Input file file_a=load(filename); for i=3:10 % Iterate through Data Column 3-10 in=file_a(:,i); % Normalise signal norm_in=in./mean(abs(in)); %++++++++++++++++++++ Signal processing +++++++++++++++++++++ % Filter the Input with ADAMO output=filter(ADAMO_b,ADAMO_a,norm_in); % Abrupt Change Detection using Wavelet Analysis timeinstant=func_wavelet(output); %timeinstant=func_mywavelet(output,'abhi'); %Plotting output
APPENDIX B: MATLAB® Implementation Codes
215
figure(1) subplot(3,1,1); plot(in) title(['Original Fault Signal, File: ' name ', Data Column: ' num2str(i)]) subplot(3,1,2); plot(output) title('ADAMO filtered signal') subplot(3,1,3); plot(timeinstant) title('Time-instants of Abrupt Changes (Method: ADAMO Filter + Wavelet)')
%++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
% Look for next data column of same file if i<10 pause close all end end % Continue for loading next file continuestr=input('Want to continue? Enter "c" : ','s'); if lower(continuestr)=='c' startstr='s'; close all; else startstr=''; end else break %exit the main while loop of file loading end end
*************************************************************************
216
APPENDIX B: MATLAB® Implementation Codes
**************************** test_mywavelet.m **************************** % Script to test Adjusted Wavelet algorithm quickly % Written By: Abhisek Ukil clear %abcreate
%creates Abhisek wavelet if not exists
% Specify the start string startstr='s'; while lower(startstr)=='s' % Load the Signal [name,path,flag]=uigetfile({'*.dat','Comtrade DATA-files (*.dat)'},'Select Input Comtrade DATA File'); if flag ~=0 % flag=1 means DAT File is Loaded, If flag=0 then Cancel %is clicked i.e. Null file filename=strcat(path,name); % Make Full Path String for Input file file_a=load(filename); for i=3:10 % Iterate through Data Column 3-10 %datacolumn=input('Data Column : '); %in=file_a(:,datacolumn); in=file_a(:,i); % Normalise signal norm_in=in./mean(abs(in)); %++++++++++++++++++ Signal processing +++++++++++++++++++++
% Abrupt Change Detection using Wavelet Analysis timeinstant=func_mywavelet(norm_in,'abhi'); %db4 %Plotting output %
figure(1)
%
subplot(2,1,1); plot(in)
% title(['Original Fault Signal, File: ' name ', Data Column: ' num2str(i)])
APPENDIX B: MATLAB® Implementation Codes
217
%
subplot(2,1,2); plot(ti)
%
title('Time-instants of Abrupt Changes (Method: Wavelet)') figure(1)
plot(in), title(['Original Fault Signal, File: ' name ', Data Column: ' num2str(i) ', Segmentation (dashed)']), hold on, plot(timeinstant.*(1.5*max(abs(in))),'k:'), hold on, plot(timeinstant.*(1.5*max(abs(in))),'k:')
%++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
% Look for next data column of same file if i<10 pause close all end end % Continue for loading next file continuestr=input('Want to continue? Enter "c" : ','s'); if lower(continuestr)=='c' startstr='s'; close all; else startstr=''; end else break %exit the main while loop of file loading end end
*************************************************************************
218
APPENDIX B: MATLAB® Implementation Codes
***************************** test_recursive.m **************************** % Script to test Recursive Identification Algorithm quickly % Written By: Abhisek Ukil clear % Specify the start string startstr='s'; while lower(startstr)=='s' % Load the Signal [name,path,flag]=uigetfile({'*.dat','Comtrade DATA-files (*.dat)'},'Select Input Comtrade DATA File'); if flag ~=0 % flag=1 means DAT File is Loaded, If flag=0 then Cancel is clicked i.e. Null file for i=3:10 % Iterate through Data Column 3-10 filename=strcat(path,name); % Make Full Path String for Input file file_a=load(filename); %datacolumn=input('Data Column : '); %in=file_a(:,datacolumn); in=file_a(:,i); % Normalise signal norm_in=in./median(abs(in)); %++++++++++++++++++ Signal processing +++++++++++++++++++++++ % Abrupt Change Detection using Recursive Identification ti=func_recursive(norm_in); %Plotting output figure(1) subplot(2,1,1); plot(in) title(['Original Fault Signal, File: ' name ', Data Column: ' num2str(i)]) subplot(2,1,2); plot(ti) title('Time-instants of Abrupt Changes (Method: Recursive Identification)')
APPENDIX B: MATLAB® Implementation Codes
219
%++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ % Look for next data column of same file if i<10 pause close all end end % Continue for loading next file continuestr=input('Want to continue? Enter "c" : ','s'); if lower(continuestr)=='c' startstr='s'; close all; else startstr=''; end else break %exit the main while loop of file loading end end
*************************************************************************
220
APPENDIX B: MATLAB® Implementation Codes
****************************** test_segment.m **************************** % Script to test Complete Segmentation Algorithm quickly % Written By: Abhisek Ukil clear % Specify the start string startstr='s'; while lower(startstr)=='s' % Load the Signal [name,path,flag]=uigetfile({'*.dat','Comtrade (*.dat)'},'Select Input Comtrade DATA File');
DATA-files
if flag ~=0 % flag=1 means DAT File is Loaded, If flag=0 then Cancel is clicked i.e. Null file filename=strcat(path,name); % Make Full Path String for Input file file_a=load(filename); for i=3:10 % Iterate through Data Column 3-10 %datacolumn=input('Data Column : '); %in=file_a(:,datacolumn); in=file_a(:,i);
%+++++++++++++++++++++Segmentation ++++++++++++++++++++++++++
% Abrupt Change Detection using Wavelet Analysis timeinstant=func_segment(in); figure(1) plot(in), title(['Original Fault Signal, File: ' name ', Data Column: ' num2str(i) ', Segmentation (dashed)']), hold on, plot(timeinstant.*(1.5*max(abs(in))),'k:'), hold on, plot(timeinstant.*(1.5*max(abs(in))),'k:') %++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
APPENDIX B: MATLAB® Implementation Codes
221
% Look for next data column of same file if i<10 pause close all end end % Continue for loading next file continuestr=input('Want to continue? Enter "c" : ','s'); if lower(continuestr)=='c' startstr='s'; close all; else startstr=''; end else break %exit the main while loop of file loading end end
*************************************************************************
222
APPENDIX B: MATLAB® Implementation Codes
***************************** test_wavelet.m **************************** % Script to test Wavelet algorithm quickly % Written By: Abhisek Ukil clear % Specify the start string startstr='s'; while lower(startstr)=='s' % Load the Signal [name,path,flag]=uigetfile({'*.dat','Comtrade DATA-files (*.dat)'},'Select Input Comtrade DATA File'); if flag ~=0 % flag=1 means DAT File is Loaded, If flag=0 then Cancel is clicked i.e. Null file filename=strcat(path,name); % Make Full Path String for Input file file_a=load(filename); for i=3:10 % Iterate through Data Column 3-10 %datacolumn=input('Data Column : '); %in=file_a(:,datacolumn); in=file_a(:,i); % Normalise signal norm_in=in./mean(abs(in)); %+++++++++++++++++++ Signal processing ++++++++++++++++++++++ % Abrupt Change Detection using Wavelet Analysis timeinstant=func_wavelet(norm_in); %timeinstant=func_mywavelet(norm_in,'db4'); %Plotting output %
figure(1)
%
subplot(2,1,1); plot(in)
% title(['Original Fault Signal, File: ' name ', Data Column: ' num2str(i)]) %
subplot(2,1,2); plot(ti)
APPENDIX B: MATLAB® Implementation Codes
%
223
title('Time-instants of Abrupt Changes (Method: Wavelet)') figure(1)
plot(in), title(['Original Fault Signal, File: ' name ', Data Column: ' num2str(i) ', Segmentation (dashed)']), hold on, plot(timeinstant.*(1.5*max(abs(in))),'k:'), hold on, plot(timeinstant.*(1.5*max(abs(in))),'k:')
%++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
% Look for next data column of same file if i<10 pause close all end end % Continue for loading next file continuestr=input('Want to continue? Enter "c" : ','s'); if lower(continuestr)=='c' startstr='s'; close all; else startstr=''; end else break %exit the main while loop of file loading end end
*************************************************************************
