Real-time Applications ofSynchrophasor-based Wide-areaMeasurement SystemsbyMatin RahmatianM.Sc., Amirkabir University of Technology (Tehran Polytechnic), 2013B.Sc., Ferdowsi University of Mashhad, 2010a thesis submitted in partial fulfillmentof the requirements for the degree ofDoctor of Philosophyinthe faculty of graduate and postdoctoral studies(Electrical and Computer Engineering )The University of British Columbia(Vancouver)March 2018c© Matin Rahmatian, 2018The following individuals certify that they have read, and recommend to the Facultyof Graduate and Postdoctoral Studies for acceptance, the dissertation entitled:Real-time Applications of Synchrophasor-based Wide-area Measurement Systemssubmitted by Matin Rahmatian in partial fulfillment of the requirements forthe degree of Doctor of Philosophyin Electrical and Computer EngineeringExamining Committee:William G. Dunford, Electrical and Computer EngineeringCo-supervisorYu Christine Chen, Electrical and Computer EngineeringCo-supervisorJose´ Mart´ı, Electrical and Computer EngineeringSupervisory Committee MemberJohn Madden, Electrical and Computer EngineeringUniversity ExaminerRyozo Nagamune, Mechanical EngineeringUniversity ExaminerAdditional Supervisory Committee Members:Supervisory Committee MemberSupervisory Committee MemberiiAbstractIn this thesis, we focus on two major real-time applications1 of modern synchrophasor-based wide-area measurement systems (WAMS), i.e., transient stability assessment (TSA)and fault detection and identification (FDI). First, we develop a tool for real-time TSAbased on automatic learning (AL) approaches. We use Classification and Regression Tree(CART) as the classification tool and Multivariate Adaptive Regression Splines (MARS) asthe regression tool. To train and validate these tools in a practical setting, we conducttest cases on the full Western Electricity Coordinated Council (WECC) system model,with emphasis on the BC Hydro (BCH) power system. While being mindful of practicalfield implementations of the proposed methods, our studies assume limited number ofphasor measurement units (PMUs) installed, in accordance with existing infrastructure inthe BCH system. The features used in our studies are carefully selected to ensure theiravailability in real time. Finally, we consider an adaptive framework for the proposedmethods to be able to comply with topological changes in the system. The trained CARTmodels are tested and show high accuracy rates, and thus, will be able to predict thetransient stability issues of the system under study following different contingencies usingthe synchrophasors obtained from limited number of PMUs in the system. Also, the MARSmodels, which are proposed to be applied for TSA for the first time, show reasonableprediction accuracy rates. Next, we investigate the possibility of an accurate real-timeFDI using synchrophasors received from PMUs installed at the two ends of a transmissionline. In our studies, we are not only content to have the synchrophasor data available, butwe also care about their accuracy, which is critical for a practical application. Therefore,we apply a new metric called goodness-of-fit (GoF), which is calculated over the time spanof measurement and can quantify the credibility of the received synchrophasors. Then,we apply the data to an FDI method to show how accurate and credible the results are.The obtained results show a reasonable relation between the GoF metric, i.e., credibilityof the measured sychrophasors, and the accuracy of the obtained results, validating thesignificance of the proposed method for real-time applications. As it is very rare to have1Throughout this thesis, real-time applications refer to “actions required within one hour or less topreserve the reliability of the bulk electric system”, as defined by North American Electric ReliabilityCorporation (NERC) [1].iiia real power system with all buses and transmission lines equipped with PMUs, we alsopropose a wide-area real-time FDI approach using a linear observer. Through this wide-area approach, we demonstrate the effectiveness of the proposed method by accuratelylocating a fault in a small test system.ivLay SummaryModern systems and equipment have been produced recently for monitoring and controlof power systems. Through these systems, measurements obtained for various electricalquantities in different points of a power system can be transferred to a central locationcalled control centre in real time. Next, some applications should be used to analyzethe received data for an efficient monitoring and control of the power system. Theseapplications should be fast, accurate and reliable. In this thesis, we investigate two ofthese applications, i.e., TSA and FDI. In TSA, we propose to use artificial intelligencetechniques to predict the stability of the power system following a contingency. In FDI,we propose to use new types of metrics to provide the credibility of the obtained faultlocation results. Finally, we propose a new method for FDI to identify and locate thefault in the power system with measurements obtained from limited number of points.The effectiveness of the proposed methods for power system operation and monitoring hasbeen demonstrated through case studies involving small and large-scale test systems.vPrefaceThis thesis has been conducted, written and prepared by the author under supervision ofProf. William G. Dunford and Prof. Yu Christine Chen. There has also been industrialco-supervision for Chapters 1–3, i.e., Giuseppe Stanciulescu from BCH (Now at PowertechLabs) for chapter 1, Dr. Ali Moshref from BBA Power Incorporation for chapter 2, andDr. Farnoosh Rahmatian from NuGrid Power Corporation for chapter 3. Also, AtefehPalizban from Powerex, BCH, gave us valuable practical advice regarding the BCH powersystem base cases and contingencies.In addition to this thesis, there is at least one paper for each chapter already pub-lished, either submitted or currently under preparation. The author of the thesis hasbeen responsible for developing the idea, formulation, simulation and composition of thepapers, aside from the first indicated conference paper in which the thesis author has beena part of a larger team, working on a project at BCH. In this case, the thesis author iscontributing to the finalization and preparation of the conference paper while completinga practicum at BCH.All of the aforementioned publications have been revised and approved by all co-authors. These publications are listed below:Chapter 1• F. Rahmatian, D. Atanackovic, M. Rahmatian, G. Stanciulescu, and M. Nag-pal, “BC Hydro Synchrophasor System for Wide Area Monitoring, Protection andControl − Functional Requirements and System Architecture Considerations,” inProceeding of 2016 CIGRE´ Canada Conference, 17–19 October, 2016.Chapter 2• M. Rahmatian, Y. C. Chen, A. Palizban, A. Moshref, and W. G. Dunford, “Tran-sient Stability Assessment Via Decision Trees and Multivariate Adaptive RegressionSplines,” Electric Power Systems Research, vol. 142, pp. 320-328, January, 2017.• M. Rahmatian, Y. C. Chen, W. G. Dunford, and A. Moshref, “Transient StabilityviAssessment of BC Hydro Power System Using Decision Trees,” in Proceeding of2016 CIGRE´ Canada Conference, 17–19 October, 2016.• M. Rahmatian, W. G. Dunford, A. Palizban, and A. Moshref, “Transient StabilityAssessment of Power Systems Through Wide-area Monitoring System,” in Proceed-ing of 2015 IEEE PES General Meeting (PESGM), 26–30 July, 2015.• M. Rahmatian, W. G. Dunford and A. Moshref, “PMU Based System ProtectionScheme,” in Proceeding of 2014 Electrical Power and Energy Conference (EPEC),12–14 November, 2014.Chapter 3• M. Rahmatian, Y. C. Chen, W. G. Dunford and F. Rahmatian “IncorporatingGoodness-of-fit Metrics to Improve Synchrophasor-based Fault Location,” IEEETransactions on Power Delivery, Available Online, doi: 10.1109/TPWRD.2018.2790410,2018.• M. Rahmatian, Y. C. Chen, W. G. Dunford, and F. Rahmatian, “SynchrophasorData Qualification Measures for Reliable Fault Location,” in Proceeding of 2017Protection, Automation and Control World Conference (pacworld) , 29–31 August,2017.Chapter 4• M. Rahmatian, Y. C. Chen, and W. G. Dunford, “Observer-based Wide-area FaultIdentification and Location” under preparation.viiTable of ContentsAbstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiiLay Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vPreface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viTable of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viiiList of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xiiList of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xivGlossary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xviiAcknowledgments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xixDedication . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xx1 Synchrophasor-based Wide-area Measurement Systems: An Introduc-tion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11.1 Motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11.2 Background . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21.2.1 Structure of a Synchrophasor-based WAMS . . . . . . . . . . . . . . . 21.2.1.1 Phasor Measurement Unit . . . . . . . . . . . . . . . . . . 61.2.1.2 Communication Network . . . . . . . . . . . . . . . . . . . 61.2.1.3 Phasor Data Concentrator . . . . . . . . . . . . . . . . . . 71.2.1.4 Synchronized Measurement Applications or Functions . . . 71.2.2 Applications of a Synchrophasor-based WAMS . . . . . . . . . . . . . 81.2.2.1 Off-line Applications . . . . . . . . . . . . . . . . . . . . . . 81.2.2.2 Real-time Applications . . . . . . . . . . . . . . . . . . . . 91.3 Research Objectives and Anticipated Impacts . . . . . . . . . . . . . . . . . 11viii1.3.1 Objective 1: Real-time TSA using Synchrophasor-based WAMS . . . . 111.3.2 Objective 2: Real-time Impedance-based FDI Using Synchrophasor-based WAMS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 131.3.3 Objective 3: Real-time Wide-area Observer-based FDI Using Synchrophasor-based WAMS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 141.4 Thesis Outline . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 152 Real-time Transient Stability Assessment Using Automatic LearningApproach . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 162.1 Background . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 162.2 Approaches for TSA . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 172.2.1 Time-Domain Approach . . . . . . . . . . . . . . . . . . . . . . . . . 182.2.2 Direct Approach . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 182.2.3 Automatic Learning Approach . . . . . . . . . . . . . . . . . . . . . 202.3 Real-time TSA Using AL Approach . . . . . . . . . . . . . . . . . . . . . . . 212.4 Real-time TSA on the BCH Power System . . . . . . . . . . . . . . . . . . . 222.5 Wide-area Dynamic Security Indices . . . . . . . . . . . . . . . . . . . . . . 232.5.1 Basic Measurement Indices . . . . . . . . . . . . . . . . . . . . . . . 232.5.2 Centre-of-Inertia-Referred Indices . . . . . . . . . . . . . . . . . . . . 232.5.3 Energy Function-Based Indices . . . . . . . . . . . . . . . . . . . . . 252.6 Data-Mining Fundamentals . . . . . . . . . . . . . . . . . . . . . . . . . . . 262.6.1 CART . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 272.6.1.1 Tree Growing . . . . . . . . . . . . . . . . . . . . . . . . . . 272.6.1.2 Tree Pruning . . . . . . . . . . . . . . . . . . . . . . . . . . 282.6.1.3 Optimal Tree Selection . . . . . . . . . . . . . . . . . . . . 302.6.2 MARS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 312.6.2.1 Forward Procedure . . . . . . . . . . . . . . . . . . . . . . 312.6.2.2 Backward Procedure . . . . . . . . . . . . . . . . . . . . . . 332.7 Training . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 342.8 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 362.8.1 Wide-area Dynamic Security Indices Results . . . . . . . . . . . . . 362.8.2 CART Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 372.8.2.1 Comparison with the Existing Methods . . . . . . . . . . . 402.8.2.2 Adaptive Training . . . . . . . . . . . . . . . . . . . . . . . 412.8.3 MARS Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 422.8.4 Sensitivity to Noise . . . . . . . . . . . . . . . . . . . . . . . . . . . . 452.9 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46ix3 Real-time Fault Detection and Identification Using Synchrophasors . 473.1 Background . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 473.2 Approaches for FDI . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 483.2.1 Impedance-based Approach . . . . . . . . . . . . . . . . . . . . . . . 493.2.2 Travelling Wave-based Approach . . . . . . . . . . . . . . . . . . . . 493.2.3 High Frequency Approach . . . . . . . . . . . . . . . . . . . . . . . . 493.2.4 Automatic Learning Approach . . . . . . . . . . . . . . . . . . . . . 503.3 Impedance-based Fault Location through Synchrophasor System . . . . . . 503.4 Preliminaries . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 513.4.1 PMU Modeling . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 513.4.2 Measurement-based Estimation of Line Impedance . . . . . . . . . . 523.4.3 Fault Location and Faulted-phase Identification . . . . . . . . . . . . 553.4.4 Motivation for a GoF Metric . . . . . . . . . . . . . . . . . . . . . . . 573.5 GoF Metric for Synchrophasor Measurements . . . . . . . . . . . . . . . . . 593.5.1 Definition of the GoF Metric . . . . . . . . . . . . . . . . . . . . . . . 593.5.2 Accounting for DC Offset . . . . . . . . . . . . . . . . . . . . . . . . 623.5.3 Effect of DC Offset . . . . . . . . . . . . . . . . . . . . . . . . . . . . 633.5.4 Effect of Fault Initiation Time . . . . . . . . . . . . . . . . . . . . . 633.5.5 Effect of CT Saturation . . . . . . . . . . . . . . . . . . . . . . . . . 663.5.6 Inclusion of GoF Metrics in the PMUs . . . . . . . . . . . . . . . . . . 703.6 Case Studies . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 703.6.1 Simulation Results for Combined Cases . . . . . . . . . . . . . . . . 723.6.2 GoF Criteria for Predicting Fault Location Accuracy . . . . . . . . . 763.7 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 784 Real-time Wide-area Observer-based Fault Detection and Identification 794.1 Background . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 794.2 Preliminaries . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 814.3 Proposed Observer . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 844.4 Case Studies . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 854.4.1 Two-bus System . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 864.4.2 Three-bus System . . . . . . . . . . . . . . . . . . . . . . . . . . . . 884.4.3 Canonical Two-area System . . . . . . . . . . . . . . . . . . . . . . 894.5 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 915 Concluding Remarks . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 925.1 Thesis Summary and Contributions . . . . . . . . . . . . . . . . . . . . . . . 925.1.1 Real-time TSA using Synchrophasor-based WAMS . . . . . . . . . . . 92x5.1.2 Real-time Impedance-based FDI Using Synchrophasor-based WAMS . 935.1.3 Real-time Wide-area Observer-based FDI . . . . . . . . . . . . . . . 945.2 Future Work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 945.2.1 Robustness of the Developed Models for TSA against MeasurementErrors, Bad Data and Missing Values . . . . . . . . . . . . . . . . . 945.2.2 Incorporating GoF Metrics to Improve Synchrophasor-based Appli-cations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 955.2.3 Wide-area Fault Location with Limited Number of PMUs . . . . . . 95Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 97A Parameters of Applied Models . . . . . . . . . . . . . . . . . . . . . . . . 107xiList of TablesTable 1.1 Communication Mediums Available for Data Transfer in Power Systems. 7Table 2.1 TSA indices selected in CART models built using data obtained in 5 cyclesafter fault clearance. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40Table 2.2 Accuracy results of obtained CART models when tested using data fromvarious time instants. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40Table 2.3 Computational times required for various CART models. . . . . . . . . . . 41Table 2.4 TSA indices selected in MARS models built using data obtained in 5 cyclesafter fault clearance. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43Table 2.5 Accuracy results of obtained MARS models when tested using data fromvarious time instants. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44Table 2.6 Computational times required for various MARS models with 15 maxi-mum BFs and 1 maximum interaction. . . . . . . . . . . . . . . . . . . . 45Table 2.7 Accuracy results of the obtained MARS models with 60 maximum BFsand 6 maximum interactions. . . . . . . . . . . . . . . . . . . . . . . . . . 45Table 2.8 Computational times required for various MARS models with 60 maxi-mum BFs and 6 maximum interactions. . . . . . . . . . . . . . . . . . . . 46Table 2.9 Accuracy results of obtained CART models when white Gaussian noise isadded to data set. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46Table 2.10 Accuracy results of obtained MARS models when white Gaussian noiseis added to data set. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46Table 3.1 Fault location and fault-impedance estimation using measurements fromthe third cycle of during-fault period. . . . . . . . . . . . . . . . . . . . . 58Table 3.2 Fault location and fault-impedance estimation using measurements fromthe first cycle of during-fault period. . . . . . . . . . . . . . . . . . . . . 58Table 3.3 Comparison of fault-location accuracy with corresponding GoF metricsfor an abcg fault occurring at 60 km from bus m with Rg = 5 Ω. Faultinitiates at the beginning of measurement window k = 1. . . . . . . . . . 74xiiTable 3.4 Comparison of fault-location accuracy with corresponding GoF metricsfor an abcg fault occurring at 60 km from bus m with Rg = 5 Ω. Faultinitiates at t = 6.94 ms, within measurement window k = 1. . . . . . . . . 74Table 4.1 Fault location results for different locations along the line (1, 2). Theresults obtained in the pre-fault period (k = 4–6) are complex whichdo not make sense for a distance value, while those obtained from theduring-fault period (k = 7–9) are real and correct. . . . . . . . . . . . . . 88Table 4.2 Fault current for different locations along the line (1, 2). The resultsobtained in the pre-fault period (k = 6) are almost zero, while thoseobtained in the during-fault period (k = 7–9) are large values. . . . . . . 88Table 4.3 Fault location results for different lines and locations in the three-bussystem. The results obtained in the pre-fault period (k = 4–6) are com-plex, which do not make sense for a distance value, while those obtainedfrom the during-fault period (k = 7–9) are correct. . . . . . . . . . . . . 89Table 4.4 Fault location results for different lines and locations in the three-bussystem. The results obtained in the pre-fault period (k = 6) are almostzero, while those obtained in the during-fault period (k = 7–9) are largevalues. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 89Table 4.5 Fault location results for different lines and locations in the canonicaltwo-area system. The results obtained in the pre-fault period (k = 4–6)are complex, while those obtained from the during-fault period (k = 7–9)are correct. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 90Table 4.6 Fault location results for different lines and locations in the canonicaltwo-area system. The results obtained in the pre-fault period (k = 6)are almost zero, while those obtained in the during-fault period (k = 7–9)are large values. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 90Table A.1 Parameters of the exciter model ESAC4A. . . . . . . . . . . . . . . . . . 107Table A.2 Parameters of the governor model GOV6. . . . . . . . . . . . . . . . . . . 107Table A.3 Parameters of the power system stabilizer model PSS1. . . . . . . . . . . 108xiiiList of FiguresFigure 1.1 A simplified conceptual architecture of a typical synchrophasor-basedWAMS showing various functional layers/tires [2]. . . . . . . . . . . . . . 4Figure 1.2 Connectivity diagram of hardware components of a typical synchrophasor-based WAMS [2]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5Figure 2.1 Different types of power system stability phenomena. Adapted from [3]. 17Figure 2.2 General automatic learning framework for TSA. . . . . . . . . . . . . . 20Figure 2.3 Plot of node impurity functions for two-class classification, as a functionof the probability p(z). Entropy-based impurity function is scaled topass through (0.5,0.5). . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29Figure 2.4 Diagram of BCH and its interconnections to BPA and AES [4]. . . . . . 35Figure 2.5 A transiently unstable case following a fault occurrence. . . . . . . . . . 36Figure 2.6 Characteristic of two COI-referred indices following different faults inthree different power flow patterns. . . . . . . . . . . . . . . . . . . . . . 37Figure 2.7 CART model trained and built using synthetic measurements obtainedin the fourth cycle after fault clearance. . . . . . . . . . . . . . . . . . . 39Figure 2.8 CART model trained based on pre-fault synthesized features obtainedfrom SCADA system. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41Figure 2.9 Adaptive framework for online TSA. . . . . . . . . . . . . . . . . . . . . 42Figure 2.10 Sensitivity analysis of MARS model trained using synthetic measure-ments obtained in the second cycle after fault clearance to maximumnumber of interactions and maximum number of BFs. . . . . . . . . . . 44Figure 3.1 General diagram of a PMU for one phase. . . . . . . . . . . . . . . . . . 52Figure 3.2 PMU model simulated in MATLAB Simulink. . . . . . . . . . . . . . . . 53xivFigure 3.3 Circuit model of three-phase transmission line (m,n) (see, e.g., [5])(a) before fault, and (b) during fault. The fault location (character-ized as the distance d from bus m) and the fault-to-ground impedances(Za(m,n),f , Zb(m,n),f , and Zc(m,n),f ) are unknown quantities. Once solved,they provide the fault location and identify the affected phases. . . . . . 54Figure 3.4 Part of two-area power system with VT, CT and PMU models simulatedin MATLAB Simulink. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55Figure 3.5 (Example 3 and 4) Phase-a actual current waveform and that recon-structed from PMU phasor measurements (a) without and (b) withDC-offset modification. Corresponding GoF values are superimposed.Fault initiates at the beginning of measurement window k = 1. . . . . . 60Figure 3.6 (Example 3 and 4) Phase-a actual voltage waveform and that recon-structed from PMU phasor measurements (a) without and (b) withDC-offset modification. Corresponding GoF values are superimposed.Fault initiates at the beginning of measurement window k = 1. . . . . . 64Figure 3.7 (Example 5) Phase-a actual current waveform and that reconstructedfrom PMU phasor measurements (a) without and (b) with DC-offsetmodification. Corresponding GoF values are superimposed. Fault initi-ates at t = 6.94 ms, within measurement window k = 1. . . . . . . . . . 65Figure 3.8 (Example 5) Phase-a actual voltage waveform and that reconstructedfrom PMU phasor measurements (a) without and (b) with DC-offsetmodification. Corresponding GoF values are superimposed. Fault initi-ates at t = 6.94 ms, within measurement window k = 1. . . . . . . . . . 66Figure 3.9 (Example 6) CT flux saturation limit is set to 2µV·s. (a) Phase-a CTmagnetizing current and core flux. (b) Phase-a actual current wave-form and that reconstructed from PMU phasor measurements without(c) with DC-offset modification. Fault initiates at the beginning ofmeasurement window k = 1. . . . . . . . . . . . . . . . . . . . . . . . . . 67Figure 3.10 (Example 6) CT flux saturation limit is set to 2µV·s. (a) Phase-a CTmagnetizing current and core flux. (b) Phase-a actual current wave-form and that reconstructed from PMU phasor measurements without(c) with DC-offset modification. Fault initiates at t = 6.94 ms, withinmeasurement window k = 1. . . . . . . . . . . . . . . . . . . . . . . . . . 68Figure 3.11 Schematic of proposed PMU. . . . . . . . . . . . . . . . . . . . . . . . . 70Figure 3.12 GoF versus error for various fault types in original form and after DCoffset remove . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71xvFigure 3.13 Fault location error versus GoF values computed for waveforms recon-structed from PMU phasor measurements obtained in measurementwindows k = 1, 2, 3 (a) without and (b) with DC-offset modification.The × markers denote cases in which fault initiates 12.5 ms (equivalentto 3/4 of the measurement window length) into measurement windowk = 1. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73Figure 3.14 Number of cases versus GoF ranges (a) in original form and (b) afterDC-offset remove for simulated cases. . . . . . . . . . . . . . . . . . . . 74Figure 3.15 Average error and SDE versus GoF ranges (a) in original form and (b)after DC-offset remove for simulated cases . . . . . . . . . . . . . . . . . 75Figure 3.16 Proposed decision tree to assess credibility of fault location. . . . . . . . 77Figure 3.17 Proportion of cases corresponding to particular ranges of fault locationerrors out of a total of (a) 445 cases with GoF > 25 dB, (b) 89 cases withGoF < 25 dB and GoF > 25 dB, and (c) 366 cases with GoF < 25 dB andGoF < 25 dB. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77Figure 4.1 General diagram of fault. . . . . . . . . . . . . . . . . . . . . . . . . . . 82Figure 4.2 Two-bus system used for testing the proposed method. . . . . . . . . . . 86Figure 4.3 Amplitude of residual current (e[k]) for two different values of µ , i.e.,µ = 0.9, and µ = 0.01, for an abcg fault imposed at 0.25` of line (1,2)with the length of ` in the two-bus system. The fault initiates immedi-ately after time instant kf = 6 and is denoted by a bold dashed line. . . 87Figure 4.4 Three-bus system used for testing the proposed method. . . . . . . . . . 89Figure 4.5 Two-area system used for testing the proposed method. . . . . . . . . . 90xviGlossaryAC alternative current. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .6AES Alberta Electricity System . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22AL automatic learning . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiiBCH BC Hydro . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiiBF basis function . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31BFs basis functions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31BPA Bonneville Power Administration . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12CART Classification and Regression Tree . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiiCOI centre-of-inertia . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12CT current transformer. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .3DC direct current . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62DT decision tree. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .12ELM ensemble learning machine . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21EMS energy management systems . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1FACTS flexible alternating current transmission systems . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9FFT Fast Fourier Transform . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51FDI fault detection and identification . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiiGoF goodness-of-fit . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiiGPS Global Positioning System . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3HVDC high voltage direct current . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .9IEDs intelligent electronic devices . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3IEEE Institute of Electrical and Electronics Engineers . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6MARS Multivariate Adaptive Regression Splines . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiiMAE mean absolute error . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44xviiMAEs mean absolute errors . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43NERC North American Electric Reliability Corporation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiiNTC nominal transfer capability . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9PDC phasor data concentrator . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3PMU phasor measurement unit . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xixPMUs phasor measurement units . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiiRF random forest . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21ROCOF rate of change of frequency . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3RTs regression trees . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12RTUs remote terminal units . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1SCADA supervisory control and data acquisition . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1SNR signal-to-noise ratio . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45SPM Special Predictive Modeler . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36SDE standard deviation of error . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76T-D time-domain . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11TSA transient stability assessment . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiiTSI transient stability index . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12TVE total vector error . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50VT voltage transformer . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3WADSIs wide-area dynamic security indices . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27WAMS wide-area measurement systems. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiiWECC Western Electricity Coordinated Council . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiixviiiAcknowledgmentsIt has been my privilege to be advised, supported, and encouraged by a lot of peoplethroughout the course of my Ph.D. studies.First and foremost, I would like to thank my dear supervisors, Prof. William G.Dunford and Prof. Yu Christine Chen, for their kind support, patience, and advice. Prof.Dunford was always encouraging and supportive. He provided me with the opportunityfor industrial experience, which eased the process of finding a job at the end of my Ph.D.program. Prof. Chen also dedicated much of her time to evaluate my work and providenovel ideas, deep insight, and invaluable feedback into my research. Her concentrationand dedication to her job will be an unforgettable lesson for me throughout my futureprofessional life.I would also like to acknowledge Dr. Ali Moshref from BBA Incorporation for the timehe spent training me on the DSATools package and for providing me with real models ofthe BCH power system. Dr. Moshref was like an older friend to me, and helped to showme the bigger picture of my life.I’m also grateful to my cousin, Dr. Farnoosh Rahmatian from NuGrid Power Corpo-ration, for his brilliant ideas and technical support on phasor measurement unit (PMU)models and WAMS specs. I could always rely on his advice and help during the past fouryears.My sincere thanks also goes to Giuseppe Stanciulescu from Powertech Labs and HelenWhittaker from BCH for their kind support in presenting Chapter 2 of this thesis to BCH.I learned a lot about BCHs road map for WAMS at the Office of Applied Technology, BCH.I would also like to thank my friends in the Electric Power and Energy Systems Groupin the Department of Electrical and Computer Engineering at The University of BritishColumbia, for their help and support: Shreya Iyer, Javier Tarazona, Hamed Ahmadi, andPouya Zadkhast. Especially, Navid Amiri was always helpful to me technically.Finally, I would like to profoundly give thanks to my wonderful and dear parents andmy beloved wife for the encouragement, support, patience and fun moments they providedduring the course of my Ph.D. studies.xixDedicationTo my dear parents, Nahid and Alireza, and beloved wife, Shirin.xxChapter 1Synchrophasor-based Wide-areaMeasurement Systems: AnIntroduction1.1 MotivationThe operation of large-scale power systems has become more challenging in recent years.Increasing interconnections with other power systems, restructuring and electricity mar-ket requirements, load growth, and increasing renewable energy resources have all madethe operation of systems more complicated. At the same time, utilities try to use pre-established capacities and infrastructures more efficiently. These facts, along with recentcascading blackouts in some power systems around the globe, have all resulted in utilitiesinvesting in modern energy management systems (EMS).Traditionally, utilities use supervisory control and data acquisition (SCADA) systemsto monitor and control power systems. In these systems, signals are obtained from remoteterminal units (RTUs), relays, and other meters, and then are transferred to a control centrethrough a communication infrastructure. This process takes 2 − 10 seconds in differentutilities [6]. Moreover, the measurements reaching the control centre from different partsof the power system would not be synchronized. The relays and RTUs at the substationsuse local clocks to tag time to the measurements, leading to inconsistent time tags at thecontrol centre. Therefore, the operators at the control centre have to use state estimatorsfor the power system and make sure all the states are obtained correctly. We should alsonote here that the most important dynamic states of a system are voltage phase angles,which are necessary for monitoring fast stability issues. As a result, the operators arenot able to monitor and control the dynamics of the power system in real-time, or even1near real-time, using SCADA-based monitoring systems. This deficiency can lead to somemajor problems in a system, specifically in peak load cases. For instance, the blackoutexperienced in the Northeastern US on August 14, 2003 could have been prevented if theoperators at the control centre had an improved observation of the system. As a result,one of the suggestions regarding this event was to use time-synchronized measurements(synchrophasors) [7]. Consequently, utilities have started shifting from SCADA systems tomodern synchrophasor-based monitoring systems called WAMS in recent years.Introducing synchrophasor-based WAMS to the electricity industry has led to new ca-pabilities in the operation of power systems. However, utilities have their own practicallimitations for utilizing this technology. In many cases, the most important concern of util-ities is their financial limitations, and establishing a modern synchrophasor-based WAMSand utilizing its capacities at the control centre can be very costly. Expensive structureand expertise from different fields would be needed and investment for training humanresources would also be indispensable. Consequently, many utilities develop road mapsand deployment plans for this new technology and define various phases to utilize it in along-term schedule. However, more research and development is still needed to make thisnew technology completely operational and express its usefulness to the utilities and EMS,even when not established completely [2].In this thesis, we first review the structure and potential applications of a modernsynchophasor-based WAMS. Then, we focus on two real-time applications of a typicalsynchrophasor-based WAMS to express the potential . First, we perform a TSA on thereal power system of BCH and show the possibility of predicting unstable cases through amodern synchrophasor-based WAMS with a limited number of PMU. The goal of this studyis to show the potential of WAMS for real-time TSA, even without a complete observability ofthe system, and compare it with the SCADA system to reveal its advantages over traditionalmethods. Next, we focus on the practical capability of the modern synchrophasor systemto create a credible fault location from a central location. This credible fault locationbased on synchrophasors at a central location is again another capability which has notbeen previously studied.1.2 Background1.2.1 Structure of a Synchrophasor-based WAMSNowadays, utilities are installing PMUs in large numbers in power systems. However, toexploit the equipment completely, this should be seen as a part of WAMS [8]. A typicalhierarchical structure of a WAMS is shown in Fig. 1.1 and Fig. 1.2, and is composed of thefollowing parts [2]:2 PMU: Field measurement devices, known as PMU, extract various prescribed quantitiesfrom the time-domain signals collected by current transformer (CT) or voltage trans-former (VT) installed at a substation; extracted quantities include phasor amplitude,phase angle, frequency, and rate of change of frequency (ROCOF). Time synchronization (or syntonization) devices (clocks): Time synchronization devicesprovide accurate time usually based on Global Positioning System (GPS) for use byPMUs in time stamping measurements at the source. The extracted quantities are allsynchronized and tagged with an accurate coordinated universal time stamp [9, 10]. Sub-phasor data concentrator (PDC): Optionally, local PDCs or sub-PDCs are used tocombine and align data from various PMUs. Communications network : Digital communication infrastructure typically uses IP-basedswitches, routers, firewalls, and fibre-optic cables to transfer data [11]. Super-PDC: PDCs are utilized for aggregating, time-aligning, and re-sampling largeamounts of synchrophasor data in central locations such as transmission system controlcenters [12]. Synchrophasor-based applications or tools: Various application software run on com-puters or intelligent electronic devices (IEDs), either locally near/at PMUs or at centrallocations (such as the control centre), fed by PDC and Super-PDCs. Examples of crucialreal-time applications that can benefit from WAMS are TSA and fault location along atransmission line. [13].Fig. 1.1 shows a simplified conceptual architecture of a WAMS infrastructure, showingvarious functional layers/tiers. PMUs are generally located in substations (or other criticallocations) where synchrophasor data is acquired. Local PDCs are proposed for aggregatesites or larger substations (sub-PDC). Use of sub-PDCs can reduce the data loss (for off-linefunctions) associated with sporadic communication failures. Sub-PDCs can also serve localsubstation applications using synchrophasors. Wide-area applications, main PDCs (super-PDCs) and archives are located in central locations. The central applications and serversare divided into EMS and non-EMS application environments to simplify security compli-ance requirements. Redundancy in each environment is not shown. Most applicationsinteracting with EMS and real-time grid control (and SCADA) are placed in the secure EMSenvironment and are fed from the super-PDC located in the same environments. Engineer-ing and monitoring applications that need to be accessible to non-control-centre staff arelocated in a separate secure non-EMS environment that is fed from different super-PDCs [2].3 Secure Non-EMS EnvironmentPMU* BDPMU* ACSub-PDCSub-PDCData Archives (External to EMS)PhasorAnalysisApplicationsStateEstimatorPrimary SecondaryContingencyAnalysisDetermine recommended system operator actions and communicate to operators/dispatchers within 2 secondsRecommendedActionsMonitoring & Control ApplicationsAlarmsPEAK-RC /WECC /External DataPhasor Data Collection TierFunctions: EMS & State Estimation Remedial Action Schemes (RAS) Dynamic Stability Indication Contingency Analysis Situational Awareness (Visualization and Alarms) Post-Disturbance Event Analysis Fault LocationFunctions: Collect accurate, high-resolution multi-purpose synchrophasor data*The PMU function may be provided by dedicated devices and/or multi-function IEDs (e.g., relays, recorders, meters, …)VisualizationSuper PDC(tied to EMS)SCADAUp to 120/s Secure EMSEnvironmentOther Synchrophasor functionsTo a common Visualization toolTo a common Visualization toolPhasor Data StorageData Aggregation TierApplication Interface TierApplication TierFunctions: Aggregate and concentrate data Local storage (e.g., 2 months)Functions: Aggregate and concentrate data Provide archival function Provide Interface for data exchange with external entitiesSubstationData ManagementAnd ApplicationsEMS-RAS(Automatic Arming) and VSA, TSA, ...Super PDC(Non-EMS)Real Time (within 100msec)EMSData ArchivesPhasor-based Dynamic Stability Analysis Data StorageOnce per 4 secondsPlanning and Engineering ApplicationsFigure 1.1: A simplified conceptual architecture of a typical synchrophasor-based WAMS showing various functional layer-s/tires [2].4Control Center CDControl Center ABOther Applications: Voltage stability indication Transient stability indication Oscillation DetectionOther Applications: Voltage stability indication Transient stability indication Oscillation DetectionExisting Gateway(Remote Access) RAIDStorageRAIDStorage60 or 120 phasors/sRAIDStorageRAIDStorage Example of a Substation with local PDC and a single communication pathRAIDStorage RAIDStorage60 or 120 phasors/s Example of Small Substation without PDC Example of a Substation with an aggregating PDC and redundant communication pathsAnother Aggregate SitePMU BDPMU BDPMU ACPMU ACPMU AC PMU AC PMU BDSubstation PDC (BD)Substation PDC (AC)Substation PDC (AC)Substation PDC (BD)Substation PDC (AC)Substation PDC (BD)LAN ACLAN BDLAN BDLAN AC60 or 120 phasors/s30-60 phasors /s60 or 120 phasors/sEMS D (CC2)Super-PDC A (CC1)Super-PDC C (CC2)EMS C (CC2)RAIDStorageRAIDStorageSuper-PDC D (CC2)RAIDStorageRAIDStorageOther Utilities and PEAK RCOther Utilities and PEAK RC60 or 120 phasors/sRouters and firewallsPMU BDSuper-PDC B (CC1)Edge routers and firewall60 or 120 phasors /sEMS A (CC1)EMS B (CC1)Not All Substations have Redundant Communication ChannelsFigure 1.2: Connectivity diagram of hardware components of a typicalsynchrophasor-based WAMS [2].The PMU fits a phasor to a time-domain signal collected by a current or voltage trans-former, and further reports the phasor amplitude, phase angle, frequency, and ROCOF [9].These quantities are all synchronized and tagged with an accurate coordinated universaltime stamp by using, for example, the reference time acquired through global positioningsystem receivers [9]. Thereafter, a wide-area communication network enables synchropha-sor measurements to be transmitted from various PMUs to PDCs and then a central loca-tion [11]. These measurements may be used by real-time application deployed at the levelof a single PMU, a local PDC, or a central decision-maker [13].Fig. 1.2 shows a connectivity diagram for hardware components of a generic syn-chrophasor system, with focus on the parts of synchrophasor systems and applicationshosted in the EMS security domain (non-EMS environment not shown). This shows re-dundancy in data acquisition, communications, and application servers. In the figure,emphasis is on having quad-redundant application servers, labeled A, B, C, and D. The Aand B systems are hosted at one control centre, while C and D are at the other grid controlcentre. Field devices can have dual redundancy (including dual communication paths):the AC PMUs will feed the A and C application servers at the two control centre locations5and the BD PMUs will feed servers B and D. Fig. 1.2 also shows that not all substations(especially smaller ones) necessarily have sub-PDCs or redundant communication links [2].Below, we discuss WAMS components with more details.1.2.1.1 Phasor Measurement UnitPower systems are basically alternative current (AC) circuits, where every signal measuredwould be sinusoidal. A sinusoidal signal can be written in the form ofx(t) = Acos(2pift+ φ) (1.1)in which A is the signal amplitude, f is the signal frequency, t is time and φ is the signalphase angle. For the phasor form of the sinusoidal signal with frequency of f at time t,one can writeX =A√2∠φ. (1.2)The function of PMU is to obtain the measured signal form the CT or VT installed at thesubstation and measure the pre-desired quantities, e.g., amplitude, phase angle, frequencyand ROCOF. Next, the PMU synchronizes the obtained phasor to coordinated universaltime (UTC) as the reference time received through GPS receivers [8]. The Institute ofElectrical and Electronics Engineers (IEEE) has issued two standards about synchrophasormeasurements and PMUs, C.37.118.1 and C.37.118.1a, in which minimum requirements forsynchrophasors and various types of PMUs are discussed [9, 10]. We will elaborate abouton the details of PMU models in Chapter 3.1.2.1.2 Communication NetworkA communication network is needed to transfer synchrophasors from the PMUs to the PDC,the application location and the control centre. PMUs are located in various substationsaround the system, far from the control centre. To have a modern WAMS, the communi-cation delay should be minimum. The communication delay depends on the propagationdelay, network bandwidth and the protocol used for data transfer. The propagation delayis highly dependable on the communication medium. Table 1.1 shows the available mediacurrently in operation in various utilities around the globe. Many utilities have showninterest in fiber optic cables in recent years to establish WAMS. There is also an IEEEstandard called C37.118.2 for synchrophasor data transfer [11, 14].6Table 1.1: Communication Mediums Available for Data Transfer in Power Systems.Communication Medium Advantage DisadvantageTelephone Line Economical, available, easy to set up. Slow propagation speed.Satellite Available. Costly, narrow bandwidth.Power Line Communication Economical. Signal decay and distortion.Microwave Easy access to remote substations. Signal decay and useless propagation.Fiber Optic Cable Wide bandwidth, Fast propagation speed, Costly.immune to electromagnetic interference.1.2.1.3 Phasor Data ConcentratorAccording to the IEEE standard C37.244, PDCs are nodes of communication networks.PDCs can be installed at a substation or at the control center of a power system, and aPDC may be built stand-alone or combined with another device [12].Various functions have been defined for a PDC in the power system. Some of the mostimportant of these functions are listed below [12]: Data aggregation: Aggregating data with preserving data and its time quality, Data forwarding : Sending data from one input to multiple outputs, Data communications: Connecting with other devices based on a serial or Ethernetnetwork, Data validation: Checking status flags of received data, time quality of the connectedPMUs, and integrity of received data, Data transfer protocol support : Supporting different protocols used by PMUs to transfersychrophasors’ data, Data transfer protocol conversion: Converting different protocols to each other, Data latency calculation: Calculation of the time between sending and receiving a datapacket in the network, Configuration management : Preparing data for local functions as well as those at thecontrol center.1.2.1.4 Synchronized Measurement Applications or FunctionsModern synchrophasor-based WAMS have made many functions feasible or more efficientin power systems, whether at the control centre or at substations. Generally, applica-tions developed based on modern WAMS can be divided into two categories, i.e., off-lineapplications and real-time applications.7Off-line applications can generally tolerate fairly long data communication latency, butrequire a high grade of data availability. On the other hand, real-time (and near-real-time)applications can range from requiring very low latency for data communication betweenPMUs and the synchrophasor application servers (e.g., TSA or system protection and real-time control) to being able to tolerate moderately long latencies (e.g., fault location andtotal transfer capability determination). Most of these real-time applications can startbeing useful even at lower application availability (and reliability) levels; they can helpimprove grid monitoring, protection, and control. Nevertheless, as they mature and be-come a critical part of EMS, they need to provide availability levels similar to the presentEMS-based state estimator applications having better than 99.95% availability with eachdown-time not to exceeding 15 to 30 minutes [2, 15].1.2.2 Applications of a Synchrophasor-based WAMS1.2.2.1 Off-line Applications Post-event analysis and data storage: Post-event analysis is defined as investigationof an event through reconstructing its sequence, identifying responses of associateddevices, and evaluating system performance. During complex events, multiple factorsmay contribute to a single event in a short interval. Modern synchrophasor-based WAMSand data storage can save a lot of time and effort in analysis. This process could takeseveral months using data obtained from the traditional SCADA system [16]. System baselining and parameterization for EMS: This application is defined as theanalysis of buck synchrophasor data and identification of nominal and off-nominal op-erating conditions of the system under monitoring. This helps engineers to obtain aninsight into the system’s behavior over a long period of time, based on real operatingconditions which would be more reliable than model-based studies and event-based ex-perience. The obtained observations can be used in the planning stage and developingoperational guidelines [16]. System protection and control planning : The data obtained through synchrophasor-based WAMS can be used to develop more advanced system-wide special protectionschemes, symtem integrity protection schems, and remedial action schemes, consideredas a part of the planning stage for system protection. Some examples include load shed-ding, generation shedding, and system controlled islanding and resynchronization. Thiswould be a large advance compared to conventional and local protection schemes [17]. Dynamic model validation: Frequent disturbances in power systems are unavoidable.Through modern synchrophasor-based WAMS, coincidental and controlled disturbances8can be captured to be utilized in validating and calibrating a dynamic model of systemequipment such as generators, exciters, governors, power system stabilizers, flexible al-ternating current transmission systems (FACTS), and high voltage direct current (HVDC)systems [18].1.2.2.2 Real-time Applications Situational awareness and visualization: Modern synchrophasor-based WAMS can pro-vide a better dynamic insight for system operators over a wide geographical area of thesystem, so that they can operate and manage the system more efficiently, anticipate po-tential problems and prevent or control the potential problems before disruption arises.The key parameter for operators here would be to have a display showing phase anglesand voltage magnitudes, in addition to key transmission lines’ flow and frequency ofcritical buses in near real-time which has become feasible through modern WAMS [16]. Enhanced state-estimation: As discussed in Section 1.1, conventional state-estimatorsare based on SCADA measurements which do not include synchronized measurements andreach to the control centre in a time range of seconds. Emerging modern synchrophasor-based WAMS has helped industry to improve hybrid state-estimators through the inclu-sion of phase angles in addition to SCADA measurements, and increase the accuracy andconvergence characteristics of state-estimators. In the case of expanding a WAMS to pro-vide a full observability of the system, state-estimators can calculate all system statesquickly and accurately through linear equations. These new types of state-estimatorsare called phasor-only state-estimators [13]. Transfer capability assessment : The current approach for real-time congestion man-agement is to compare actual flow on a transmission line with the nominal transfercapability (NTC), obtained through off-line studies. The NTC is the minimum valueamong thermal limit of the line, voltage stability and transient stability limits. Ther-mal rating for a transmission line is considered a fixed value during summer or winterseasons. Some variables such as ambient temperature, wind speed and solar heating areused as input to calculate a conservative value for transmission line conductor capacity.In a SCADA-based EMS, the congestion management tool uses the state-estimator resultsto control the level of power flow in a line to maintain the system reliability and lowersystem operation costs. Modern synchrophasor-based WAMS can improve performanceof the congestion management tool by providing real-time data of power flow [13]. Oscillation monitoring and control : Power oscillations may occur in long transmissionlines connecting two different areas or transferring power within an area from remotegeneration units to load centres. Synchrophasor-based WAMS can provide facilities for9real-time detection of oscillations, identification of mode shapes throughout intelligentprocessing and alarm the operator. Some oscillations may be controllable through activeequipment such as FACTS, brake insertions or phase shifting transformers [16]. System protection and real-time control : Currently, event-driven discrete-time controlschemes, e.g. system protection schems or remedial action schemes, are installed inpower systems. A modern synchrophasor-based WAMS makes continuous feedback con-trol systems using synchronized data technically feasible in power systems such as real-time control of generators, HVDC systems, FACTS and other devices [17]. Voltage security assessment and control : Monitoring voltage-related conditions, i.e.,voltage profile, voltage sensitivities and reactive power reserves using synchrophasorshelps system operators to observe voltage security issues in real time. Furthermore,synchrophasor-assisted software tools can be applied to take automated or manual ac-tions to control the system and prevent probable catastrophes. These control actionsinclude adjusting set point or switching in/out reactive power compensation devicessuch as exciters, static VAR compensators, synchronous condensers, shunts and capac-itor banks [13]. Frequency stability monitoring and control : Modern synchrophasor-based WAMSs pro-vide the capability of real-time system-wide frequency monitoring and trending andgive the operator a situational awareness on the seriousness of the event and systemfrequency excursions following an event. Frequency deviations exceeding certain limitscan also initiate control practices or set points of governor response or the capability ofthe automatic generation control in a short period of time [16]. TSA and control : Synchronized phase angles provided through WAMS can be used forassessment of transient stability, also called angular stability. Synchrophasor-assistedapplications can be set up at the control centre to monitor phase angle difference at keysubstations of the system, such as large generating units, two ends of important powercorridors and interconnections with neighbouring utilities, to alarm the operators inreal-time in case of angle differences exceeding certain limits. In the next step, certaincorrective control actions can be provided in real-time, either automatically or manually,to maintain the system’s angular stability [13, 15]. FDI: Accurate real-time FDI in the transmission system is very beneficial for systemoperators and customers in terms of reducing power system outage, restoration timeand costs and improving the system’s reliability. Real-time synchrophasors providedthrough WAMS can be used by an application at a central location for fault locationbased on impedance approach with an acceptable accuracy and confidence-level, with10which the crew can target a limited location range and save a lot of time. Since faultlocation is a fast phenomenon and is usually cleared after a few cycles in a transmissionlevel, P-class synchrophasors should be used by the respective application [19]. Intelligent islanding and Resynchronization: Controlled system islanding or separationis considered as a last resort to prevent a power outage in an uncontrolled fashion.To avoid a wide-spread blackout following severe disturbances, the control centre mayseparate the transmission system in a controlled manner. But, it should be ensured that,in each island, the remaining load and generation are sufficiently balanced that powercan be transferred through the remaining transmission system. Intelligent islandingcan also make the restoration and resynchronization process easier in comparison to acomplete blackout [17]. Renewables integration: Using synchronized data provided by modern WAMS and newassociated data storage capabilities, the generational and operational patterns of vari-able renewable generation units, e.g. wind and solar, and their interactions with conven-tional generation units can be studied precisely for planning and integration purposes.Moreover, through real-time facilities provided by WAMS, renewable resources can bemonitored and controlled more efficiently [18]. Real-time data sharing with neighbouring utilities: Through modern WAMS, criticalmeasurements can be shared between neighbouring utilities such as data related tointerconnections and big power plants[16].1.3 Research Objectives and Anticipated ImpactsReal-time applications of synchrophasor-based WAMS have received attention from re-searchers and engineers. Applications at a central location should be able to deliver theresults in near real-time and with a high accuracy rate. Having this in mind, we havepicked two examples of real-time applications, i.e., TSA, and FDI, and have proposed newalgorithms that are capable of responding both fast and accurately based on the receivedsynchrophasors. In both applications, we have attempted to solve a problem which is alsoof important concern for industry.1.3.1 Objective 1: Real-time TSA using Synchrophasor-based WAMSGiven the obvious utility of TSA, various approaches to accomplish this have been ex-plored in the literature that can be divided into three main categories: time-domain (T-D)approaches, direct approaches, and AL approaches. The main drawback of the T-D ap-proaches is that they are computationally burdensome [20]. On the other hand, in direct11approaches, energy functions are difficult to construct for practical large-scale systems [21].In light of these drawbacks, in this thesis, we pursue an AL approach. Recent work hasrevealed AL approaches to be promising for real-time TSA [22–27]. Unlike T-D and directapproaches, AL approaches have been successfully tested on real large-scale power systems[26, 28]. For example, in [25], a comprehensive database of the Hydro-Que´bec power sys-tem along with relevant contingency cases are constructed and a decision tree (DT) is usedto classify stable vs. unstable scenarios. Another classifier is the random forest, which hasbeen shown to operate efficiently in the presence of small changes in network topology [29].In [22], regression trees (RTs) are used to indicate the number of overloaded lines and buseswith voltage magnitude violations following a contingency. In [23], application of multi-ple DTs leads to improved reliability for determining stable and unstable scenarios. Thesupport vector machine method is used as the classifier in [27], where energy-based powersystem features are applied. Further, high-accuracy prediction rates are demonstrated,albeit only for small-scale systems (New England 39-bus power system model), in [27].In [24] and [30], in order to tailor for online applications, an ensemble learning machine isused to decrease the time of training and decision-making.In Chapter 2, we focus on the practical implementation of classification and regressionmodels to assess transient stability in real time using PMU measurements obtained fromthe system. The classification model is developed using CART, which is a well-known DTalgorithm [31]. The most important advatange of DT algorithms is their transparency andtheir simiplicity for interpretation. We then use the developed model to classify test casesas either stable or unstable with a high degree of accuracy. In addition to CART, whichreturns discrete-valued stable versus unstable classifications, we propose to use MARS [32]in conjunction with a continuous-valued transient stability index (TSI) to assess the sever-ity of the contingency and the system’s proximity to instability. The MARS method isalso transparent and easy to interpret. Application of this method can also be usefulto determine probable corrective actions. Development of both CART and MARS modelsproceeds as follows. First, we identify system characteristics that are pertinent to TSA,including voltage and current phasors, deviations from centre-of-inertia (COI) angle andspeed, and potential- and kinetic-energy related quantities. Next, these characteristicquantities, called “indices” or “features”, are evaluated using synthetic measurements ob-tained from time-domain simulations of the full WECC system. These evaluated quantitiesare used to build classification and regression models, which are then tested and veri-fied on the WECC system, with emphasis on transient stability issues in the BCH powersystem that may emerge due to high power exchange rates with the Bonneville PowerAdministration (BPA) system. We also show the necessity and practicality of updatingthe obtained models after topological changes or in the presence of high penetration ofrenewable resources.12To validate the proposed CART and MARS methods in a practical setting, we conducttest cases on the full WECC system model, with emphasis on the BCH power system. Itis worth noting that this is the same detailed model that BCH uses for its planning andoperations studies. In light of practical field implementations of the proposed methods,our studies assume limited number of PMUs installed, in accordance with existing infras-tructure in the BCH system. Similarly, we assume only large generating plants of theneighbouring regions are equipped with PMUs. As a consequence of the limitation in thenumber of installed PMUs, some features, such as the location of the fault highlightedin [33], may not be readily available. Moreover, certain features proposed in literaturemay not be available in real time, such as the duration of the fault and the time to loss ofsynchronism considered in [34]. Since this chapter focuses on the practical implementationof CART and MARS methods, features used in our studies are carefully selected to ensurethey are available in real time, while being cognizant of existing infrastructure limitations.In addition to considering practical limitations, we also take advantage of recent advancesin PMU technology and assume that synchrophasor measurements are available at thecontrol centre through WAMS at every electrical cycle [16].1.3.2 Objective 2: Real-time Impedance-based FDI UsingSynchrophasor-based WAMSIn this section, we focus on practical issues that a typical synchrophasor-based FDI methodencounters. Given the potential of PMU measurements to realize and improve fault loca-tion, numerous methods have been proposed to solve this problem in known literature [35–40], all categorized as impedance-based approaches. There are two main categories ofimpedance-based fault location methods based on PMU measurements, i.e., single-end[35, 36] and double- or multi-end methods [37, 38]. Single-end methods utilize measure-ments obtained from one bus and do not need to obtian any data from remote buses.Therefore, they have high reliability, but the obtained results do not have enough accuracy.On the other hand, communication technologies available for power systems monitoringhave improved over time. Therefore, it has become possible to use measurements from twoends of transmission lines or multi ends in the case of multi-terminal lines, to obtain moreaccurate results. Although the two- or multi-end algorithms result in higher accuracy,unsynchronized measurements obtained from remote buses may cause some problems thatare taken into account in the corresponding literature [39, 40]. Utilizing PMU measure-ments, i.e., synchrophasors, would solve such problems [41], which is a clear advantage ofmodern synchrophsor-based systems.The total fault occurrence procedure is divided into pre-fault, during-fault and post-fault periods. PMU measurements are expected to be precise in the pre-fault period since ittakes place in a steady-state operating condition. But as the fault initiates, the PMU may13not be able to measure the signal perfectly as it enters into transient-state conditions.The transient-state lasts until the line disconnection or until the first few cycles afterthe fault clearance time. It should be noted that the applied inputs to a fault locationmethod use current and/or voltage synchrophasors of the during-fault period. It shouldbe investigated if the fault application can rely on the synchrophasors.Although numerous approaches have been proposed for using phasor measurementsfor fault location [38–40, 42–44], they do not study the effect of measurement accuracyon the fault location error. Particularly, PMUs return phasor quantities that assume thesystem is operating in a sinusoidal steady state, which may not be the case immediatelyafter fault initiation, when electrical transients are expected.Owing to the importance of accurate real-time applications, some metrics have beenbenchmarked to quantify measurement error in both steady and dynamic states, such astotal vector, frequency, and ROCOF errors [9]. However, there is still no guarantee toobtain satisfactory results for various applications under these criteria. Therefore, theimpact of PMU measurement errors on various applications is still under investigation byboth researchers and industry [45, 46]. GoF, as a new metric is the only criterion calculatingthe measurement error over a time frame, increasing its fidelity as a benchmark for variousapplications.In Chapter 3, inspired by the fault-location method proposed in [42], we extend thisto an impedance-based fault location for a three-phase circuit for simultaneous FDI andincorporate PMU models to the simulations. Then, we highlight the need for the GoFmetric for real-time applications, specifically with respect to fault location and show thesensitivity of fault location accuracy to GoF. For this purpose, we need to model theproblem in detail using electromagnetic transients-based models in MATLAB Simulink.We also need to model the PMUs at the ends of the line to simulate the possible errorsincorporated in synchrophasors in the during-fault period. We will show a large GoFimplies accurate fault location while a small GoF implies inaccurate fault location.1.3.3 Objective 3: Real-time Wide-area Observer-based FDI UsingSynchrophasor-based WAMSThere are various approaches for fault location in literature, including impedance-based,travelling wave-based, high frequency, and automatic learning approaches [47]. The methodwe proposed to investigate in 1.3.2 is categorized under impedance-based approaches andapplies double-end measurements for fault location. This means that the faulted line isalready detected by relays and synchrophasors are just applied for fault location. However,in some recent studies, FDI is considered as a system-wide problem [44, 48–52], in whichthe faulted line can be detected prior to the fault location on that line and independentof relays’ operation. In such cases, FDI is not just based on the measurements obtained14from one- or double-end measurements, but also on measurements obtained from variousterminals of the system for identification and location the fault. We assume that we re-ceive data at the control system from PMUs in the system, while they may not have beeninstalled at all system buses. This system-wide FDI is a potential real-time application ofWAMS.In Chapter 4, we propose a new wide-area fault location method based on the networkimpedance matrix and a linear observer. We test the proposed method on some systemsmodeled in DSATools. These models are based on positive sequence for large-scale powersystems. We assume that synthesized sychrophasors are received at the control centreor the application location from various PMUs in real time and the method calculate theresults with the same rate. The proposed method is able to identify and locate the faultsimultaneously as well as in real-time, and is not based on a hierarchical structure. Sincethis method uses a sequential set of measurements instead of two sets, it is less prone tofalse alarms and mis-detections. Moreover, the applied data to the proposed algorithmdoes not need any pre-processing.1.4 Thesis OutlineThis thesis is composed of five chapters. In the current chapter, i.e., Chapter 1, anintroduction of modern synchrophasor-based WAMS and its functionalities is provided. Inaddition, research objectives with a brief discussion about the two real-time applicationsof WAMS under study are presented. Application of modern synchrophasor-based WAMSin real-time TSA is discussed in Chapter 2. The AL approach used for the study andthe related classification and regression tools, i.e., CART and MARS, are also discussed.In Chapter 3, the effect of GoF as a new metric is studied on the credibility of faultlocation results. Chapter 4 provides a new fault location method based on observer-basedtechniques. Finally, our research conclusion and the proposed future work are provided inChapter 5.15Chapter 2Real-time Transient StabilityAssessment Using AutomaticLearning Approach2.1 BackgroundReal-time dynamic security assessment and control is a challenging issue in modern electricpower systems. Uncertainties from renewable generations, increased electricity demand,and recent blackouts have revealed the need for better tools to monitor stability. Muchresearch has been conducted in this field to enable power systems to operate closer to theirlimits and withstand the occurrence of unexpected contingencies.With the addition of new devices and technologies to modern power systems, real-timeWAMS has become more feasible. WAMS can provide real-time synchronized measurements,vulnerability assessments, self-healing, and adaptive reconfiguration of the power system.PMUs have a key role in providing real-time synchronized measurements. The samplingrate for data acquisition is in the milliseconds range. This data can be used for post-contingency dynamic vulnerability assessment, which is the preliminary step for operationof system protection schemes [41].Power system security is generally defined as the robustness of a power system tooperate at an equilibrium point in normal and perturbed conditions. Power system securityis divided into types, i.e., “static security” and “dynamic security” [3].Many of the recent studies have been conducted in the field of “dynamic security”, or,“stability” of power systems. Power system stability is known as the ability of a powersystem to remain in a stable equilibrium state at normal operating conditions, and alsoto reach an acceptable equilibrium state after being perturbed by a disturbance [5].16Figure 2.1: Different types of power system stability phenomena. Adapted from [3].Different parameters should be considered to classify the stability phenomena in powersystems such as the size of the disturbance, the time span considered for the study of powersystem stability issues, and the nature of the resultant instability [3].In Fig. 2.1 , different types of power system stability phenomena are shown. Transientstability or angular stability is a sub-category of dynamic security. Based on the afore-mentioned important parameters, transient stability is categorized as the study of anglestability for a short-period of time under large disturbances occurring in the system.The important concern in transient stability studies is to determine if the system isable to maintain its synchronous operation when it is subjected to a large disturbance. Inother words, when a large disturbance occurs in a power system, machine rotor angles aredeviated from their normal operating state and hence, the synchronism may be lost [53].Transient instability usually happens in a few seconds (less than a minute) after thedisturbance and is known as the fastest type of instability of power systems. Furthermore,as a consequence of large disturbance, many nonlinear characteristics and models shouldbe considered. Hence, differential equations along with algebraic constraints needed to beused to model the system, in addition to algebraic constraints.2.2 Approaches for TSATSA studies in power systems can be categorized into three main groups, i.e., time-domain17approaches, direct approaches, and automatic learning approaches. These approaches arediscussed in the next parts of this section.2.2.1 Time-Domain ApproachIn this approach, system robustness against a disturbance is assessed through the time-domain solutions of the following equations as [3]x˙ = f(x, y, p) (2.1)0 = g(x, y, p) (2.2)where x is the vector of dynamic states, y is the vector of static states, and p includes theparameters such as loading level, generator output, and interface flows. The vector x in(2.1) is composed of dynamic state variables, e.g., generator rotor angles, generator rotorspeeds, and field voltages, and the vector of static states y includes voltage magnitudesand angles [21].(2.1) includes differential equations used to model the elements acting dynamically inpower systems such as generators, motors, their controls, and etc., while, (2.2) is composedof algebraic equations used to model static parts of the system such as static loads.Through these equations, dynamic behaviour of the system is modelled for the during-fault period (almost 100 ms) and post-fault period. The post-fault simulation period maybe quite different considering the modelling details, but it does not usually pass beyond15 s for full detailed modelling [3].This approach has numerous advantages. For instance, it is capable of providinginformation about relevant parameters of power systems during the dynamic evolution ofthe system with time. Furthermore, if precise and detailed models are applied through theequations, it can provide accurate results about the system dynamic behaviour followinga disturbance.However, this approach cannot provide any overview about the system’s stability mar-gins, and this margin can play an important role preventive and corrective control andprotection. Moreover, this approach is not fast enough for real-time TSA of a power system[3]. Recently, new methods have been proposed to accelerate the calculation process andmake this approach suitable for real-time applications [54].2.2.2 Direct ApproachIn the direct approach, an energy function is defined to estimate the amount of criticalenergy which can be absorbed by the system before fault clearance. The hope is that thesystem will be able to restore to a stable state after fault clearance [55].18The main advantage of the energy function method is lower integration burden. Inthe time-domain approach, both during-fault and post-fault systems should be integratedfor transient TSA, while in direct methods, only the fault-on, i.e., during-fault, path needsto be integrated to determine if the initial point of the post-fault trajectory lies withinthe stability region of an equilibrium point. If the energy function value is lower than thecritical amount of the energy function, the system will be stable. This approach is alsocapable of providing a quantitative measure of the system stability margin [21].Direct methods for stability analysis include the following steps [21],1. Numerical simulation of the fault-on trajectory,2. Computation of the initial point for the post-fault system,3. Construction of the energy function for the post-fault system,4. Computation of the energy function value for the initial point of the post-fault system,5. Computation of the critical energy for the fault-on trajectory,6. If the value obtained at step 4 is smaller than the values determined at step 5, thesystem will be stable. Otherwise, it may be unstable.If the initial condition of the post-fault system is inside the stability region of thestable equilibrium point, then it can be claimed that the post-fault trajectory is stable.Therefore, the stability region is a basic concept in studies related to the direct method[21].Energy functions in general are divided into two categories, i.e., analytical and numer-ical. Analytical energy functions can be obtained only for the lossless systems [21]. Theyare path independent and no integration related to the angle can be seen in their formula-tion. In contrast, numerical energy functions include path-dependent integral terms andcan also be obtained for systems with losses, as well.To derive a numerical energy function, first, the analytical terms are constructed withconsideration of the system’s lossy components. Second, the numerical energy terms whichare composed of path-dependent integral terms are added. The result is a numerical energyfunction composed of analytic energy and path-dependent terms related to the loss of thenetwork [21].A complete numerical energy function is very complicated and composed of manydifferent terms. A simplified numerical energy function for an n-machine power systemcan be obtained as [21]Ws =n∑i=112Miω2i −n∑i=1∫ θiθSiPmidδi +n∑i=1∫ θiθSiPeidδi +n∑i=1∫ θiθSiDiωidδi (2.3)The first term represents the kinetic energy, while the second, third and fourth termsrepresent the mechanical, electrical and damping energies. We will be inspired by these19Data Set Generation Stable/unstable cases Probable ScenariosFeature Selection/Extraction Principal component analysis (PCA) Cross-correlation (CC) Physical interpretable indicesDeveloping a Classification/Regression Tool Artificial Neural-Networks (ANN) Decision tress (DTs) Support vector machine (SVM)Figure 2.2: General automatic learning framework for TSA.equations to obtain useful indices in in the next sections.A complete energy function associated with power flow equations and rotor circuitequations includes more path-dependent integrals and is more complicated to calculate.As a result, it will be difficult to apply an energy function approach for real-time TSA inlarge scale systems. Indeed, the available studies are based on simplified energy functionapproaches [56].2.2.3 Automatic Learning ApproachIn this approach, first, the dynamic response of the system in post-disturbance period issimulated in the time-domain. The output of the process is stored in a data set. Thesesimulation cases may be stable or unstable. Next, dimension reduction techniques areused. Different feature selection and feature extraction techniques may be applied toobtain better discrimination capability between different system contingencies. Thesefeatures may be physical features derived from the simulations or may be obtained basedon mathematical theories such as principal component analysis, cross correlation, etc. Theresultant feature sets extracted from the obtained simulation cases will be saved in a database to be used in the following step.The final step is the application of classification or regression tools. These classificationtools are “trained” to provide the capability of automatic prediction of unstable cases,immediately after fault clearance, before the system reaches an unstable stage. This earlyprediction provides the feasibility for corrective control actions. Furthermore, regressionmethods may be applied. Regression methods are applied when a continuous variable isconsidered as the output of the simulation instead of discrete classes. In other words, theidea is to “train” a model which maps the input features to a continuous output variable.The framework including all of the above steps is shown in Fig. 2.2.202.3 Real-time TSA Using AL ApproachAs discussed in Section 2.2, the time-domain approaches, direct approaches, and finally,automatic learning approaches are three main approaches of TSA in the respective lit-erature, while the latter has recently begun garnering more attention in the literature[22–25, 29, 30, 34].A comprehensive database of the Hydro-Que´bec power system is constructed in [25].Indices or features named “wide-area severity indices” and in addition to voltage criteriaindices in the time-domain are applied here. These features result in a higher discrim-inability for different cases. In the next step, a DT is applied as the classifier. Somefeatures related to center of inertia are applied in [34]. A random forest (RF) classifierwhich is an ensemble of DTs is used. The correlation of features is also considered to elimi-nate repetitive features which may increase the processing burden during the classificationstage without accuracy increase. In [29], it is shown that RF can operate efficiently in thepresence of small changes in the network topology and various dynamics of the powersystem.Important post-contingency security issues, including voltage magnitude violation,thermal limit violation, voltage stability, and transient stability are considered simul-taneously in [22] . The focus of [22] is not just on transient stability; both static anddynamic security are considered using different indices. For example, thermal limit vi-olation and voltage magnitude violation are considered in static security, while voltagestability and transient stability are in the category of dynamic security assessment. Inthis paper, development of a prediction tool to obtain the severity of a contingency isalso under consideration. Hence, RTs have been applied in this study. RTs are appliedto indicate the number of overloaded lines and buses with voltage magnitude violationsfollowing a contingency. A higher number shows a more severe situation. In [23], thedata-mining section is under focus. Application of multiple DTs while also boosting thetrees leads to improved reliability for the proposed method. The algorithm is updatedevery day to incorporate the probable reconfigurations of the system. This paper coversthe most significant disadvantage of model-free methods, i.e., compatibility with networkreconfigurations after a while.ensemble learning machine (ELM) is used in [30] and [24] to decrease the time oftraining and decision making and prepare the schemes for online applications. A singleELM is applied in [30] while an ensemble of ELM prediction tools is assessed in [24].Selection of an appropriate tool for TSA has always challenged researchers, and thus,various classifiers have been compared in [57]. These classifiers have been divided into twomajor groups; first, black-box predictors such as RFs, neural networks, and support vectormachines; and second, comprehensible, transparent and interpretable predictors such as21DT and Fuzzy-DT. Black-box predictors can reach to the high accuracy rates necessary formodern power grids, while the accuracy rate of interpretable classification tools expressingtransparent rules decreases. However, many system operators and managers in controlcentres prefer transparent and interpretable models, because of their ease in decision-making and forensic uses [57]. On the other hand, using discriminative indices and featurescan increase the accuracy rate of any predictor, including transparent ones.2.4 Real-time TSA on the BCH Power SystemIn this chapter, we have selected an AL approach for real-time TSA. We apply the proposedmethod to the WECC power system. This system includes 17724 buses, 22005 transmissionlines and transformers, and a 170000 MW generation capacity [58]. Next, we will focus onthe transient stability issues of the BCH power system as a part of the WECC system. Thispower system is connected to BPA from the southern part and exchanges up to 2500 MWpower with the United States Utilities. The interconnection is called the Northern Intertie.Therefore, any severe fault leading to a disconnection of the Northern Intertie can leadto transient stability issues in the power system of BCH. The diagram of the BCH powersystem is shown in Section 2.7. Here, we have considered that PMUs are installed at themajor power plants and critical substations of BCH, BPA and also, the Alberta ElectricitySystem (AES) as the second power system interconnected to the BCH system. The obtainedsynchrophasor measurements can be sent to the control centres of the systems througha communication infrastructure. A portion of these assets are now available at the BCHsystem.To conduct this research, we have defined and simulated 48 different credible contingen-cies and have considered 3 different topologies of the system affecting the Northern Intertiepower transfer. In total, we have created a data set composed of 9474 cases, about 12%of which are unstable. This huge data set certifies the reliability of the obtained models.We propose and select different indices for wide-area dynamic security assessment.The features are divided into three main categories; first, features based on basic timesynchronized measurements obtained from WAMS; second, indices obtained through simplecalculations considering the COI; and third, the kinetic and potential energy functions formajor power plants calculated as the energy function-based indices [59].The indices obtained for the first five cycles after the fault clearance are consideredto “train” five different DTs as the classification tools. Since DTs are interpretable clas-sification tools, their most discriminative features can also be selected after developingthe models [60]. Moreover, the indices of the first five cycles after the fault clearance canbe used to develop five MARS models to obtain close form equations of transient stabilitymargins based on the proposed features.222.5 Wide-area Dynamic Security IndicesMany techniques have been explored in the context of TSA, including energy functions[21, 61, 62]. Recent work has revealed prediction tools to be promising for online TSA[22–25, 29, 34, 63]. With regard to this, numerous metrics, features, and indices have beenproposed to be applied to classifiers, including pure real-time measurements obtained fromthe power system and more complicated indices or features resulting from computationsinvolving such measurements [22–25, 29, 34, 63]. In this section, we describe TSA indicesused in our work, which can be categorized as (i) pure measurement indices, (ii) COI-referred indices, and (iii) energy function-based indices.2.5.1 Basic Measurement IndicesThe most basic indices used in literature are direct time-synchronized measurements of,e.g., generator rotor angles, rotor speeds, bus voltage phase angles and magnitude, andcurrent flows [22–24, 34]. Other basic indices are obtained from simple calculations in-volving time-synchronized measurements, e.g., active and reactive generation, active andreactive load, and voltage drop across transmission lines [22–24].In this work, based on the practical availability of measurements from the WAMS, weselect several synchronized measurements from the pre- and post-fault systems and usethem as TSA indices. In any real power system, major power plants with the highestgeneration are candidates for PMU installation and real-time monitoring [64]. So, withregard to pre-fault measurements, we use the total active power generation of major powerplants in each region, the active power generation of each individual major power plant,the voltage phase angles of buses connected to major power plants and those of interfacebuses, the active power flow across inter-ties connecting the regions, and the active powerflow along critical long transmission lines. Under post-fault conditions, we use only time-synchronized measurements of voltage phase angles at the buses connected to major powerplants and at buses that interface between two regions.2.5.2 Centre-of-Inertia-Referred IndicesSince deviation from the COI angle can indicate system stress, below, we introduce COI-referred response signals. Consider an interconnected power system in which the monitoredregion hasK−1 neighbouring regions. Assume we monitor region k which containsN busesand NG major synchronous generators to be monitored. Denote, by Pmi(t) and Pei(t),the mechanical power input [p.u.] and the electrical power output [p.u.], respectively, ofgenerator i at time t. Let δi(t) and ωi(t) denote the rotor angle [rad] and speed [rad/s],respectively, of generator i at time t; and let Mi denote the inertia constant [s2/rad] of23generator i. Then, the COI angle of the monitored region is defined as [65]δCOI,k(t) :=∑NGi=1Miδi(t)Mtotal, (2.4)where Mtotal =∑NGi=1Mi; and the COI rotor speed is analogously defined as [65]ωCOI,k(t) :=∑NGi=1Miωi(t)Mtotal. (2.5)Using (2.4) and (2.5), we define the COI-referenced angle and speed of generator i inregion k as δ˜i(t) = δi(t) − δCOI,k(t) and ω˜i(t) = ωi(t) − ωCOI,k(t), respectively. Withthese notations in place, we next describe COI-referred indices that are applied to our TSAalgorithm.First, we compare the pre- and post-fault conditions and propose the following centredeviation index:∆δ(t) =δCOI,k(t)− δpre-faultCOI,kδpre-faultCOI,k, (2.6)where δpre-faultCOI,k represents the pre-fault COI angle of the monitored region. Furthermore,it has been shown that, based on energy function concepts, products involving post-faultpower mismatches and angle and speed deviations are good indicators of system stress [66].Thus, we also apply the following indices [25, 66]:αk(t) =NG∑i=1fi(t)ω˜i(t), (2.7)where fi(t) = Pmi(t)− Pei(t)− MiMtotalPCOI,k(t), with PCOI,k(t) =∑NGi=1(Pmi(t)− Pei(t)),βk(t) =NG∑i=1fi(t)δ˜i(t), (2.8)andγk(t) =NG∑i=1ω˜i(t)(δ˜i(t)− δ˜pre-faulti ), (2.9)where δ˜pre-faulti denotes the pre-fault COI-referenced angle of generator i. In addition tothe indices described in (2.6)–(2.9), we modify γk(t) in (2.9) to reflect the behaviour ofthe COI angle and speed, and include the following index:σk(t) = ωCOI,k(t)(δCOI,k(t)− δpre-faultCOI,k ). (2.10)24Not captured by indices in (2.6)–(2.10), we may wish to penalize generators whoserotor angles deviate further from the COI angle. Thus, aimed at an index that is sensitiveto loss of synchronism, we make use of the integral square generator angle index, which isexpressed as [67]J =1T∫ T01MtotalNG∑i=1Mi(δi(t)− δCOI,k(t))2dt, (2.11)where T represents the observation time window. Note that, in (2.11), the index is nor-malized with respect to both Mtotal as well as T . Thus, in the case of disconnection ofany generator, its contribution to both Mtotal and δCOI,k(t) must be removed from (2.11).To assess transient stability, it is not only important to consider characteristics of themonitored region, but to also account for interactions between the monitored region andits neighbours. To this end, we consider the centre of angle stability index, expressedas [59]∆δk,j(t) =∣∣∣∣∣δCOI,k(t)− δCOI,j(t)δpre-faultCOI,k + δpre-faultCOI,j∣∣∣∣∣× 100, (2.12)which can be interpreted as the angle deviation between regions k and j where k 6= j.Rotor angles and speeds are very important indicators in power system stability studies[5]. The differences of bus angle and frequency of two buses are combined in [68] tocalculate the rms-coherency between those buses in order to partition a system. Inspiredby this, to assess the coherency between a region and its neighbouring regions following afault, we consider COI angle and rotor speed of region k and j as a criterion and calculatethe rms-coherency between them, as follows:ψk,j =√1T∫ T0((δCOI,k(t)− δCOI,j(t))2+(ωCOI,k(t)− ωCOI,j(t))2)dt (2.13)where T represents the observation time window.To summarize, we consider COI-referred indices (2.6)–(2.13) in real-time TSA.2.5.3 Energy Function-Based IndicesAssuming the remainder of the system can be modelled as an infinite bus, the dynamicsof generator i can be expressed as [65]Mid2δidt2= Pmi − Pei. (2.14)25From (2.14), we recover the kinetic energy function Vi of generator i as [65]Vi(t) =12Mi(dδidt)2. (2.15)Similarly, the potential energy function Ui of generator i can be expressed asUi(t) =∫ t0(−Pmi + Pei(τ)) dδidτdτ. (2.16)In Section 2.7, our simulations assume that Pei(τ) measurements can be obtained fromPMUs, from which Ui(t) can be numerically computed.Based on (2.16), as an indicator of transient stability in region k, we propose to applythe following total potential energy index:Utotal,k(t) =NG∑i=1Ui(t), (2.17)where Ui(t) is computed numerically.Similarly, from (2.15), we obtain the total kinetic energy in region k asVtotal,k(t) =NG∑i=1Vi(t), (2.18)which is used as an index in our study. As discussed in Section 2.5.2, the interactionsbetween the monitored region and its neighbours are also critical for transient stabilityassessment. To this end, we propose the kinetic energy separation index, formulated as∆Vk,j(t) = |Vtotal,k(t)− Vtotal,j(t)|, (2.19)to be applied to to the proposed TSA algorithms.In summary, we consider energy-based indices in (2.17)–(2.19) in training CART andMARS models for real-time TSA.2.6 Data-Mining FundamentalsDTs, as the most important transparent prediction tool, were first used in security assess-ment in [69]. This classification tool has been also used to obtain the critical clearing timeof contingency [70], voltage security assessment [71] and controlled separation studies [72].MARS is a highly automated tool for regression analysis. In RTs, a piecewise constantapproach is used for regression and; therefore, there are jumps and discontinuities. But,these constants are replaced with piecewise linear continuous approaches in MARS. Con-26sequently, smooth curves and surfaces are obtained in regression which can be stated bysimple, closed form and simple equations. MARS algorithm is only used as a regressionmodel. It automatically chooses the variables to be used in an optimal manner. Throughthe algorithm, the interactions between the variables are detected and the model is self-tested to protect against overfitting [32].In this section, we describe CART as the classification method and MARS as the regres-sion method used in the chapter. In order to build and test the models in these methods,S time-domain simulation cases are obtained by combining different contingencies andoperating conditions. For each case i = 1, ..., S, we consider R features or wide-area dy-namic security indices (WADSIs). As described in section 2.5, in the context of this work,WADSIs include pure measurement, COI-referred, and simplified energy function-based in-dices. Collect feature values for each simulation case in the S×R matrix X = [xij ], wherexij represents the value of feature j for case i.2.6.1 CARTIn order to build the CART, denote the target variable for case i by yi, which takes the valueof either 0 or 1, representing unstable or stable case, respectively. With these notations inplace, we represent the CART data set by D = {(eTi X, yi), i = 1, 2, ..., S}, where ei denotesa column vector of all zeros except with the ith entry equal to 1.To build a classification or regression model, the data set D is divided into two subsets.First, the training set, denoted by T consisting of ST simulation cases, is used to buildthe CART. The remainder partition D\T , called the test set, is used to ensure the validityof the CART built using T . Our objective is to construct a tree from the training set andsubsequently, use it for prediction of new unseen contingency simulation cases. Such atree is constructed in three stages: (i) tree growing, (ii) tree pruning, and (iii) optimaltree selection. These stages are summarized below. Interested readers may refer to [31]for more details.2.6.1.1 Tree GrowingBeginning from a root node, a feature or index, called the “splitting feature”, is selectedand a value for that feature, called the “splitting value”, is chosen to divide the traininginto two subsets. The result is two child nodes, either internal or terminal. An internalchild node becomes a parent node, at which another splitting feature and correspondingsplitting value are selected. This sequential splitting procedure eventually leads to someterminal nodes called “leaf nodes”. Each leaf node is classified as either stable or unstable.A path from the root node to a leaf node leads to a decision region in the feature set relatedto that terminal node. This region can be characterized by a sequence of yes/no questions.27In this manner, the decision tree is simple and transparent to interpret. The root node inthe tree structure includes the entire training set, but just a portion of cases would passthrough the lower nodes [31]. The subset of training data corresponding to a parent nodeis partitioned for the two child nodes to reduce the number of misclassified cases at eachchild node.A node is considered “pure” when majority of cases belonging to that node are asso-ciated with one class label, either stable or unstable. The partitioning process describedabove should lead to increase in purity of the nodes’ corresponding data. To accomplishthis aim, a node impurity function is applied, discussed below.Node Impurity Function: The node impurity function is defined in such a way to bemaximized if the cases belonging to a node include an equal number of stable and unstablecases, making it difficult to label the node with any confidence. Conversely, the impurityfunction is minimized if all cases belonging to the node are either stable or unstable, sothat the node can be labelled as “stable” or “unstable” with high accuracy [73]. Basedon the intuition above, we describe two possibilities for the node impurity function, oneof which is based on the Gini diversity index and the other on the Entropy function.Recall that, in the context of this work, there are only two classes: stable (1) andunstable (0). Denote the probability of stable cases corresponding to node z as p(z) =Pr(1|z). Then, the impurity function based on the Gini diversity index at node z is [73]i(z) = 2p(z)(1− p(z)). (2.20)On the other hand, the impurity function for node z, based on the Entropy function fortwo classes, is [73]i(z) = −p(z) log p(z)− (1− p(z)) log(1− p(z)). (2.21)Indeed, both candidate impurity functions in (2.20)–(2.21) are maximized with p = 0.5and minimized with p = 0 or 1, as expected, which is shown in Fog. 2.3.2.6.1.2 Tree PruningIntuitively, the DT’s prediction accuracy rate should increase as the number of the nodesincreases. However, if the growing process is continued until no two simulation casesbelong to the same terminal node, the tree would likely be overfitting the training set,and consequently, would not be a good classifier for the test set. As the tree grows, aftersome degree of complexity, the tree’s misclassification rate may increase. Thus, the treeobtained from the process described in section 2.6.1.1 must be pruned.We begin by defining variables that are relevant to the pruning process. Let Z represent280 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 100.10.20.30.40.5pvalue EntropyGini indexFigure 2.3: Plot of node impurity functions for two-class classification, as a functionof the probability p(z). Entropy-based impurity function is scaled to passthrough (0.5,0.5).the classification tree obtained from the tree-growing process described in section 2.6.1.1.A node z′ is a descendent of node z if there exists a connected path down the tree leadingfrom z to z′. A branch Zz of Z consists of the node z and all descendants of z, and branchZz has z as its root node. Pruning a branch Zz from a tree Z consists of deleting, fromZ, all descendants of z, while retaining node z, which becomes a leaf node of the prunedtree. Suppose Z ′ is obtained from Z by successively pruning off branches, then Z ′ Z isa pruned subtree of Z. The pruning process identifies branches that contribute least tothe model accuracy rate and eliminates them. In order to do so, the misclassification rateneeds to be estimated. One of the most prevalent estimation methods of misclassificationrate is the “resubstitution estimate”. At node z, the resubstitution estimate, denoted byr(z), for a two-class data set, is calculated as [73]r(z) = 1−max (p(z), 1− p(z)) = min (p(z), 1− p(z)) . (2.22)where p(z) = Pr(1|z).Now, let Z˜ = {z1, z2, ..., zL} denote the set of L(Z) terminal nodes of classificationtree Z and pi(zl) the proportion of all simulation cases that belong to node zl. Based on(2.22), the resubstitution estimate of the misclassification rate for Z is expressed asR(Z) =L(Z)∑l=1r(zl)pi(zl). (2.23)We apply the method of minimum cost-complexity pruning to create a set of subtreesfrom the original tree. In this method, the resubstitution estimate in (2.23) is modified to29include a complexity parameter, α ≥ 0. Then, for any subtree Z ′ ≺ Z the cost-complexitypruning rate isRα(Z′) = R(Z ′) + αL(Z ′). (2.24)In (2.24), αL(Z ′) can be interpreted as a penalty term for the size of subtree Z ′, which isnot considered in R(Z ′).For each value of α, the subtree Z∗(α) of Z that minimizes Rα(Z ′) is selected as follows:Z∗(α) = arg minZ′ZRα(Z′). (2.25)The value of α determines the size of the tree. If α is small, then the penalty termin (2.24) would be small, causing R(Z ′) to dominate. In turn, the size of the minimizingsubtree Z∗(α) would be large. Conversely, for sufficiently large values of α, Z∗(α) wouldconsist of the root node z only. It is worth noting that although α is defined on theinterval [0,∞), the number of subtrees of Z is finite. Thus, for a sequence of complexityparameters 0 < α0 < α1 < ... < αM , nested subtrees are obtained asZ Z∗(α0) Z∗(α1) ... Z∗(αM ). (2.26)2.6.1.3 Optimal Tree SelectionFrom the set of subtrees generated in the pruning procedure, we select the subtree withthe smallest misclassification rate. An independent test set is used to estimate the errorrates of the subtrees Z∗(α0), ..., Z∗(αM ). The subtree with minimum misclassification rateis selected as the final CART model [73].In the testing process, each of the simulation cases in the test set traverses each prunedtree to be classified into one of the classes. Since the class label for each observation ispredetermined, the misclassification rate for subtree Z∗(αm), corresponding to complexityparameter α = αm, can be expressed asRts(Z∗(αm)) = R(Z∗(αm)). (2.27)In other words, the resubstitution estimate calculated for the independent test set is usedas the misclassification rate. In the case of identical costs of misclassification for all classes,Rts(Z∗(αm)) is simply calculated as the proportion of misclassified observations to all thecases in Z∗(αm). The obtained misclassification rates are then used to select the bestpruned tree, Z∗, asZ∗ = arg minm=0,...,MRts(Z∗(αm)). (2.28)302.6.2 MARSSimilar to section 2.6.1, we consider the feature variable w = [w1, ..., wR]T , all of which aredescribed in section 2.5, and S simulation cases. Recall that the S×R matrix X containsthe values of each feature j = 1, ..., R, associated with each case i = 1, ..., S.Instead of considering a discrete output of 0 or 1, as for CART in section 2.6.1, here,we use MARS to map the feature variable w to a continuous quantity, which is indicativeof system stability, particularly we use the angle-based TSI defined asTSI :=360−∆δmax360 + ∆δmax× 100 = f(w) + , (2.29)where f(w) : RR → R. The goal of MARS is to fit to a function of features f(w) thisoutput observation while minimizing error . In (2.29), ∆δmax denotes the maximumangle separation between any two single generators in a region. TSI > 0 and TSI < 0represent stable and unstable conditions, respectively [74]. For ease of notation, let yidenote the TSI metric computed in simulation case i.The MARS method results in a model of the formf(w) =Q∑q=0βqhq(w), (2.30)where hq(w) denotes a basis function (BF) or a product of two or more basis functions (BFs)(we elaborate on these later), βq represents a constant multiplicative coefficient, Q denotesthe number of additive terms in the final model, and for case i, w = [xi1, xi2, ..., xiR].Similar to the CART algorithm, the MARS method first creates an overfit model and thenprunes it back. This process is composed of two stages—forward and backward, which wesummarize below [73]. Interested readers may refer to [32] for details.2.6.2.1 Forward ProcedureThe MARS method first collects pairs of BFs for each feature wj . The BF may be a constantvalue or a hinge function. A hinge function takes the form of max(0, wj−xij) for wj ≥ xij ,and max(0, xij − wj) for wj < xij , where the constant xij is called a “knot” and is thevalue of feature j for a particular simulation case i in the training set. The pair of BFs forthe jth feature is expressed as the setCj = {max(0, wj − υj),max(0, υj − wj)}, (2.31)where υj ∈ {x1j , x2j , ..., xST j}. Given R features or WADSIs, 2R BFs in the set C =C1 ∪ C2 ∪ ... ∪ CR are defined for the training set.31With the set of BFs C in place, let R = {1, ..., R} and let the model obtained at the endof the forward procedure be f(w) =∑Qmaxq=0 βqhq(w) with Qmax denoting the total numberof the added terms in the model. The forward procedure is iterative and begins with aconstant function f0(w) = h0(w) = c which is the mean of the training cases’ output tominimize the error in (2.29).At each iteration q = 1, ..., Qmax, for each j ∈ R and simulation case r = 1, ..., ST inthe traning set, we formhj,rq = βj,rq1 hq−1max(0, wj − xrj) + βj,rq2 hq−1max(0, xrj − wj). (2.32)from which we formf j,rq (w) = fq−1(w) + hj,rq (w). (2.33)In 2.32, the coefficients βj,rq1 and βj,rq2 are calculated based on least-squares criterion asβj,rq1 , βj,rq2 = arg min(ST∑i=1∣∣∣∣yi − f j,rq (w)∣∣w=[xi1,...,xiR]T∣∣∣∣2). (2.34)Next, from the pool of candidate functions f j,rq , we would choose the one that minimizesleast-squares error as follows:j?, r? = arg minj∈R,r=1,...,ST(ST∑i=1∣∣∣∣yi − f j,rq (w)∣∣w=[xi1,...,xiR]T∣∣∣∣2), (2.35)Finally, we set fq(w) = fq−1(w) + hj?,r?q (w) and remove feature j? from consideration inthe next iteration, i.e., set R = R \ j? and C = C \ Cj?.The iterative forward procedure is continued until the change in the residual squarederror isST∑i=1∣∣∣∣yi − fq(w)∣∣w=[xi1,...,xiR]T∣∣∣∣2 < 0, (2.36)where 0 > 0 is predetermined, or until the preset maximum iterations Qmax is reached.Algorithm 1 summarizes the pseudo code of the forward procedure [32].This forward procedure is hierarchical, i.e., the terms from the collection C are mul-tiplied by the terms already involved in the model to build a multiway product. In thisway, we can identify higher-order interactions between the features [73].32Algorithm 1Input: C, R = {1, ..., R}, and the output of simulation cases yi, i = 1, 2, ..., ST .Output: The overfit model of f(w) =∑Qmaxq=0 βqhq(w) and Qmax.1: Initialize. Set the initial function f0(w) = c =: h0(w) and set the counter q = 0.2: while∑STi=1∣∣∣∣yi − f j,rq (w)∣∣w=[xi1,...,xiR]T∣∣∣∣2 > do3: Set. q ← q + 14: for j ∈ R do5: for r = 1 to ST do6: form. hj,rq (w) = βj,rq1 hq−1max(0, wj − xrj) + βj,rq2 hq−1max(0, xrj − wj).7: form.f j,rq (w) = fq−1(w) + hj,rq (w).8: Calculate. βj,rq1 , βj,rq2 = arg min(∑STi=1∣∣∣∣yi − f j,rq (w)∣∣w=[xi1,...,xiR]T∣∣∣∣2).9: Set. hj,rq (w)=βj,rq1 hq−1max(0, wj − xrj) + βj,rq2 hq−1max(0, xrj − wj).10: end for11: end for12: Select. j?, r? = arg minj∈R,r=1,...,ST(∑STi=1∣∣∣∣yi − f j,rq (w)∣∣w=[xi1,...,xiR]T∣∣∣∣2)13: Set. fq(w) = fq−1(w) + hj?,r?q (w).14: Set. C ← C \ Cj?15: Set. R ← R \ j?16: Set. Qmax = q17: Set. hq = fq(w)18: end while2.6.2.2 Backward ProcedureThe large model obtained at the end of the forward process typically overfits the data set,so, a backward detection procedure is necessary. In this procedure, terms in fQmax areremoved to successively improve the accuracy of the final resulting MARS model [73].Algorithm 2 shows the pseudo code of the backward procedure. In this iterativeprocedure, we begin with fQmax(w) obtained from the forward procedure and initializefˆλ+1(w) = fQmax(w). At each iteration, λ = Qmax, ..., 1, we form fˆ`λ(w) by removing the`th term from fˆλ+1(w), where ` = 0, ..., λ. Among the λ + 1 functions, we select fˆ`?λ (w),where`? = arg min`=0,...,λ(ST∑i=1∣∣∣∣yi − fˆ`(w)∣∣w=[xi1,...,xiR]T∣∣∣∣2). (2.37)At the end of iteration λ, set fˆλ+1(w) = fˆ`?λ (w) and repeat the process above. Finally,from the pool of candidate fˆλ functions, λ = 1, ..., Qmax, the final MARS model is chosen33Algorithm 2Input: The overfit model f(w) =∑Qmaxq=0 βqhq(w) and the output of simulation cases yi,i = 1, 2, ..., ST .Output: fˆλ?(w) =∑Qq=0 βqhq(w)1: Initialize. fˆλ+1(w) = fˆQmax =∑Qmaxq=0 βqhq(w).2: for λ = Qmax to 1 do3: for l = 0 to λ do4: Form. fˆ lλ(w) by removing lth term of fˆλ+1(w)5: end for6: Select. l? = arg minl=0,...,λ(∑STi=1∣∣∣∣yi − fˆl(w)∣∣w=[xi1,...,xiR]T∣∣∣∣2).7: Set. fˆλ+1(w) = fˆl?λ (w)8: end for9: Select. λ? = arg minλ=1,...,Qmax(GCV(λ) =∑STi=1(yi−fˆλ(w))2(1−M(λ)/ST )2)as fˆλ?(w) whereλ? = arg minλ=1,...,Qmax(GCV(λ) =∑STi=1(yi − fˆλ(w))2(1−M(λ)/ST )2), (2.38)M(λ) can be interpreted as a penalty factor for greater model complexity. Specifically,suppose the model includes F linearly independent BFs, then M(λ) is obtained as [73]M(λ) = F + c1(Qmax − 1), (2.39)where the penalty coefficient c1 is a user-defined constant. Interested reader may refer to[32] for more details.2.7 TrainingIn this work, we consider the WECC bulk transmission system model mentioned previously,from summer 2013, which includes 17724 buses, 22005 transmission lines and transformers,and 170000 MW generation capacity [58]. As shown in Fig. 2.4, we monitor the BCH powersystem, which is connected to two neighbouring regions: the AES and the BPA systems. Asmentioned, the interconnection between BCH and BPA is called the Northern Intertie andconsists of a 500-kV double-circuit and a 230-kV transmission line, transferring over morethan 2500 MW of active power exchange between BCH and BPA. BCH is also connectedto AES on the east, via two interties, a 500-kV and a 230-kV transmission line, enablingthem to exchange up to 800 MW active power.According to the North American Synchrophasor Initiative (NASPI), PMUs should34Figure 2.4: Diagram of BCH and its interconnections to BPA and AES [4].be installed at major power plants with high generation capacity as well as along criticalcorridors and bottlenecks [64]. The present study assumes that PMU measurements areavailable at major power plants with high generation capacity in the three power systems(5 in BCH, 15 in BPA, and 7 in AES), the interties between these three systems, andalso a major corridor transferring power southward from northern BC. Currently, PMUsare installed at those 5 power plants in BCH, i.e., GMS, PCN, KMO, REV, and MCA.Consequently, power flows, angles, voltage magnitudes and frequencies at these locationslisted above are accessible at the control centres.To cover a wide range of operating conditions, 9474 time-domain simulation cases,we include 3 possible topologies, 200 WECC power-flow base cases from the summer of2013, and 48 credible N − k (k = 1, ..., 4) contingencies. We conduct these simulations inTSAT [74] and find 1098 cases to be unstable. The simulations assume that a contingencyoccurs at t = 100 ms and lasts for 3 cycles, i.e., 50 ms, before it is cleared. Based onthe simulations, the WADSIs or features are obtained for different contingencies in variouspower flow patterns. The synthetic time-synchronized measurements in the present studyare obtained from these time domain simulations.Selection of the appropriate time instant for data acquisition is critical in TSA. Differenttime spans after fault clearance are considered for synchronized measurements sampling inthe literature. For example, in [29] and [57], 150 ms and 300 ms after the fault clearancetime and in [75] are considered, as well as 5–6 cycles after the clearance time in [76].35Figure 2.5: A transiently unstable case following a fault occurrence.Finally, in [77] and [63], it is considered that the PMU measurements are updated at eachcycle (16.67 ms).In this work, synthetic synchronized measurements are obtained in each of the 5 cyclesafter the fault clearance, which is possible with P-class PMUs [10]. Consequently, the timeconsidered for sampling will be less than 100 ms.The synchronized measurements of the first 5 cycles after the fault clearance time areused to train 5 DTs via the procedure described in section 2.6.1. 70% of the data setis used for training the DTs. The remaining 30% is applied for testing. The syntheticmeasurements are also used to obtain the transient stability thresholds of the system inthe first 5 cycles of the post-fault instants, using the MARS algorithm described in section2.6.2. We make use of existing commercial software Special Predictive Modeler (SPM)[78, 79] to build the CART and MARS models from the training set data.2.8 ResultsIn this section, we first show some of the characteristics of the WADSIs obtained from time-domain simulations are shown. Next, the results obtained from CART and MARS methodswill be presented and discussed [80].2.8.1 Wide-area Dynamic Security Indices ResultsFigure 2.5 shows an unstable simulation case. In this figure, it can be observed that therotor angles are diverging and divided into three groups after the fault clearance. Each ofthese three groups belong to one of the regions under the study, i.e., BPA, BCH and AES. Itcan be inferred from this figure, that the coherent regions have been considered correctlyin our study.The simulation results for some of the WADSIs are shown in the Fig. 2.6. In these36(a) α(t)BCH ,(b) ∆δ(t)BCH .Figure 2.6: Characteristic of two COI-referred indices following different faults inthree different power flow patterns.figures, three different power flow cases have been selected, two of them for the activepower transfer from BCH to BPA and one for the active power transfer from BPA to BCH.Some separate contingencies are also selected which can lead to stable and unstable cases.In Fig. 2.6, the proposed WADSIs based on the COI are plotted for BCH and the during-fault and the post-fault periods have been assigned to further bold the indices. Here, thecurves are independent of the power flow cases or the operating points and are capable ofdiscriminating between the stable and unstable cases.2.8.2 CART ResultsAs mentioned in Section 2.7, five CART models are trained using synthetic measurementsof each of the first five cycles after fault clearance, along with pre-fault measurementsdescribed in Section 2.5.1. The CART model obtained from the fourth-cycle measurementsis shown in Fig. 2.7, which is typical of resulting trees from the other four cycles. In37Table 2.1, we summarize the selected features for the top three levels in each of the fiveCART models. The root node in all five models is the index αBCH . The index ∆VBCH,BPAemerges in the second level in three of the five models. Moreover, the indices σBCH , γBCH ,∆δBCH,BPA, and Vtotal,BPA are selected in the second level. All the indices selected in thefirst three levels are either COI-referred or energy function-based ones. This observationindicates the importance of these indices in TSA. Moreover, indices related to the BCHsystem play an essential role in the trees’ structures.We report accuracy rates of the resulting CART models in Table 2.2, in which the CARTmodel built using the indices from each cycle is tested against the data obtained at allfive cycles. In other words, the diagonal entries in Table 2.2 indicate the accuracy of eachCART model when fed by data from the expected sampling time, while the off-diagonalentries indicate the model accuracy when data from an earlier or later time are used. Foraccuracy rates other than those shown on the diagonal, the CART models are tested usingthe whole data set of the other cycle (instead of only the test set). We note that, inTable 2.2, the accuracy rates of the corresponding cycles (diagonal elements) are about99%. When the trained CART model is tested against data obtained from one cycle beforeor after the corresponding cycle, the accuracy rates are around 90% and in some casesclose to 95%. These off-diagonal entries are relevant for cases in which fault-clearancetime is not detected properly by the PMUs and the data is fed into the wrong tree. Also,considering the off-diagonal entries in Table 2.2, the decrease in accuracy is particularlynoticeable when the training data belongs to cycles 4 and 5, and the testing data belongsto the earlier cycles 1 and 2. Thus, our results suggest that the obtained models are moresensitive when they are trained using data obtained from the later cycles and tested withdata from earlier cycles.The idea of training a sequential set of CART models over the time frame of five cyclescan also be interpreted as a ensemble of DTs distributed over time. Therefore, the overalloutput of the models can be based on voting result of all the five trained CART models,which increases the reliability of the final result compared to a single CART model.The computational times required to obtain the CART models are reported in Table 2.3.As evident from the table, the time required to obtain a CART model using the data setfrom the same cycle (diagonal entries) is around 12 s, and it increases to around 19 s in thecase of using data from different cycles (off-diagonal entries). This increase is the resultof the increase in the size of training and testing data sets. Note that the commercialsoftware SPM provides the total time required to train and test the models.38αBCHClass Cases Percent 0 764 11.5% 1 5902 88.5% Total Cases = 6666γBCHClass Cases Percent 0 21 11.5% 1 5478 88.5% Total Cases = 5499γBPAClass Cases Percent 0 21 3.9% 1 515 96.1% Total Cases = 536Class = 1Class Cases Percent 0 0 0% 1 4963 100% Total Cases = 4963Class = 1Class Cases Percent 0 0 0.0% 1 494 100% Total Cases = 494Class = 0Class Cases Percent 0 21 50% 1 21 50% Total Cases = 42γBCH 0.0 γBCH>0.0αBCH 0.06γBPA 0.0 γBPA>0.0ΔVBCH, BPAClass Cases Percent 0 743 63.7% 1 424 36.3% Total Cases = 1167αBCHClass Cases Percent 0 591 92.2% 1 50 7.8% Total Cases = 641βBPA Class Cases Percent 0 152 28.9% 1 374 71.1% Total Cases = 526αBPAClass Cases Percent 0 26 39.4% 1 40 60.6% Total Cases = 66Class = 0Class Cases Percent 0 565 98.3% 1 10 1.7% Total Cases = 575ΔVBCH, BPA 0.05ΔVBCH, BPA>0.05αBCH 0.08αBCH>0.08αBCH>0.06βAESClass Cases Percent 0 19 5.8% 1 308 94.2% Total Cases = 327Pg, BCHClass Cases Percent 0 133 66.8% 1 66 33.2% Total Cases = 199βBPA 0.27βBPA>0.27Pg, BCH 10860.49Pg, BCH >10860.49Class = 1Class Cases Percent 0 0 0% 1 48 100% Total Cases = 48Class = 0Class Cases Percent 0 133 88.1% 1 18 11.9% Total Cases = 151βAES 0.0Class = 0Class Cases Percent 0 15 83.3% 1 3 16.7% Total Cases = 18αBPAClass Cases Percent 0 4 1.3% 1 305 98.7% Total Cases = 309βAES >0.0Class = 1Class Cases Percent 0 0 0.0% 1 304 100% Total Cases = 304Class = 0Class Cases Percent 0 4 80% 1 1 20% Total Cases = 5αBPA 0.05 αBPA >0.05αBPA>0.03Class = 0Class Cases Percent 0 23 100% 1 0.0 0.0% Total Cases = 23αBPA 0.03PG5BPA Class Cases Percent 0 3 7% 1 40 93% Total Cases = 43PG5BPA>4067.71PG5BPA 4067.71Class = 1Class Cases Percent 0 3 100% 1 0 0.0% Total Cases = 3Class = 0Class Cases Percent 0 0 0.0% 1 40 100% Total Cases = 40Figure 2.7: CART model trained and built using synthetic measurements obtained in the fourth cycle after fault clearance.39Table 2.1: TSA indices selected in CART models built using data obtained in 5 cyclesafter fault clearance.Level CART 1 CART 2 CART 3 CART 4 CART 51 αBCH αBCH αBCH αBCH αBCH2σBCH ∆VBCH,BPA, ∆VBCH,BPA. ∆VBCH,BPA, JAES ,Vtotal,BPA. γBCH . ∆δBCH,BPA.3αBCH , αBCH , αBCH , αBCH , ∆δBCH ,Vtotal,BCH , βAES . βAES . βBPA, αBPA,βAES . γBPA. ∆VBCH,BPA.Table 2.2: Accuracy results of obtained CART models when tested using data fromvarious time instants.Test DataTrained CART ModelCART 1 CART 2 CART 3 CART 4 CART 5Cycle 1 99.50% 93.98% 80.15% 85.49% 74.87%Cycle 2 93.38% 99.18% 84.58% 84.58% 82.38%Cycle 3 96.78% 96.62% 99.61% 91.85% 90.39%Cycle 4 96.89% 96.71% 95.52% 98.97% 92.90%Cycle 5 96.41% 94.46% 94.41% 96.04% 99.43%2.8.2.1 Comparison with the Existing MethodsHere, we compare our CART results in Table 2.2 with existing studies in the literature. In[28] and [33], authors demonstrate accuracy rates of 97–98% using only pre-fault features.However, all probable faults are pre-defined and used as features to train and test DTmodels. Similarly, in [81], the fault location is used as a feature to obtain a predictionaccuracy rate of 97%. Unsurprisingly, the fault type and location are very effective featuresfor assessing transient stability. However, in practical implementation with a limitednumber of PMUs installed, it is impossible to determine the type and location of allcredible faults in real time. As shown in Table 2.2, models developed in this chapterlead to similar prediction accuracy rates without knowledge of the particular fault typeand location. In the case of using the pre-fault features in our study, the accuracy ratedecreases to 55.41%. The obtained DT is shown in Fig. 2.8. However, it can be observedthat the more important pre-fault features are power generation of the BCH system, active-power transfer between BCH and BPA, and active-power transfer from the northern BC tothe lower mainland. The obtained features and the thresholds can be used as a guide fortransient stability limits for the system operators. We note that all of these features canbe obtained from the SCADA system.In [25] and [34], authors use post-fault features obtained 1–2 s after fault-clearance.This time is reduced to 150–300 ms in [57]. The accuracy rates reported in these works arein the range of 93% with simple DTs and 99% with Ensemble DTs. However, the durationof simulation until normal end (in stable cases) or loss of synchronism (in unstable cases)is used as a feature, which is not available in real-time TSA. Models developed in thischapter take advantage of new PMU technologies that are capable of transmitting post-40Table 2.3: Computational times required for various CART models.Test DataComputation Time [s]CART 1 CART 2 CART 3 CART 4 CART 5Cycle 1 12 16 17 21 20Cycle 2 18 12 17 21 19Cycle 3 19 17 12 20 20Cycle 4 18 18 19 12 19Cycle 5 17 19 19 19 11fault measurements immediately after fault-clearance. As shown in Table 2.2, predictionaccuracy rates from simple DTs are comparable to prior work using ensemble DTs withoutknowledge of offline simulation duration.Pg,BCHClass Cases Percent 0 764 11.5% 1 5902 88.5% Total Cases = 6666γBCHClass Cases Percent 0 2 0.1% 1 2091 99.9% Total Cases = 2093Pg,BCH≤10860.49PNS,BCHClass Cases Percent 0 762 16.7% 1 3811 83.3% Total Cases = 4573αBCHClass Cases Percent 0 408 23.4% 1 1338 76.6% Total Cases = 1746PING-CUST Class Cases Percent 0 354 12.5% 1 2473 87.5% Total Cases = 2827PNS,BCH≤2741.02 PNS,BCH>2741.02Pg,BCH>10860.49βAESClass Cases Percent 0 249 14.9% 1 1417 85.1% Total Cases = 1666Pg, BCHClass Cases Percent 0 105 9.0% 1 1056 91.0% Total Cases = 1161PING-CUST≤2296 PING-CUST>2296Figure 2.8: CART model trained based on pre-fault synthesized features obtainedfrom SCADA system.2.8.2.2 Adaptive TrainingAs mentioned in Section 2.7, we train the CART model based on the synthetic measure-ments obtained from three different topologies. This can be viewed as having updatedthe model with two topology updates. In order to showcase the necessity of updating theCART model when new topologies arise, we train a CART model using only data from onetopology and subsequently test this model using data obtained from another topology. Inthis case, the prediction accuracy rate decreases to 94.62%. This result clearly shows thenecessity for an adaptive framework with the capability of updating the classification orregression model after topological changes, which is shown in Fig. 2.9. Such a frameworkis also important in the case of high penetration of renewable resources that leads to morefrequent changes in operating conditions. To account for these aspects, the trained modelscan be updated using the adaptive framework described in [33].41Preparing a data set based on base cases and predefined contingenciesTraining 5 cascading DTs using the pre- and post-fault features Initial developed CART modelsOffline TrainingIncorporate new training cases obtained after topology changes into data setTraining 5 cascading DTs Refined DT modelsUpdate of trained modelsTime-synchronized features obtained from PMUs TSA resultsFigure 2.9: Adaptive framework for online TSA.2.8.3 MARS ResultsSimilar to CART models, five different MARS models are obtained for the first five cyclesafter fault clearance. Also in this case, since there are synthetic measurements obtainedfrom three different topologies, we can assume that we have updated the obtained modelbased on the first topology for two more topologies. First, we select a simple model with amaximum of 15 possible BFs and no interactions between the indices, i.e., multiplicationof different BFs. As an example, the model obtained for the first cycle is shown below:42MARS 1 = −19.54− 1.925 ·max(0, 7.709− σBCH)− 97.74 ·max(0, αBCH − 0.319)+ 60.04 ·max(0, 0.319− αBCH)− 7.497 ·max(0, δBCH,BPAintercon,1 − 6.219)+ 0.665 ·max(0, 6.219− δBCH,BPAintercon,1 )− 10397 ·max(0, γBCH − 0.0138)+ 890.3 ·max(0, 0.0138− γBCH)+ 9.315 ·max(0, βBCH + 1.138)− 60.89 ·max(0,−1.138− βBCH)+ 50.73 ·max(0, αBPA + 0.050)+ 3453 ·max(0,−0.050− αBPA)− 22.41 ·max(0, σBPA − 8.368)− 0.057 ·max(0, PG2BPA − 1339.1). (2.40)Table 2.4: TSA indices selected in MARS models built using data obtained in 5 cyclesafter fault clearance.Model Transient Stability Assessment IndicesMARS 1 αBCH , βBCH , γBCH , σBCH , αBPA, δBCH,BPAintercon,1 , PG2BPA.MARS 2 αBCH , βBCH , σBCH , ωCOI,BCH , Vtotal,BCH , δBCH,BPAintercon,1 , PG3BPA.MARS 3 αBCH , σBCH , ωCOI,BCH , Vtotal,BCH , Utotal,BCH , δBCH,BPAintercon,1 , ψBCH,BPA.MARS 4 αBCH , βBCH , σBCH , Utotal,BCH , ∆VBCH,BPA, βAES , δBCH,BPAintercon,1 , δG3,PreBPA .MARS 5 αBCH , βBCH , γBCH , σBCH , δBCH,BPAintercon,1 , δBCH,BPAintercon,2 .The TSA indices selected in all five MARS models are shown in Table 2.4. The COI-referred and energy function-based indices have the most participation in the models,especially those related to BCH. Also, δBCH,BPAintercon,1 , the relative angle of the 500 kV busconnecting BCH to BPA, emerges in all five models. On the whole, we conclude that COI-referred and energy function-based indices from Sections 2.5.2 and 2.5.3 act as importantTSA characteristics in both classification and regression, more than direct measurementsfrom Section 2.5.1.The accuracy of the MARS model is evaluated by comparing the TSI values obtainedfrom MARS models on the benchmark computed via (2.29). The mean absolute errors(MAEs) obtained for each of the five models based on the test set data are reported inTable 2.5. Also, MAEs are reported for cases in which the test data from a cycle otherthan the cycle used to train the MARS model is fed to it. As is evident from Table 2.5, the431 2 3 4 5 6 71.522.533.544.55Maximum Number of InteractionsMean Absolute Error Maximum 15 BFsMaximum 30 BFsMaximum 45 BFsMaximum 60 BFs15 20 25 30 35 40 45 50 55 60 651.522.533.544.55Maximum Number of Basis FunctionsMean Absolute Error Maximum 1 InteractionsMaximum 2 InteractionsMaximum 3 InteractionsMaximum 7 InteractionsFigure 2.10: Sensitivity analysis of MARS model trained using synthetic measure-ments obtained in the second cycle after fault clearance to maximum numberof interactions and maximum number of BFs.Table 2.5: Accuracy results of obtained MARS models when tested using data fromvarious time instants.Test DataTrained MARS ModelMARS 1 MARS 2 MARS 3 MARS 4 MARS 5Cycle 1 5.323 10.72 9.167 10.04 9.786Cycle 2 5.326 4.998 6.944 7.714 7.645Cycle 3 8.102 7.678 5.333 6.479 7.034Cycle 4 11.16 9.657 6.747 5.026 6.672Cycle 5 15.40 12.30 7.030 7.544 4.834error is around 6 for all the models when the test data is from the correct cycle. Whenthe test data is from one cycle prior or after, mean absolute error (MAE) usually increasesto not more than 8. Similar to the CART models, the five sequential trained MARS modelsover the time frame of five cycles of the post-contingency period can be considered as anensemble of MARS models and the average output of all models is more reliable. Accordingto Table 2.6, the computational time required to build MARS models with 15 maximumBFs and 1 maximum interaction is approximately 8 s using the data set from the same44Table 2.6: Computational times required for various MARS models with 15 maxi-mum BFs and 1 maximum interaction.Test DataComputation Time [s]MARS 1 MARS 2 MARS 3 MARS 4 MARS 5Cycle 1 8 13 14 14 14Cycle 2 13 7 14 14 13Cycle 3 13 14 7 14 14Cycle 4 13 14 14 7 13Cycle 5 13 14 14 14 8Table 2.7: Accuracy results of the obtained MARS models with 60 maximum BFsand 6 maximum interactions.Model MARS 1 MARS 2 MARS 3 MARS 4 MARS 5MAE 2.479 1.784 1.725 2.004 1.982cycle (diagonal entries). However, the computational time increases to about 14 s in thecase of using training and testing data from different cycles (off-diagonal entries). Note,again, that SPM provides the total time required for the training and testing stages.In an effort to reduce MAEs in the MARS model, we have the option of increasing themaximum number of BFs and the maximum number of interactions. Increasing bothof these parameters in the models makes them more complex, but reduced MAE canbe achieved. We vary the maximum number of BFs in the range of 15 to 65 and themaximum number of interactions in the range of 1 to 7 and obtain MARS models for eachcase. In Fig. 2.10, we plot the sensitivity of model accuracy to the maximum numberof interactions and BFs for the MARS model trained using data obtained from the secondcycle after fault clearance. According to Fig. 2.10, as the maximum number of interactionsincreases to 3 or more and the maximum number of BFs increases to 50 or more, the MAEsdecrease to less than 2, significantly lower than the MAEs resulting from the simple modelsreported in Table 2.5. In Table 2.7, we report MAEs for MARS models obtained with 60BFs and 6 interactions. Also, referring to Table 2.8, we find that the computational timerequired to obtain the MARS models with 60 maximum BFs and 6 maximum interactionsis approximately 50 s.2.8.4 Sensitivity to NoiseIn order to assess the proposed methods’ sensitivity to noisy measurements, we add whiteGaussian noise to the data set with different levels of signal-to-noise ratio (SNR). Theobtained results are shown in Tables 2.9 and 2.10. As expected, as the SNR decreases, theaccuracy rates in the CART models decrease and the MAEs in the MARS models increase.However, the models are quite robust against injected measurement noise as compared torates reported in the existing literature. For example, in [57], accuracy rates decrease to81% when models are trained and tested in the presence of noise with 30 dB SNR.45Table 2.8: Computational times required for various MARS models with 60 maxi-mum BFs and 6 maximum interactions.Model MARS 1 MARS 2 MARS 3 MARS 4 MARS 5Computational53 53 49 50 50Time [s]Table 2.9: Accuracy results of obtained CART models when white Gaussian noise isadded to data set.SNR [dB] CART 1 CART 2 CART 3 CART 4 CART 5No noise 99.50% 99.18% 99.61% 98.92% 99.43%40 99.00% 98.93% 98.72% 98.58% 99.04%30 98.93% 98.11% 98.33% 98.65% 98.47%20 98.43% 98.08% 97.01% 98.33% 97.40%10 96.72% 96.79% 95.98% 97.93% 96.55%Table 2.10: Accuracy results of obtained MARS models when white Gaussian noiseis added to data set.SNR [dB] MARS 1 MARS 2 MARS 3 MARS 4 MARS 5No noise 5.323 4.998 5.333 5.026 4.83440 6.297 5.464 5.347 5.309 5.31530 6.015 5.501 5.566 5.328 5.49320 7.357 7.249 6.710 6.577 5.97510 8.768 8.405 7.878 7.569 7.4992.9 SummaryUsing PMU measurements obtained from the system, we develop practical CART and MARSmodels to assess TSA in real time. In order to train these models, we apply direct measure-ment, COI-referred, and energy function-based indices, which are chosen based on theiravailability in real time. Via case studies involving the full WECC system, we show thatthe trained CART and MARS models predict stable and unstable cases with high accuracyrates and provide a transient stability margin of the system based on the applied features,respectively. We observe that COI-referred and energy function-based features have moreimpact in the classification and regression models. Moreover, we study the sensitivity ofthe trained models to the data obtained from earlier or later sampling times. Avenuesfor future work include analyzing the robustness of the obtained models in the cases ofmissing or bad data, which are probable occurrences in a practical WAMS. Another com-pelling topic for future work is the use of trained models to assess the fidelity or qualityof real-time measurements.46Chapter 3Real-time Fault Detection andIdentification UsingSynchrophasors3.1 BackgroundThe transmission grid transfers large amounts of electric power from generators to loadcentres via a complex network of transmission lines. From time to time, these lines expe-rience faults caused by various effects, such as storms, lightning, and insulation aging andbreakdown [47]. In order to maintain system-wide power availability and quality, mainte-nance crews must swiftly repair damages to and further restore the faulted line. Vital toaccomplishing these tasks are real-time tools that can accurately pinpoint the fault loca-tion along the line of interest in a timely manner, which can be effectively implementedby a synchrophasor-based WAMS.Through modern WAMS, synchrophasors are available at the control centre of powersystems. Capability of a precise fault location at the control centre or at any other cen-tral location decreases restoration time significantly, but despite all of the comprehensivestudies and various proposed methods, the effect of synchrophasor accuracy on the preci-sion of fault location has not been investigated profoundly in the corresponding literature.To have a better grasp of this problem, we should consider fault occurrence and removalprocedures more precisely.By expanding PMUs throughout the power system and establishing new WAMSs incontrol centres, it has become more feasible to run various applications using synchropha-sors. For example, fast and accurate fault location is not only critical to power systemrestoration, but is also highly applicable to new transient stability algorithms.47The use of time synchronized voltage and current measurements can serve a number offunctions in electric power systems. Further, there is growing interest in using synchropha-sor systems that provide voltage and current phasors time-tagged to an absolute (common)time, to enable a plurality of real-time and non-real-time applications. In recent years,there has been a tremendous amount of research, development, and pilot demonstrationsof various applications using synchrophasor data. Further, it is clear that for productiongrade applications, reliable performance is critical. IEEE Std. C37.118-2005 [9] and its laterversions (e.g., IEEE Std C37.118.1-2011 [11] and IEEE Std C37.118.1a-2014 [10]) specifycertain quality flags and performance requirements for PMUs; however, these performancerequirements may not be sufficient for certain applications. For example, the accuracy ofa PMU during transient conditions on a faulted line is not well-defined. In fact, describinga waveform which is not a simple periodic sinusoidal is not possible using the basic phasordefinition. Accordingly, using synchrophasor measurements during a fault condition maylead to unreliable results as measurements may not be fair representations of the voltageand current waveforms on the line.One of the applications of interest for synchrophasor systems is fault location. Multi-terminal fault location using synchronized measurements at each end of the line has beenshown to be an effective way to calculate fault location [42, 43], helping speed up thetime to re-energization. When using synchrophasor data for this purpose, however, theaccuracy of the fault location is very much dependent on the accuracy of the synchrophasordata provided by the PMUs during the fault. As mentioned above, the synchrophasor datameasured by the PMU may or may not be useful (may misrepresent the waveform).In this chapter, we review fault location approaches, focus on impedance-based faultlocation methods applying sychrophasors, and elaborate on the use of the GoF metric asintroduced in [82] as a quantitative indication of the usefulness of PMU synchrophasordata for fault location application. We demonstrate how providing the GoF together withsynchrophasor data from a PMU can help improve the reliability and trustworthiness of asynchrophasor-based fault location application3.2 Approaches for FDIFollowing a fault occurrence on a transmission line, the line may be disconnected by pro-tective relays and circuit breakers. In this case, the maintenance crew should be sentto the fault location to investigate and repair any damages. Therefore, fast and accuratefault location is critical for transmission line restoration, power system reliability, and eco-nomic efficiency. Notably, more precise fault location has become feasible through moderncommunication systems and interchanging data at the ends of transmission lines [47].Generally, fault location approaches are divided into the following categories:48 Impedance-based approaches Travelling wave-based approaches High frequency approaches based on voltage and current waveforms generated by thefault AL approaches3.2.1 Impedance-based ApproachThe impedance-based approach, known as the simplest approach for fault location, usesphasors obtained from voltage and current signals at the terminals of a faulted line, plusthe line parameters. In this approach, the impedance of one segment of the faulted lineis calculated, allowing a calculation of the distance between one terminal and the faultlocation. This approach is widely considered, since it is practical and relatively inexpensiveto implement using already available instruments in substations [47].3.2.2 Travelling Wave-based ApproachThe travelling wave-based approach is also widely considered in literature. This approachis based on the correlation between the backward and forward voltage or current travellingwaves from the fault location to the faulted line terminals. A high frequency voltage orcurrent travelling wave is induced by the fault. In this approach, the propagation speedof the high frequency wave in the line is a critical element in fault location calculations.Additionally, the sampling rate in the measurement window should be high in this ap-proach, even though difficulties may arise in distinguishing the correct travelling wave.These issues make the travelling wave approach complicated and expensive, although theresults obtained through this approach have been proven accurate [47].3.2.3 High Frequency ApproachThe high frequency approach is in fact an extension on the travelling wave-based approach,based on high frequency components of voltage and current signals after fault incipience.When a fault occurs on a line, it generates voltage and current waveforms which covera wide range of frequencies. As the related frequency of the travelling wave componentsincreases, propagation speed increases and amplitude decreases. Therefore, different com-ponents reach the line terminal at various times. This method is based on a successiveidentification of the arrival of these different components at the measurement location [47].493.2.4 Automatic Learning ApproachThis last approach is based on AL methods. Similar to the discussions in 2.2.3, an extensivedata set should be generated here, covering all possible cases in terms of fault type, faultlocation, fault impedance, and fault inception time. Next, discriminating features shouldbe selected/extracted to make the pattern recognition process easier for the final stage oftraining a classification tool, such as an artificial neural network or DT [47].3.3 Impedance-based Fault Location throughSynchrophasor SystemExisting impedance-based methods can be broadly categorized into (i) single-end meth-ods and (ii) double- or multi-end methods. Single-end methods (see, e.g., [35, 36]) utilizemeasurements obtained from only one end of the faulted line and do not need data fromremote buses. These methods may experience inaccurate fault-location results. On theother hand, methods that utilize measurements from both ends (or multiple ends in caseof multi-terminal lines) of the faulted line are generally more accurate [37, 38]. Neverthe-less, unsynchronized measurements obtained from remote buses may cause fault-locationerrors, which is a shortcoming that has been identified in the respective literature [39, 40].An enabling technology that can address this challenge is the PMU, which provides syn-chrophasor measurements [41].Recognizing the potential for synchrophasors to reduce fault-location errors arisingfrom unsynchronized measurements, numerous approaches have been proposed to takeadvantage of their capabilities [42–44]. However, these studies do not address the impactof phasor measurement errors on the resulting fault-location accuracy or credibility. Atypical metric to quantify errors in the PMU output is the total vector error (TVE), whichmeasures the difference between the information from a PMU that describes a phasor andthe true phasor itself [9, 10]. Such a metric is relevant only if the actual voltage and currentsignal of interest is in a sinusoidal steady state. However, immediately following suddenevents, such as fault initiation and clearance, voltages and currents experience electricaltransients that render the sinusoidal steady-state assumption held by the TVE metric in-valid. To address this shortcoming, the GoF metric measures the difference between thetime-domain waveform reconstructed from phasor measurement and the actual voltage orcurrent signal [82]. By explicitly comparing time-domain signals, the GoF metric is appli-cable under both sinusoidal steady-state and electrically transient operating conditions.This is particularly relevant for real-time, reliability-critical applications—such as faultlocation—that use synchrophasors while the system experiences electrical transients.Given that double-end fault-location methods are easily implementable with moderncommunication infrastructure in today’s power systems, we focus on such a method similar50to that described in [42] and extend the method to locate (possibly unsymmetrical) faultsalong three-phase transmission lines using synchrophasors. Next, via a numerical exam-ple, we highlight that the accuracy of fault-location results depends on how well voltageor current waveforms can be approximated as phasors. In order to quantify the mismatchbetween actual time-domain waveforms and their phasor representations during electricaltransients, we use the GoF metric described in [82], as well as a variant. With GoF metricsin hand, we perform numerous simulations for various fault conditions and locations alonga transmission line to investigate the relationship between the accuracy of synchrophasor-based fault-location results and the GoF values for corresponding phasors. Through ex-tensive numerical case studies, we observe that synchrophasors with GoF values abovecertain thresholds provide sufficient confidence levels in the accuracy of fault-location re-sults. Thus, we recommend to include GoF metrics as additional prescribed quantities forPMUs to transmit to the fault-location application. Such a modification may also help toassess the accuracy of results from other synchrophasor-based applications that require ahigh level of confidence (as is often the case in production-grade real-time applications inthe electric power industry).3.4 PreliminariesFollowing the release of [9], synchrophasor algorithms have been further investigated withthe aim of improving the basic model in [9] to meet the static and dynamic performancerequirements assigned in industry standards (see, e.g., [83, 84]). Synchrophasors obtainedfrom PMUs can be categorized into two performance classes: measurement (M-class) andprotection (P-class) [10]. Since the focus of this study is on a specific application of PMUsfor system protection, we elect to use the basic P-class PMU model proposed in Annex Cof [9], from which measurements are extracted. However, we note that the choice of amore detailed PMU model would not affect our analysis or results significantly. In thisstudy, PMU measurements are used to estimate steady-state transmission-line impedancesand to locate unsymmetric faults. Through the fault-location application, we motivatethe need to incorporate GoF metrics for PMU measurements.3.4.1 PMU ModelingWe model a PMU in MATLAB Simulink, using a simple model as the contribution of thiswork is not focused on PMU modeling. The block diagram of the model for one phase isshown in Fig. 3.1. As shown in the figure, Fast Fourier Transform (FFT) is used in thefirst step to obtain the amplitude and phase at a frequency of 60 Hz1. Then, the phasevalue is further processed to obtain the final value of phase and frequency. The obtained1FFT is an efficient algorithm for computation of discrete Fourier transform.51FFTDerivativeMeanIntegrationInitial PhaseCT/VT Output+60AmplitudePhaseFrequency++-Figure 3.1: General diagram of a PMU for one phase.values are reported once in a cycle, i.e., every 16.666 ms. Since we are not using ROCOFin this study, we have not considered it in our model.While many PMUs may only measure positive sequence values, our PMU model measuresthe voltage and current variables of each phase separately. The voltage phasor of phasex at bus l and time k is shown as V xj [k] where x ∈ {a, b, c}. Also, the current phasor ofphase x out of buses j and l through line (j, l) at time k is shown by Ix(j,l)[k] and Ix(l,j)[k],respectively. Based on these phasors, we are able to regenerate the phase voltage andcurrent signals in each cycle and compare them with their original measured signal in thenext steps of our study. Finally, the model is implemented in MATLAB Simulink. TheSimulink block diagram for one of the three phases is illustrated in Fig. 3.2.3.4.2 Measurement-based Estimation of Line ImpedanceSteady-state transmission-line impedance values are prone to variations of up to 40% un-der various ambient and loading conditions [85]. In order to ensure accurate fault-locationresults, the pre-fault steady-state line impedance for each phase must be updated peri-odically to match evolving operating conditions. Inspired by the parameter-free methodproposed in [42], in which measurements in the pre-fault period are used instead of relyingon a priori knowledge of line parameters, below we estimate steady-state lumped line pa-rameters assuming three-phase balanced impedances. We assume that P-class PMUs areinstalled at buses m and n, which are connected via the transmission line of interest, andthat per-phase bus-voltage and current-flow measurements are available.Suppose PMU phasor measurements are reported at intervals of ∆t, i.e., at timestk = k∆t, k = 1, 2, . . . Denote measured voltages (as phasors) at bus m and time tkfor phases a, b, and c as V am[k], Vbm[k], and Vcm[k], respectively. Further denote themeasured currents (as phasors) through line (m,n) (assume positive flow from bus mto n measured at bus m) at time tk in phases a, b, and c as Ia(m,n)[k], Ib(m,n)[k], and521Vabc1Mag_a-K-Gain32Pha_a-K-Gain12FreqInMean(Variable Frequency) -K-Gain13FreqInMean(Variable Frequency)93freq_rate_a4freq_a60 Constant31s Integrator1-K- Gain17-K-Gain18TransportDelay TransportDelay1ScopeScope15Mag_b6Pha_b7freq_rate_b8freq_bScope2Scope39Mag_c10Pha_c11freq_rate_c12freq_c|u|∠uFourier-K-Gain1-K-Gain2-K-Gain10-K-Gain5FreqInMean(Variable Frequency)1 -K-Gain6FreqInMean(Variable Frequency)260 Constant11s Integrator2-K- Gain7-K-Gain8TransportDelay2 TransportDelay3Scope4Scope5Scope6Scope7|u|∠uFourier1-K-Gain4-K-Gain9-K-Gain11-K-Gain15FreqInMean(Variable Frequency)3 -K-Gain16FreqInMean(Variable Frequency)460 Constant21s Integrator3-K- Gain19-K-Gain20TransportDelay4 TransportDelay5Scope8Scope9Scope10Scope11|u|∠uFourier2-K-Gain14-K-Gain21Zero-OrderHold2Zero-OrderHold1Zero-OrderHold3Zero-OrderHold4Zero-OrderHold5Zero-OrderHold6Zero-OrderHold7Zero-OrderHold8Zero-OrderHold9Zero-OrderHold10Zero-OrderHold11Zero-OrderHold12Figure 3.2: PMU model simulated in MATLAB Simulink.Ic(m,n)[k], respectively. Collect per-phase voltage measurements at bus m into Vm[k] =[V am[k], Vbm[k], Vcm[k]]T; also collect per-phase measurements of currents in line (m,n) intoI(m,n)[k] = [Ia(m,n)[k], Ib(m,n)[k], Ic(m,n)[k]]T. In symmetrical components, the voltage at busm is defined as V sm[k] := T−1Vm[k], where T is the so-called symmetrical components trans-formation matrix [86]. Similarly, for currents through line (m,n), Is(m,n)[k] = T−1I(m,n)[k].With these definitions in place, consider a three-phase line with balanced impedances,for which the circuit in symmetrical compo ents is decoupled into three independentimpedance networks, i.e., [86]V sm[k]− V sn [k] = diag(Zs(m,n))Is(m,n)[k], (3.1)where the vector Zs(m,n) contains the positive-, negative-, and zero-sequence networkimpedances, and diag(Zs(m,n)) denotes the diagonal matrix formed with diagonal entriescomposed of entries of Zs(m,n)[k]. Then, using per-phase voltage and current measurementsobtained at time tk, we can compute1Zs(m,n)[k] = diag(Is(m,n)[k])−1(V sm[k]− V sn [k]) (3.2)= diag(T−1I(m,n)[k])−1T−1(Vm[k]− Vn[k]),where the second equality above results by immediate application of the symmetricalcomponents transformation. From (3.2), we can straightforwardly obtain an up-to-dateestimate of the abc-frame three-phase impedance matrix asZ(m,n)[k] = Tdiag(Zs(m,n)[k])T−1. (3.3)1We use the fact that diag(x)y = diag(y)x, where x and y are vectors of appropriate dimension.53Z(m,n)anmI ),(bnmI ),(cnmI ),(amVbmVcmVanVbnVcnVamnI ),(bmnI ),(cmnI ),((a) Z(m,n)diag(d) Z(m,n)diag(lI3-d)amVbmVcmVanmI ),(bnmI ),(cnmI ),(anVbnVcnVamnI ),(bmnI ),(cmnI ),(afnmZ ),,(bfnmZ ),,(cfnmZ ),,((b)Figure 3.3: Circuit model of three-phase transmission line (m,n) (see, e.g., [5]) (a)before fault, and (b) during fault. The fault location (characterized as thedistance d from bus m) and the fault-to-ground impedances (Za(m,n),f , Zb(m,n),f ,and Zc(m,n),f ) are unknown quantities. Once solved, they provide the faultlocation and identify the affected phases.The ideas presented above are summarized in Fig. 3.3a. Since we consider the caseof balanced impedances, the diagonal entries of Z(m,n)[k] (representing self impedances)have the same value, and the off-diagonal entries (representing mutual impedances) arealso equal to each other.Example 1. In this example, we consider the canonical two-area test power system (see,e.g., [5]). Particularly, we are interested in one of the two identical 230-kV three-phasebalanced transmission lines of length 220 km connecting the two area via buses m and n. Atransposed distributed line model is used to simulate the transmission line. The PMUs aremodelled in the MATLAB Simulink environment using the FFT function at the nominal60-Hz frequency to obtain phasor magnitudes and phase angles. Then, these values arereported once per electrical cycle, i.e., ∆t = 16.667 ms. Four PMU models are connected toboth ends of the test transmission line using CT and VT models, part of which is shown inFig. 3.4. Then, via (3.2), we use the voltage- and current-phasor measurements obtainedat steady state to compute the positive-, negative-, and zero-sequence impedance valuesas Zs(m,n) = [116.88∠84.243◦, 116.88∠84.243◦, 617.59∠55.004◦]T Ω. Indeed, the resultingsymmetrical impedance values precisely match those obtained based on the predefinedparameter values for the line model. The abc-frame three-phase impedance matrix can beeasily obtained using (3.3) as Z(m,n), in which the diagonal elements are 276.51∠62.931◦Ω54413 MW---->Discrete,Ts = 1e-06 s.3 B1230e3V0.9916pu15.06deg.3 B2230e3V1.003pu-11.43deg.+Line 2(220 km)+Line 1b(110 km)+Line 1a(110 km)A B CA B CFaultABCabcBrk1ABCabcB2ABCabcB1ABCArea 2 A B CArea 1ABCabc Brk21 2+ +Saturable Transformer1 2+ +Saturable Transformer11 2+ +Saturable Transformer2v+-Voltage Measurement4v+-Voltage Measurement5v+-Voltage Measurement612++Saturable Transformer312++Saturable Transformer412++Saturable Transformer5v+ -Voltage Measurement7v+ -Voltage Measurement8v+ -Voltage Measurement9Scope5Switch+1ohmSwitch1+1ohm1Switch2+1ohm2Scope20Scope19Scope17Scope18Scope14Scope15Scope26Scope27VabcMag_aPha_afreq_rate_afreq_aMag_bPha_bfreq_rate_bfreq_bMag_cPha_cfreq_rate_cfreq_cSubsystemScope22Scope23Scope24Scope25Scope30Scope31Scope32Scope33Scope35Scope34Scope8Scope16Scope39Scope38Scope36Scope37VabcMag_aPha_afreq_rate_afreq_aMag_bPha_bfreq_rate_bfreq_bMag_cPha_cfreq_rate_cfreq_cSubsystem1v+-Voltage Measurement3A2v+-Voltage MeasurementA2v+-Voltage Measurement1A2abc|u|∠uSequence Analyzer2 Scope40i+ -Current Measurement2A2i+ -Current Measurement3A2abc|u| ∠uSequence Analyzer3Scope41Scope42|u|∠uFourier6 Scope45Scope43|u|∠uFourier7 Scope44Scope46|u|∠uFourier8 Scope47v+-Voltage Measurement11A2 Scope48|u|∠uFourier9 Scope52v+-Voltage Measurement22A2 Scope53|u|∠uFourier10 Scope49v+-Voltage Measurement33A2 Scope51|u|∠uFourier11 Scope501 2+ +Saturable Transformer61 2+ +Saturable Transformer71 2+ +Saturable Transformer8v+-Voltage Measurement1v+-Voltage Measurement2v+-Voltage Measurement3Scope54Scope55Scope56Scope61Scope62VabcMag_aPha_afreq_rate_afreq_aMag_bPha_bfreq_rate_bfreq_bMag_cPha_cfreq_rate_cfreq_cSubsystem2Scope57Scope58Scope59Scope60Scope63Scope64Scope65Scope6612++Saturable Transformer912++Saturable Transformer1012++Saturable Transformer11v+ -Voltage Measurement10v+ -Voltage Measurement11v+ -Voltage Measurement12Scope67Switch3+1ohm3Switch4+1ohm4Switch5+1ohm5Scope72Scope71Scope69Scope70Scope74Scope73Scope79Scope68Scope78Scope77Scope75Scope76VabcMag_aPha_afreq_rate_afreq_aMag_bPha_bfreq_rate_bfreq_bMag_cPha_cfreq_rate_cfreq_cSubsystem3i+ -Current Measurement1A2Figure 3.4: Part of two-area power system with VT, CT and PMU models simulatedin MATLAB Simulink.and the off-diagonal elements are 172.91∠48.669◦Ω. 3.4.3 Fault Location and Faulted-phase IdentificationConsider a three-phase transmission line (m,n) with total length ` and steady-stateimpedance matrix estimated as Z(m,n) from (3.3), with an unknown fault imposed, asshown in Fig. 3.3b. In this section, we simultaneously locate the fault and estimate thefault-to-ground impedance using synchrophasors obtained at buses m and n. The esti-mated fault location and impedance values also indicate the phases that are affected bythe fault. To this end, denote the unknown distance from bus m to the fault location attime tk as d[k] = [da[k], db[k], dc[k]]T, where da[k], db[k], and dc[k] represent the distancesin phases a, b, and c, respectively. Also, denote the unknown fault-to-ground impedance attime tk as Z(m,n),f [k] = [Za(m,n),f [k], Zb(m,n),f [k], Zc(m,n),f [k]]T, where Za(m,n),f [k], Zb(m,n),f [k],and Zc(m,n),f [k] represent fault-to-ground impedances in phases a, b, and c, respectively.The variables d[k] and Z(m,n),f [k] are unknowns to be solved in this problem. Based on55the circuit in Fig. 3.3b, we can express three-phase KVL equations from bus m to n attime tk asVm[k]− Vn[k] = Z(m,n) (diag(d[k])) I(m,n)[k]− Z(m,n)diag(`13 − d[k])I(n,m)[k]= Z(m,n)diag(I(m,n)[k] + I(n,m)[k])d[k]− `Z(m,n)I(n,m)[k]=: A1[k]d[k]− c[k], (3.4)where 13 = [1, 1, 1]T. Similarly, we can write another KVL equation at time tk from busm to ground asVm[k] = Z(m,n)diag(d[k])I(m,n)[k]+ diag(Z(m,n),f [k])(I(m,n)[k] + I(n,m)[k])= Z(m,n)diag(I(m,n)[k])d[k]+ diag(I(m,n)[k] + I(n,m)[k])Z(m,n),f [k]=: A2[k]d[k] +A3[k]Z(m,n),f [k]. (3.5)Next, combine (3.4) and (3.5) in matrix form as[Vm[k]− Vn[k] + c[k]Vm[k]]=[A1[k] 0A2[k] A3[k]][d[k]Z(m,n),f [k]]. (3.6)At each time tk, the system of equations in (3.6) can be solved simultaneously to obtainestimates of the fault location and the fault-to-ground impedances. As we will show inExample 2, we can locate the unknown and unsymmetric fault via the solution of (3.6).Note that the solvability of (3.6) can be ensured by checking that A1[k] and A3[k] areinvertible. The matrix A1[k] is invertible as long as the entries of Zs(m,n) are nonzero, andthe matrix A3[k] is invertible as long as the current in each phase of the transmission lineis nonzero.Example 2. In this example, we again consider the transmission line connecting buses mand n from Example 1. To illustrate the proposed fault location method summarizedas (3.6), we choose simulation cases covering five types of faults: three phase to ground(abcg), double phase to ground (abg), single phase to ground (ag), three phase (abc), anddouble phase (ab). In each case, unbeknownst to operators, the fault occurs at a distanceof 60 km from bus m with fault-to-ground resistance 5 Ω at time t = 0 s, and the faultsustains for three cycles before circuit breakers are tripped. Phasor measurements are56collected at times tk = k∆t, k = 1, 2, 3, i.e., at the end of each electrical cycle during thefault.We first report and discuss results using measurements from the third cycle, i.e., inmeasurement window 2∆t < t ≤ 3∆t. The unknown variables d[3] and Z(m,n),f [3] aresolved via (3.6) and results are reported in Table 3.1, based on which we make the followingobservations: Fault location: For all five fault types, the estimated distances are nearly identical tothe actual location with errors of less than 0.1% normalized to the line length of 220 km. Faulted-phase identification: For phases that remain intact, the estimated distances frombus m are unreasonably large, i.e., they are significantly greater than the line length.The corresponding fault-to-ground impedances have extremely large magnitudes, i.e.,they represent open circuits. These observations are indicative of fault signature andhelp to identify the faulted phases.On the other hand, Table 3.2 reports fault location and impedance estimates obtainedwith synchrophasor measurements from the first cycle, i.e., 0 < t ≤ ∆t, of the during-faultperiod. By comparing Tables 3.1 and 3.2, we find that the results from the former aremuch more accurate than those in the latter. Of particular concern is that errors in theestimated fault location grows from less than 0.1% to more than 2.0%, a sizeable increase.3.4.4 Motivation for a GoF MetricDuring faults or other fast transients, phasors extracted from PMUs may not accuratelyrepresent current/voltage waveforms as shown in Fig. 3.5a, where the current waveformreconstructed from phasor measurements is depicted by the dash-dot trace, while theoriginal current signal is represented by the solid trace. We observe that the currentwaveform reconstructed from phasor measurements reasonably matches the actual signalin the third cycle during the fault, but not the first. Such discrepancies can lead tosignificant errors in applications that use these measurements, as we reveal in Example 2with respect to fault location. Hence, the question arises as to how the fault-locationapplication can assess the credibility of its results. As we show next, this problem canbe tackled by appropriately evaluating the “goodness” of the fit between the signal thatis reconstructed from the phasor measurements and the actual signal based on which thephasor was obtained. Such metrics quantify the accuracy of sychrophasors, which can inturn be related to the accuracy of fault-location results. In this way, these metrics can helpindicate whether or not fault-location results obtained using particular synchrophasors arecredible [80].57Table 3.1: Fault location and fault-impedance estimation using measurements from the third cycle of during-fault period.Fault Typeda[3] db[3] dc[3] Za(m,n),f[3] Zb(m,n),f[3] Zc(m,n),f [3]km % Error km % Error km % Errorabcg 60.060 0.027 60.080 0.036 60.007 0.003 4.332∠3.607◦ 5.183∠6.687◦ 4.945∠0.470◦abg 59.978 0.010 59.995 0.002 2483.7 − 5.683∠− 10.112◦ 5.020∠33.005◦ 1.803× 107∠− 96.918◦ag 59.985 0.007 4044.4 − 1186.4 − 4.921∠0.943◦ 5.370× 106∠54.002◦ 1.010× 107∠− 18.311◦abc 60.059 0.027 60.080 0.036 60.007 0.003 4.332∠3.607◦ 5.183∠6.688◦ 4.945∠0.470◦ab 59.824 0.080 59.864 0.062 2414.4 − 22.088∠− 11.559◦ 19.038∠160.902◦ 7.330× 106∠− 159.306◦Table 3.2: Fault location and fault-impedance estimation using measurements from the first cycle of during-fault period.Fault Typeda[1] db[1] dc[1] Za(m,n),f[1] Zb(m,n),f[1] Zc(m,n),f [1]km % Error km % Error km % Errorabcg 64.999 2.223 61.023 0.465 57.978 0.919 30.766∠5.365◦ 20.541∠34.746◦ 46.341∠38.725◦abg 63.619 0.736 63.290 1.463 13425 − 30.810∠− 6.595◦ 21.816∠45.561◦ 4.869× 107∠− 145.285◦ag 62.997 1.362 15275 − 14420 − 61.283∠10.380◦ 1.354× 108∠58.313◦ 5.729× 107∠− 60.366◦abc 64.998 2.272 61.023 0.465 57.978 0.919 30.766∠5.365◦ 20.541∠34.746◦ 46.341∠38.725◦ab 63.875 1.761 63.932 1.787 2673.3 − 38.005∠− 16.708◦ 17.498∠88.404◦ 7.616× 106∠− 166.954◦583.5 GoF Metric for Synchrophasor MeasurementsAlong a measurement window of length ∆t, the PMU conducts FFTs on time-domain per-phase voltage and current waveforms and produces corresponding phasor representationswith amplitude and phase components [9]. However, as highlighted in Section 3.4.4, suchphasor representations may not accurately reflect the time-domain waveforms if there arefast transients in the voltages and currents of interest. In such a case, the real-time appli-cation would benefit from additional quantities that describe the quality (i.e,. accuracy)of the synchrophasors. To this end, the so-called goodness-of-fit metric evaluates the mis-match between the actual voltage or current waveform and the corresponding waveformreconstructed from their respective phasor measurements obtained from PMUs [82]. In thissection, we define this metric (and a variant), highlight their key properties, demonstratetheir usage via several numerical examples, and finally propose them as additional outputsto be provided by PMUs.3.5.1 Definition of the GoF MetricDenote the actual a-, b-, and c-phase voltage waveforms at bus m at time t by v̂am(t), v̂bm(t),and v̂cm(t), respectively. Similarly, let îa(m,n)(t), îb(m,n)(t), and îc(m,n)(t) denote the actualcurrents in line (m,n). Within the measurement window (k − 1)∆t < t ≤ k∆t or simplyk, a PMU at bus m samples the actual per-phase voltage and current waveforms at a rateof N = 166 samples-per-measurement-window (i.e., sampling time of ∆τ = 100µs) [83].Then, the PMU conducts FFTs on the collected samples to produce corresponding pha-sor representations at time tk = k∆t, which are sent to the application [87]. Denotethe resultant a-, b-, and c-phase voltage phasors as V am[k] = |V am[k]|∠φam[k], V bm[k] =|V bm[k]|∠φbm[k], and V cm[k] = |V cm[k]|∠φcm[k], respectively, with steady-state frequenciesωam[k], ωbm[k], and ωcm[k], respectively. Further define |Vm[k]| = [|V am[k]|, |V bm[k]|, |V cm[k]|]T,φm[k] = [φam[k], φbm[k], φcm[k]]T and ωm[k] = [ωam[k], ωbm[k], ωcm[k]]T. Notation for measuredphasor representations of per-phase current flows are established similarly.With the above notation in place, we develop GoF-related concepts with respect to a-phase quantities, but note that b- and c-phase quantities would be evaluated analogously.Using the phasor magnitude, frequency, and phase information from a PMU, we can recovera corresponding reconstructed time-domain signal. For example, the reconstructed a-phasevoltage and current waveforms are expressed asvam(t) = |V am[k]| · cos(ωam[k]t+ φam[k]), (3.7)ia(m,n)(t) = |Ia(m,n)[k]| · cos(ωam[k]t+ θa(m,n)[k]), (3.8)where (k − 1)∆t < t ≤ k∆t. The reconstructed signals in (3.7) and (3.8) differ from their59-0.08 -0.06 -0.04 -0.02 0 0.02 0.04 0.06-500005000050100150200(a)-0.08 -0.06 -0.04 -0.02 0 0.02 0.04 0.06-500005000050100150200(b)Figure 3.5: (Example 3 and 4) Phase-a actual current waveform and that recon-structed from PMU phasor measurements (a) without and (b) with DC-offsetmodification. Corresponding GoF values are superimposed. Fault initiates atthe beginning of measurement window k = 1.corresponding actual waveforms by∆vam(t) = v̂am(t)− vam(t), (3.9)∆ia(m,n)(t) = îa(m,n)(t)− ia(m,n)(t). (3.10)Such discrepancies can be represented by the GoF as a single quantity over the measurementwindow k. As an example, the GoF metric (in [dB]) for the a-phase voltage and current60are computed as [82]ϕam[k] = 20 log|V am[k]|√1N−MN∑j=1(∆vam(tk−1 + j∆τ))2 , (3.11)ψa(m,n)[k] = 20 log|Ia(m,n)[k]|√1N−MN∑j=1(∆ia(m,n)(tk−1 + j∆τ))2 , (3.12)where N represents the number of samples in one measurement window (166 samples here)and M represents the number of estimated parameters to reconstruct the waveform. Inour setting, M = 3, where the estimated parameters are amplitude, phase, and frequency.Finally, ∆τ denotes the PMU internal sampling time in s. The GoF metric, as definedin (3.11) and (3.12), evaluates to a large value when the actual and reconstructed signalsare well matched. The logarithmic function amplifies the range of GoF values in the casewith greater mismatch between actual and reconstructed signals.Next, via a numerical example, we illustrate how the GoF metric quantifies the mis-match between the PMU phasor measurement and the corresponding actual time-domainsignal.Example 3. Consider the transmission line of interest in Fig. 3.3b from Example 2. Att = 0 s, we initiate a three-phase-to-ground fault with the fault resistance of 5 Ω on theline (m,n) (that is, da = db = dc = 60 km), which sustains for three cycles. Based onsimulated actual waveforms and those reconstructed from PMU measurements, we evaluateGoFs in pre-fault, during-fault and post-fault operating conditions. Figure 3.5a containsthe actual a-phase current waveform îa(m,n)(t) (solid trace) and the corresponding oneia(m,n)(t) reconstructed from phasor measurements (dash-dot trace). The GoF quantitiesare evaluated in each measurement window, which is marked by circles. Based on a visualexamination of Fig. 3.5a, we make the following observations: Pre-fault period (k ≤ 0): The waveform reconstructed using PMU measurements matchesthe actual waveform obtained from simulations. In this period where the mismatch isminimal, evaluated GoFs ψa(m,n)[k] ≈ 50 dB, k ≤ 0. During-fault period (1 ≤ k ≤ 3): Immediately following the fault initiation time t = 0 s,the mismatch between the actual and reconstructed waveforms grows, so that the GoFevaluated in the first during-fault cycle is ψa(m,n)[1] = 7.5 dB. As the transients settle inmeasurement windows k = 2, 3, the match between the two waveforms improves, and61the GoFs correspondingly increase to ψa(m,n)[k] = 19 dB for k = 2, and ψa(m,n)[k] = 31 dBfor k = 3. Post-fault period (k > 3): After 3 cycles, i.e., at t = 50 ms, the fault is cleared byopening circuit breakers at both ends of the faulted line at near zero current. Followingthe disconnection of one line, the current through that line decreases sharply. Again,such a sudden change causes the mismatch between the waveforms to increase, and inturn leads to lower GoF in measurement window k = 4.Due to the nature of the system and the faults under consideration, voltages turn out tobe well regulated throughout the fault period as shown in Fig. 3.6a. Consequently, theactual and reconstructed waveforms may match quite well. As such, we do not see roomdiscuss about voltage waveforms in detail like what we did for current waveforms. 3.5.2 Accounting for DC OffsetIn Example 3, we observe that while there is nearly perfect match between the recon-structed and the actual waveforms in sinusoidal steady state, the mismatch is quite notablein the during-fault period, when the system experiences electrical transients. Recall thatreconstructed signals (see, e.g., dash-dot trace in Fig. 3.5a) are recovered from the ampli-tude, phase, and frequency values, which are obtained via the PMU, using (3.7) and (3.8).However, the three aforementioned parameters do not encode information about the di-rect current (DC) offset of the actual waveform, which may be nonzero during electricaltransients. As a result, the GoF metric may be decrease in value simply due to a DC offset,and not errors in the phasor information reported. For example, this intuition is evidentvia visual inspection of Fig. 3.5a for k = 1. Thus, we devote some time to investigate theeffects of such DC offsets on the GoF metric.The DC offset can be easily computed by the PMU at time tk for the previous measure-ment window (k − 1)∆t < t ≤ k∆t. As an example, DC offsets for a-phase voltage andcurrent are evaluated as:∆vam[k] =1NN∑j=1∆vam(tk−1 + j∆τ), (3.13)∆ia(m,n)[k] =1NN∑j=1∆ia(m,n)(tk−1 + j∆τ), (3.14)respectively. These can be interpreted as the mean value of the difference between theactual and reconstructed signals in measurement window k. Now, using the DC-offsetinformation in conjunction with PMU phasor measurements (i.e., amplitude, phase, andfrequency), we hope that the reconstructed waveforms would more closely match the actual62ones. For example, for (k−1)∆t < t ≤ k∆t, the reconstructed a-phase voltage and currentwaveforms can be expressed asvam(t) = vam(t) + ∆vam[k] (3.15)ia(m,n)(t) = ia(m,n)(t) + ∆ia(m,n)[k], (3.16)respectively, where vam(t) and ia(m,n)(t) are obtained from (3.7) and (3.8), respectively.Accordingly, the mismatch between the actual waveforms and the ones reconstructed fromsynchrophasors (accounting for DC offset) are computed as v̂am(t)− vam(t) and îa(m,n)(t)−ia(m,n)(t). Finally, these can be substituted into (3.11) and (3.12) to obtain updated GoFmetrics ϕam[k] and ψa(m,n)[k] for the a-phase voltage and current, respectively.3.5.3 Effect of DC OffsetExample 4. Consider the same system and fault scenario as in Example 3, but in ad-dition to |Ia(m,n)[k]|, θa(m,n)[k], and ωam[k], assume that the DC offset obtained from (3.14)is available. The reconstructed current waveform is superimposed onto the actual one inFig. 3.5b. In this case, in contrast to Fig. 3.5a, the reconstructed and actual waveformsmatch much more closely, with the greatest improvement in the during-fault period, espe-cially in measurement window k = 1. Correspondingly, as compared with Example 3, GoFvalues are greater in the during-fault period, as depicted in Fig. 3.5b. Specifically, in theproblematic measurement window k = 1 from Example 3, the GoF increases from 7.5 dBto 20.5 dB.The effect of DC-offset remove on voltage waveforms is shown in Fig. 3.6b. As it isobvious in the figure, unlike the current waveforms, the effect of DC-offset on the mismatchis not considerable for voltage waveforms and has not been considered in the next stagesof the study. Previously, by making available the DC offset of the actual signal obtained from (3.13)and (3.14), we significantly improved the match between the actual waveform and thatreconstructed from PMU measurements, as illustrated in Example 4. On the other hand,as we show next via Examples 5 and 6, respectively, the problem of matching the recon-structed waveform to the actual one is further complicated if (i) the fault initiation timedoes not coincide with the measurement window boundary, or (ii) the current transformer(CT) saturates due to large line currents.3.5.4 Effect of Fault Initiation TimeExample 5. In this example, we consider the same system and fault scenario as in Ex-amples 3 and 4, except the fault initiation time is set to t = 6.94 ms, which is close to the63-0.08 -0.06 -0.04 -0.02 0 0.02 0.04 0.06-2-10123105050100150200250300(a)-0.08 -0.06 -0.04 -0.02 0 0.02 0.04 0.06-2-10123105050100150200250300(b)Figure 3.6: (Example 3 and 4) Phase-a actual voltage waveform and that recon-structed from PMU phasor measurements (a) without and (b) with DC-offsetmodification. Corresponding GoF values are superimposed. Fault initiates atthe beginning of measurement window k = 1.centre of the measurement window k = 1, as shown in Fig. 3.7. Evident from a visualinspection of Fig. 3.7a, in which the reconstructed waveform is obtained using only cur-rent phasor estimates, the substantial mismatch in measurement window k = 1 cannot becompletely obviated by a simple DC offset. This observation is substantiated in Fig. 3.7,where the reconstructed waveform incorporates the DC offset in the actual signal. Specif-ically, in measurement window k = 1, the GoF with the DC-offset modification is 9.7 dB64-0.08 -0.06 -0.04 -0.02 0 0.02 0.04 0.06-500005000050100150200(a)-0.08 -0.06 -0.04 -0.02 0 0.02 0.04 0.06-500005000050100150200(b)Figure 3.7: (Example 5) Phase-a actual current waveform and that reconstructedfrom PMU phasor measurements (a) without and (b) with DC-offset mod-ification. Corresponding GoF values are superimposed. Fault initiates att = 6.94 ms, within measurement window k = 1.as compared to 4.4 dB in Fig. 3.7a. Indeed, while the improvement in the match is lesssignificant than that observed in Example 4, by considering the DC offset, the mismatchis still reduced.Similar to the effect of DC offset on voltage waveforms, the effect of fault initiationtime on voltage waveforms is negligible as illustrated in Fig. 3.8. 65-0.08 -0.06 -0.04 -0.02 0 0.02 0.04 0.06-2-10123105050100150200250300(a)-0.08 -0.06 -0.04 -0.02 0 0.02 0.04 0.06-2-10123105050100150200250300(b)Figure 3.8: (Example 5) Phase-a actual voltage waveform and that reconstructedfrom PMU phasor measurements (a) without and (b) with DC-offset mod-ification. Corresponding GoF values are superimposed. Fault initiates att = 6.94 ms, within measurement window k = 1.3.5.5 Effect of CT SaturationExample 6. The measurement accuracy of monitoring and protection devices such asPMUs and relays depends on the quality of analog signals obtained from CTs [88]. Exam-ples 3–5 assumed an ideal CT model without saturation. In reality, due to the nonlinearcharacteristic of its core, the CT is vulnerable to the saturation phenomenon when line-current magnitudes are much higher than their steady-state values [89]. In this example,66-0.08 -0.06 -0.04 -0.02 0 0.02 0.04 0.06-101104-505(a)-0.08 -0.06 -0.04 -0.02 0 0.02 0.04 0.06-500005000050100150200(b)-0.08 -0.06 -0.04 -0.02 0 0.02 0.04 0.06-500005000050100150200(c)Figure 3.9: (Example 6) CT flux saturation limit is set to 2µV·s. (a) Phase-a CTmagnetizing current and core flux. (b) Phase-a actual current waveform andthat reconstructed from PMU phasor measurements without (c) with DC-offset modification. Fault initiates at the beginning of measurement windowk = 1.67-0.08 -0.06 -0.04 -0.02 0 0.02 0.04 0.06-101104-505(a)-0.08 -0.06 -0.04 -0.02 0 0.02 0.04 0.06-500005000050100150200(b)-0.08 -0.06 -0.04 -0.02 0 0.02 0.04 0.06-500005000050100150200(c)Figure 3.10: (Example 6) CT flux saturation limit is set to 2µV·s. (a) Phase-aCT magnetizing current and core flux. (b) Phase-a actual current waveformand that reconstructed from PMU phasor measurements without (c) withDC-offset modification. Fault initiates at t = 6.94 ms, within measurementwindow k = 1.68we consider the same scenario as in Examples 4 and 5, except a non-ideal CT model isused and the core flux saturation limit is set to 2µV·s. As shown in Fig. 3.9a, follow-ing the increase of the core flux to an amplitude more than 2µV·s, the CT saturatesand the magnetizing current increases sharply. Then, as shown in Fig. 3.9c, the satu-ration phenomenon leads to non-sinusoidal current waveforms (as seen by the PMU) inthe during-fault period. Again, we note that the mismatch between the original and thereconstructed waveforms cannot be removed by a simple DC offset. We have shown theresults of saturation in for fault initiation time t = 6.94 ms in Fig. 3.10. The results donot show in siginificant difference with Fig. 3.9. Via Examples 4–6, we observe that the GoF metric quantifies the error of the reportedsynchrophasor, and GoF values may be improved when the reconstructed signal considersthe mean DC offset in a measurement window (k − 1)∆t < t ≤ k∆t. Thus, we hypoth-esize that GoF metrics can serve as a predictor for the performance quality of real-timeapplications, such as fault location.Finally, we note that some of the distortions in typical substation waveform readingsfrom, e.g., digital fault recorders, are due to measurement noise. In such a case, GoF alsohelps to quantify the match between the original and regenerated current waveforms, butfor the particular application in this paper (fault location), signal levels (fault currents)are generally high and result in fairly good signal-to-noise ratio. Accordingly, we focusour attnetion on waveform distortion and DC offset instead of noise-related issues. Forfuture work in, e.g., state estimation, we intend to fully investigate the impact of noise onGoF metrics and application accuracy.As a conclusion, trying to use synchrophasor data in practical real-world applications,we have observed that the GoF can be low for various reasons including: DC offsets (decaying DC) during the early cycles of faults with long time constants, Amplitude jumps in the first cycle of fault, Distortions due to CT saturation.Some other reasons can be: Distorted waveforms during high-impedance faults, Noise in the measured signal, especially at low currents (low signal-to-noise ratio), Distorted waveforms, particularly during the first or last cycle of faults.We also observed that the usefulness of synchrophasor data with lower GoF is differentdepending on the reason for GoF being low, and in particular, if the DC offset is the only69Standard PMUGoFApplicationSynchrophasorDC offset modification & calculationGoF evaluationDC offsetGoFActual voltage or current waveformGoFProposed PMUFigure 3.11: Schematic of proposed PMU.reason for GoF being low, the synchrophasor data may still be quite useful (serve the func-tion of interest). This is perhaps intrinsically obvious as the synchrophasor calculationalgorithms (such as Fourier Transforms) are very effective at filtering out the DC compo-nent and identifying the magnitude and phase of the power frequency component of thesignal even in the presence of significant DC. Nonetheless, the generic GoF metric given byequation (3.12) does not discriminate between the cases where GoF is low due to presenceof a DC offset, as opposed to other possible sources stated above.3.5.6 Inclusion of GoF Metrics in the PMUsBased on the discussion above, we propose to include the GoF and its DC-offset-modifiedvariant as prescribed quantities for the PMU to transmit to real-time applications. Inthe proposed structure (marked as the dashed trace in Fig. 3.11), in addition to thephasor information (amplitude, phase, and frequency) reported by the standard PMU, theGoF and GoF (which denotes the GoF-variant that quantifies the mismatch between theoriginal signal and the one reconstructed accounting for DC offset) are also computed andtransmitted. It is worth noting that we do not advocate to transmit the DC offset itselfto the application, because our goal is not to reconstruct the actual voltage and currentwaveforms at the application, but rather to determine whether or not real-time applicationresults are credible. Next, via numerical case studies, we illustrate how the GoF metricscan help to predict the performance quality of the fault-location application.3.6 Case StudiesIn this section, we demonstrate the advantage of incorporating GoF metrics into PMUs byapplying them to the fault-location application. In Example 2, we considered five fault700 10 20 30 40 50 60024681012Fault location error (%)25 30 35 40 45 50 55 6000.20.40.60.81(a) abcg fault in original form0 10 20 30 40 50 60024681012Fault location error (%)25 30 35 40 45 50 55 60Magnified window00.20.40.60.81(b) abcg fault after DC offset remove-10 0 10 20 30 40 50 60024681012Fault location error (%)25 30 35 40 45 50 55 6000.20.40.60.81(c) abg fault in original form-10 0 10 20 30 40 50 60024681012Fault location error (%)25 30 35 40 45 50 55 6000.20.40.60.81(d) abg fault after DC offset remove0 10 20 30 40 50 6001234Fault location error (%)25 30 35 40 45 50 55 6000.10.20.3(e) ag fault in original form0 10 20 30 40 50 6001234Fault location error (%)25 30 35 40 45 50 55 600.10.30.5(f) ag fault after DC offset remove-10 0 10 20 30 40 50 60024681012Fault location error (%)25 30 35 40 45 50 55 6000.20.40.60.81(g) ab fault in original form-10 0 10 20 30 40 50 60024681012Fault location error (%)25 30 35 40 45 50 55 6000.20.40.60.81(h) ab fault after DC offset removeFigure 3.12: GoF versus error for various fault types in original form and after DCoffset removetypes: abcg, abg, ag, abc, and ab. Here, in order to include all fault conditions in oursimulations, we also consider variations in (i) fault resistance, (ii) fault location alongthe line, and (iii) fault initiation time with respect to measurement window boundaries.Particularly, for each fault type, simulations cover fault resistances in the range of 5 −−30 Ω, fault locations of 10, 60, 110, 160, and 210 km from bus m, and 12 fault initiation71times equally spaced within a measurement window. Thus, in total, we conduct 660simulations for each of the five fault types. In each simulation, we assume that the faultlasts for three electrical cycles, which provides three synchrophasors for the during-faultperiod. Each set of synchrophasors (corresponding to a “case”) is used to locate the faultand estimate the fault-to-ground impedance via (3.6).We show the obtained fault location error versus the average GoF of two PMUs mea-suring current of phase a at buses m and n on line (m,n) in Fig. 3.12 in the original formand after DC-offset remove. The obtained results are for various types of fault and dif-ferent measurement windows k = 1, 2, and 3 of the during-fault period are discriminatedwith various colors. It can be observed from Fig. 3.12, the highest fault location errorsare dominantly obtained in measurement window k = 1, while the error rate decreases inmeasurement windows k = 2 and k = 3 for all different types of fault. The results arereasonable considering Tables 3.1 and 3.2. However, the maximum fault location errorvaries based on the fault type, i.e., while ab has the highest, abc and abcg come afterwardsand the minimum values belongs to abg and ag, respectively.Next, we combine the obtained results for all the five fault types to make more generalobservations. But first, we reduce the fault resistances to 5 Ω, 20 Ω, and 30 Ω, and pick4 fault initiation times equally spaced within a measurement window. Thus, we pick 60simulations for each of the five fault types and show them in Fig. 3.13.3.6.1 Simulation Results for Combined CasesSimilar to Fig. 3.12, for each case, we obtain synchrophasor measurements of I(m,n) andI(n,m) at buses m and n, respectively, as well as their corresponding GoF values. Then,for each pair of current measurements obtained at either end of the line, we take theaverage of the two GoF values. As shown in Fig. 3.13a, we subsequently plot the faultlocation error in each case against the average GoF of corresponding measurements. Faultlocation results obtained using synchrophasors from measurement windows k = 1, 2, 3of the during-fault period are marked in red, blue, and green colours, respectively. Viavisual inspection of Fig. 3.13a, we observe that the greater fault location errors (> 1%)are predominantly associated with synchrophasors obtained measurement window k = 1,with low GoF values (< 25 dB). Here, thresholds of 1% and 25 dB are considered tofacilitate subsequent discussions, but they may be tailored for each application scenariobased on the user’s tolerance for accuracy. Moreover, it is worth noting that within eachfault type considered, cases in which the fault initiates 12.5 ms (equivalent to 3/4 of themeasurement window length) into measurement window k = 1 are associated with thehighest fault-location errors (marked as red × in Fig. 3.13a). In these, the PMU does nothave enough information from the during-fault period, so the resulting phasor leads toa poor fit between the time-domain waveform and that reconstructed from the phasor,720 10 20 30 40 50 60024681012Fault location error (%)25 30 35 40 45 50 55 6000.20.40.60.81(a)0 10 20 30 40 50 60024681012Fault location error (%)25 30 35 40 45 50 55 6000.20.40.60.81(b)Figure 3.13: Fault location error versus GoF values computed for waveforms recon-structed from PMU phasor measurements obtained in measurement windowsk = 1, 2, 3 (a) without and (b) with DC-offset modification. The × mark-ers denote cases in which fault initiates 12.5 ms (equivalent to 3/4 of themeasurement window length) into measurement window k = 1.as previously highlighted in Fig. 3.7a. On the other hand, synchrophasors obtained inmeasurement window k = 3 (and some from k = 2) are associated with higher GoF values(greater than 25 dB), which indicate good fit. Corresponding fault-location errors are lessthan 1%. These results echo the phenomenon observed in Example 2 for fault-to-groundimpedance of 5 Ω and fault location of 60 km from bus m.Recall from Example 4 that a simple constant DC-offset modification can significantlyimprove the GoF value for the reconstructed waveform. Thus, we also compute the GoFvalues after DC-offset modification (we call this GoF variant GoF). Similar to Fig. 3.13a,we plot fault location error for each case against its corresponding GoF value in Fig. 3.13b.730 10 20 30 40 50 60020406080100120(a)0 10 20 30 40 50 60050100150200(b)Figure 3.14: Number of cases versus GoF ranges (a) in original form and (b) afterDC-offset remove for simulated cases.Table 3.3: Comparison of fault-location accuracy with corresponding GoF metricsfor an abcg fault occurring at 60 km from bus m with Rg = 5 Ω. Fault initiatesat the beginning of measurement window k = 1.First cycle Second cycle Third cycle% Error GoF GoF % Error GoF GoF % Error GoF GoFWithout saturation 0.6463 6.6909 19.239 0.0438 19.029 31.685 0.1375 32.978 45.808With saturation 21.628 4.0000 4.1050 4.2160 10.150 10.370 1.9305 17.990 18.195Table 3.4: Comparison of fault-location accuracy with corresponding GoF metricsfor an abcg fault occurring at 60 km from bus m with Rg = 5 Ω. Fault initiatesat t = 6.94 ms, within measurement window k = 1.First cycle Second cycle Third cycle% Error GoF GoF % Error GoF GoF % Error GoF GoFWithout saturation 2.1869 1.6994 6.9052 0.3558 12.525 24.041 0.1961 24.574 35.496With saturation 21.015 −0.830 1.5550 8.2592 6.8750 6.8900 3.3985 10.980 11.035740 10 20 30 40 50 6000.511.52Fault location error (%)(a)0 10 20 30 40 50 6000.511.52Fault location error (%)(b)Figure 3.15: Average error and SDE versus GoF ranges (a) in original form and (b)after DC-offset remove for simulated cases .Using this modified GoF metric, many of the cases with acceptable fault location errors(less than 1%) are shifted to the right with GoF values greater than 25 dB. Particularly,in Fig. 3.13a, there are 445 cases with GoF values greater than 25 dB, all of which corre-spond to sufficiently accurate fault location errors of less than 1%. On the other hand,in Fig. 3.13b, 89 more cases with sufficiently accurate fault location results are identifiedbased on the criterion that their corresponding GoF values are greater than 25 dB. Thus,we conclude that indeed, the modified GoF metric is able to better distinguish cases whosefault location results are sufficiently accurate.One can also see that the GoF does not get much better for some red dots as the DCoffset is considered in Fig. 3.13b, most likely because DC offset was not the main reasonfor poor GoF, rather distorted (transitioned) waveform was. For the time window afterthe fault initiation (blue dots in Fig. 3.13), removing DC offset substantially improves GoF,since the main contribution to lower GoF is DC offset (during the fault). The fourth timewindow results are not shown here and they are very similar to the first time window as75they include a distorted waveform as we transition from fault waveform to regular or nosignal (interrupted fault current) during this time (We had fault duration of 3 cycles inall cases).The effect of DC offset modification on shifting the reliable results to the higher rangesof GoF is more obvious in Fig. 3.14 through histograms of the number of cases versus GoFranges for cases shown in Fig. 3.13. A visual comparison of Fig. 3.14a and 3.14b showsthat using GoF one can identify more cases (synchrophasor data during fault) that havehigher GoF and are more suitable for fault location with a high level of accuracy and a highlevel of confidence. For example, Fig. 3.14 shows that using synchrophasor data with GoFor GoF larger than 35 dB will result in fault location accuracy better than 0.3% (within660 m (about 2000 ft) along a 220 km (about 140 mile) transmission line), and Fig. 3.14shows that using GoF (as opposed to the basic GoF) in a software application will yield ahigher number of cases where goodness-of-fit better than 35 dB is available. One can evenchoose GoF or GoF larger than 40 dB to obtain even more accurate (< 0.1% error) faultlocation with a very high level of confidence [90].Finally, in Fig. 3.15, we calculate the average error and the standard deviation of errorin 5 dB steps of GoF. As we can see, not only the average error, but also standard deviationof error (SDE) is less in GoF values higher than 25 dB, either in original form or after DCoffset modification. This fact also proves the reliability of the points falling in this rangefor the application or operator monitoring the fault and location results from the faultand without any certain information about the fault initiation time in a very short timeafter the fault occurrence.Remark 1 (Effect of Current-transformer Saturation). Tables 3.3 and 3.4 report faultlocation accuracy using synchrophasors obtained in measurement windows k = 1, 2, and 3,as well as their corresponding GoF and GoF values. Comparing the rows in Table 3.3 and3.4 we observe that GoF and GoF values are similarly low when CT saturation occurs.This is because the mismatch between the original and reconstructed waveforms is notdue to DC offset, as concluded in Example 6. Even though GoF and GoF values aresimilar, their low values still provide a good predictor for uncredible fault-location results.In measurement windows k = 1, 2, and 3, fault-location errors are quite high, which arereflected in low GoF and GoF values. In fact, the closeness of the two values indicatesthat the reason for low GoF is not DC offset, but the waveform is actually distorted. 3.6.2 GoF Criteria for Predicting Fault Location AccuracyAfter a fault occurs, as shown in Fig. 3.11, synchrophasor measurements and correspondingGoF or GoF values are received by the fault-location application, which utilizes (i) thesynchrophasors to determine the location of the fault, and (ii) the GoF or GoF values to76Fault Location ResultConsiderGoF CredibleInconclusiveGoF < 25 dB GoF ≥ 25 dBGoF < 25 dBCredibleGoF ≥ 25 dBFigure 3.16: Proposed decision tree to assess credibility of fault location.0 0.2 0.4 0.6 0.8 1 1.2020406080(a)0 0.2 0.4 0.6 0.8 1 1.2010203040(b)0 1 2 3 4 5 6 7 8 9 10 11 120204060(c)Figure 3.17: Proportion of cases corresponding to particular ranges of fault locationerrors out of a total of (a) 445 cases with GoF > 25 dB, (b) 89 cases withGoF < 25 dB and GoF > 25 dB, and (c) 366 cases with GoF < 25 dB andGoF < 25 dB.77assess whether or not the fault-location results are sufficiently credible. Based on thediscussion in Section 3.6.1, we note that a GoF value of 25 dB is a reasonable demarcationin this case to determine whether or not fault location results are sufficiently credible.Using the decision process outlined in Fig. 3.16, we first consider cases with GoF >25 dB and conclude that, as shown in Fig. 3.17a, they predominantly correspond to faultlocation errors of less than 0.5%, and 73.03% of those cases correspond to errors of lessthan 0.1%. If GoF is not available because the PMU is not equipped to obtain the DCoffset of time-domain measurement waveforms, the assessment would end here. However,in the case that GoF is available to the application, we note that almost all the caseswith GoF < 25 dB but GoF > 25 dB correspond to fault location errors of less than 1% asshown in Fig. 3.17b, and 86.52% of the cases have errors of less than 0.5%. Finally, asshown in Fig. 3.17c, cases with GoF < 25 dB and GoF < 25 dB have fault location errors ina wide range (as high as nearly 12%). By observing Fig. 3.17, we conclude that cases withGoF > 25 dB or those with GoF > 25 dB return credible fault-location results, but caseswith GoF < 25 dB or GoF < 25 dB are inconclusive. We use term “inconclusive” to avoidinterpretation of the above observations as a linear correlation between the GoF values andthe fault location errors, since there are some cases with small fault location errors, whilealso having small GoF values.3.7 SummaryIn this chapter, we validated the applicability of this metric for fault location applica-tion. Owing to the importance of accurate real-time applications, some metrics have beenbenchmarked to quantify measurement error in both steady and dynamic states, such asTVE, frequency, and ROCOF errors [9]. But, there is still no guarantee to get satisfactoryresults for various applications under these criteria. Therefore, the impact of PMU mea-surement errors on various applications is still under investigation by both researchers andindustry [45, 46]. In this chapter, a new metric called GoF has been applied to illustrate thecredibility of fault location results. Inspired by the fault-location method proposed in [42],we extended it to impedance-based fault location for three-phase circuit for simultaneousFDI. Then, we highlighted the need for GoF metric for real-time applications, specificallywith respect to fault location and show the sensitivity of fault location accuracy withrespect to GoF. Large GoF implies accurate fault location while inaccurate fault locationimplies small GoF. GoF, as a new metric is the only criterion calculating the measurementerror over a time frame, increasing its fidelity as a benchmark for various applications.78Chapter 4Real-time Wide-areaObserver-based Fault Detectionand Identification4.1 BackgroundAs discussed in Chapter 3, the approaches proposed for FDI can be categorized into thefollowing: impedance-based approaches, travelling wave-based approaches, high frequencyapproaches, and finally, AL-based approaches [47]. The method we investigated in Chap-ter 3 is categorized under impedance-based approaches and applies double-end measure-ments for fault location along a particular line. On the other hand, in recent studies, FDIis considered as a system-wide problem [44, 48–52], in which measurements are obtainedfrom various nodes across the system, and the goal is to identify the line that is faulted,as well as the location of the fault along that line.Wide-area FDI methods have received more attention in recent years [44, 48–52, 91], dueto three main reasons. First, PMUs and other digital fault recorders may not be availableat all terminals of the system, which makes the use of the system impedance or admittancematrix necessary to obtain all the voltage and current phasors in the system. Second, thestatus of protective relays and circuit-breakers have been traditionally used for FDI. But,undetected problems may occur in these devices, making protection systems vulnerable tofailures and even leading to power system blackouts. Finally, by expanding the modernsynchrophasor-based WAMS in power systems, a capacity has emerged to accomplish FDIindependently from the protection system. For example, consider the FDI method weinvestigated in Chapter 3. We assumed that the current and voltage synchrophasors areavailable at the both ends of the transmission line. In that problem, the task was to address79the reduced problem of locating a fault along a line with synchrophasor measurementsavailable at both ends; the diagnosis of the faulted line was not needed. However, thismay not be the case when we look at the problem in a large-scale system, that is, thefaulted line must be identified prior to location of the fault along that line. Also, if PMUswere not installed at the both ends of a line of interest, the method proposed in Chapter 3would not be applicable.Different wide-area FDI methods have been proposed in the literature. A wide-area FDImethod is proposed in [44] based on sparse estimation technique. The authors reformu-late the fault location as a sparse estimation problem. The proposed method is tested onrelatively large systems and is capable of identifying the fault type, even in teed-circuit con-nections. But the proposed sparse estimation technique is based on iterative calculations,which is burdensome in case of large power systems and in presence of noise in the mea-surements. Therefore, the authors propose to decompose the sparse estimation probleminto sub-problems as future work. The method proposed in [49], uses the bus-impedancematrix and measurements obtained from various buses in the power system. However,it does not discuss the robustness of the proposed FDI method against deficiencies in atypical WAMS and its PMUs, such as bad data, misshing values, and measurement errors.In [50], authors match the synchrophasors obtained from the WAMS with predefined faultlocation factors for every bus, computed based on the fault distance from the consideredbus. Next, the homogeneity of the proposed factor is monitored in the during-fault periodto identify the fault region. Finally, the exact location of the fault is determined. But theproposed method uses an interative two-stage algorithm and investigates many buses andlines at each stage, which makes it unsuitable for real-time applications. In [52], systemimpedance matrix is used along some dispersed PMUs for wide-area FDI. The proposedmethod is hierarchical, i.e., first fault inception is detected, followed by diagnosis of faultzone and faulted line. But this three-stage method would be computationally burdensomefor a real-time application. The method presented in [51] is based on solving an over-determined system to obtain the best fit of current sources for the faulted lines. Here,both positive- and negative-sequence networks are analyzed to find the fault location viaan optimization problem, satisfying equality constraints between both networks. However,this method is not tested on large cases.In this chapter, we model the power system by state-space equations which are validin both pre-fault and during-fault periods and consider the faults occurring in the systemas sudden changes in the input of this model. We then design an observer-based FDIfilter to simultaneously identify the faulted line and the fault location along that line.Observers are commonly used for dynamic state estimation in power systems [92, 93].The application of observers, such as linear quadratic estimators and the Kalman filter,has made it possible to use real-time simulators and feed measurements obtained from a80system to estimate fast power system dynamics [94, 95]. In light of this, we propose awide-area FDI method based on the network admittance matrix and the linear Luenbergerobserver. The proposed method simultaneously identifies and locates the fault by ensuringthe error residual between the actual system and the observer dynamics. Moreover, theproposed method uses sequentially sampled measurements over time, instead of just oneset of measurements from each of pre-fault and during-fault scenarios, as in the case inexsiting methods. Consequently, it is more robust against false alarms and misdetections.The proposed method also has the potential to be extended to three-phase system modelsand identify various types of faults. However, in the remainder of this chapter, we restrictour models and simulations to three-phase-to-ground faults in positive-sequence networksusing per-phase analysis.4.2 PreliminariesSuppose PMU phasor measurements are reported at intervals of ∆t, i.e., at times tk =k∆t, k = 1, 2, . . . Here, we assume each sampling interval is one electrical cycle, i.e.,∆t = 16.666 ms. For an N -bus power system in pre-fault period, let I[k] ∈ CN andV [k] ∈ CN denote system terminals’ current and voltage vectors, respectively, sampled attk = k∆t. Also, define ∆V [k] = V [k + 1]− V [k] and ∆I[k] = I[k + 1]− I[k]. Via Ohm’slaw, we can write∆I[k] = Y∆V [k], (4.1)where Y ∈ CN×N is the network admittance matrix.Next, suppose a three-phase-to-ground fault occurs between time instants kf and kf+1along line (m,n) with admittance ymn and length `, at distance d away from bus m,as shown in Fig. 4.1. The fault divides the line into two sections connected to busesm and n with admittances ymmn =ymnx , and ynmn =ymn1−x , respectively, where x =dl .Let If [k] and Vf [k] denote the fault current and voltage at time instant k. Further, let∆If [k] = If [k+1]−If [k] and ∆Vf [k] = Vf [k+1]−Vf [k]. Then, similar to 4.1, by applyingthe Ohm’s law to the faulted network and considering the fault location as bus N + 1, wecan write[∆I[k]∆If [k]]=[Y˜ yyT ŷ][∆V [k]∆Vf [k]], k ≥ kf . (4.2)In (4.2), entries of Y˜ are the same as Y except1 [Y˜ ]mn = 0, [Y˜ ]nm = 0, [Y˜ ]mm =1We indicate the entry of Y˜ in row m and column n by [Y˜ ]mn.81yny mFigure 4.1: General diagram of fault.[Y ]mm − ymn + ymmn, and [Y˜ ]nn = [Y ]nn − ymn + ynmn since there is no direct connectionbetween buses m and n in the faulted network. Furthermore, entries of y are 0 exceptym = −ymmn and yn = −ynmn. Notably, the locations of the nonzero entries of y correspondto the buses at either end of the faulted line. Finally, ŷ = ymmn + ynmn.We can write the second equation in (4.2) as∆If [k] = yT∆V [k] + ŷ∆Vf [k], (4.3)from which we solve∆Vf [k] =1ŷ(∆If [k]− yT∆V [k]). (4.4)Substituting (4.4) into the first equation in (4.2), we obtain∆I[k] = Y˜∆V [k] + y∆Vf [k]= Y˜∆V [k] + y∆If [k]ŷ− 1ŷyyT∆V [k]=(Y˜ − 1ŷyyT)∆V [k] +∆If [k]ŷy. (4.5)As mentioned above, entries of Y˜ and Y are the same, except [Y˜ ]mn = 0, [Y˜ ]nm = 0,[Y˜ ]mm = [Y ]mm − ymn + ymmn, and [Y˜ ]nn = [Y ]nn − ymn + ynmn. On the other hand, theentries of the matrix 1ŷyyT are zero except:[1ŷyyT]mm=x(1− x)ymny2mnx2=(1− x)xymn, (4.6)[1ŷyyT]mn=x(1− x)ymny2mn(1− x)x = ymn, (4.7)[1ŷyyT]nm=x(1− x)ymny2mn(1− x)x = ymn, (4.8)82[1ŷyyT]nn=x(1− x)ymny2mn(1− x)2 =x(1− x)ymn. (4.9)Therefore, in (4.5), relevant entries in Y˜ − 1ŷyyT can be simplified as follows:[Y˜ ]mm −[1ŷyyT]mm= [Y ]mm − ymn + ymmn −(1− x)xymn = [Y ]mm, (4.10)[Y˜ ]mn −[1ŷyyT]mn= 0− ymn = [Y ]mn, (4.11)[Y˜ ]nm −[1ŷyyT]nm= 0− ymn = [Y ]nm, (4.12)[Y˜ ]nn −[1ŷyyT]nn= [Y ]nn − ymn + ynmn −x(1− x)ymn = [Y ]nn, (4.13)and soY˜ − 1ŷyyT = Y. (4.14)Substitution of (4.14) into (4.5) results in∆I[k] = Y∆V [k] +∆If [k]ŷy. (4.15)Since ŷ = ymmn + ynmn =ymnx +ymn1−x =ymnx(1−x) , for the during-fault period, we can simplify(4.15) as∆I[k] = Y∆V [k] +∆If [k]ymnx(1− x)y, k ≥ kf . (4.16)Finally, by unwrapping (4.1), we arrive at the following recurrence relation that de-scribes how the current injections evolve in the pre-fault period,I[k + 1] = I[k] + Y∆V [k], k < kf . (4.17)Similarly, by unwrapping (4.16) in the during-fault period, we getI[k + 1] = I[k] + Y∆V [k] +∆If [k]ymnx(1− x)y, k ≥ kf . (4.18)In summary, in the pre-fault period, current dynamics are governed by (4.17), and inthe during-fault period, by (4.18) with nonzero ∆If [k]. We note that (4.18) differs from(4.17) by one additional term. This vector quantity has nonzero values only in entriescorresponding to the two buses on either end of the faulted line. Furthermore, embedded83in these nonzero entries are x, which represents the fault location along the faulted line,and If , which denotes the fault current. Our goal is to design an observer-based faultdetection and identification filter that takes measurements ∆V [k] and ∆I[k] as inputs andproduces a residual e[k] with the following properties:(i) When no fault has occurred, i.e., ∀ k < kf , ∆If [k] = 0, the filter residual asymptot-ically converges to zero, i.e., limk→∞ e[k] = 0.(ii) When a fault occurs between time instants kf and kf + 1, i.e., for some fault alongline (m,n), ∆If [k] 6= 0, k ≥ kf , then as k → ∞, the filter residual tends to alignwith y, the nonzero entries of which help to identify the faulted line.(iii) Once the faulted line is identified via the nonzero entries of the filter residual in (ii)above, the magnitudes of the nonzero entries help to solve for the location of thefault along the identified line, as well as the fault current.If the filter design satisfies these properties, then fault detection can be achieved byobserving the filter residual e[k]. Furthermore, if each fault under consideration has aunique signature, then that fault can be identified by the observing the residual, as wellas the fault current.4.3 Proposed ObserverWith the desired properties listed above, in this section, we propose the following observeras a candidate for the FDI filter:Î[k + 1] = Î[k] + Y∆V [k] + L(I[k]− Î[k]), (4.19)where Î[k] denotes the estimated current at time instant k and L = (1−µ)diag(1N ), with−1 < µ < 1 and 1N denotes the vector of 1’s of N .Let e[k] = I[k] − Î[k] denote the error residual between the system and observerdynamics. Then, if no fault has occurred, for all k > 0, the residual dynamics computedwith the pre-fault system dynamics in (4.17) followse[k + 1] = µe[k]. (4.20)The solution to the difference equation in (4.20) ise[k] = µke[0], (4.21)and since −1 < µ < 1, it follows that limk→∞ e[k] = 0, as desired.84If, on the other hand, a three-phase-to-ground fault occurs along line (m,n) betweentime instants kf and kf + 1, then the residual dynamics computed with the during-faultsystem dynamics in (4.18) followe[k + 1] = µe[k] + ∆If [k]x(1− x)ymny, k ≥ kf , (4.22)the closed-form solution of which ise[k] = µke[kf ] +k−1∑r=0µk−r−1∆If [r]x(1− x)ymny, k ≥ kf . (4.23)Since the first term in (4.23) vanishes as k →∞, e[k] tends to align with y, as desired.Here, we made two key observations. First, with appropriate choice of µ (i.e., smalland positive), we find via numerical case studies that the first term vanishes with a fewmeasurement samples. Moreover, the magnitude of the mth and nth entries of e[k] arefunctions of x (fault location) and ∆If [k] (fault current). With this in mind, at eachtime instant in the during-fault period, we can solve for unknowns x and ∆If [k] from thefollowing equations:em[k + 1] = µem[k] + ∆If [k](1− x), (4.24)en[k + 1] = µen[k] + ∆If [k]x, (4.25)where em[k] and en[k] denote the mth and nth entries of e[k] (similarly for em[k + 1] anden[k + 1] in e[k + 1]), which is the error residual obtained between the actual-system andobserver dynamics.4.4 Case StudiesIn this section, we validate the proposed method for wide-area FDI via case studies in-volving a two-bus, a three-bus, and the canonical two-area systems. All simulations areconducted using the DSATools package [74], and all imposed faults are three-phase-to-ground (abcg) type. The proposed method can be extended to identify other fault types;and we reserve this for future work. The simulations generate synthetic voltage- andcurrent-phasor measurements once per electrical cycle, which are used as inputs to theproposed FDI filter. In practical implementation, these measurements would be availablethrough WAMS.85G1ymnI1V1P2+jQ2V2I2Figure 4.2: Two-bus system used for testing the proposed method.4.4.1 Two-bus SystemWe begin with the simplest system, i.e., a two-bus system, as depicted in Fig. 4.2. Thesystem is composed of a generator at bus 1, serving a load at bus 2, with P2 = 9.67 p.u. andQ2 = −1.06 p.u., using base power and voltage of 100 MVA and 230 kV, respectively. Theline connecting buses 1 and 2 has impedance 0.002 + j0.02 p.u. The model chosen for thegenerator is a synchronous machine with transient open-circuit time constant T′do = 8 s,T′qo = 0.4 s, sub-transient open-circuit time constant T′′do = 0.03 s, T′′qo = 0.05 s, inertiaconstant H = 6.5 s, synchronous reactance Xd = 1.8 p.u., Xq = 1.7 p.u., transientreactance X′d = 0.3 p.u., X′q = 0.55 p.u., sub-transient reactance X′′d = 0.25 p.u., andfinally, leakage reactance Xl = 0.2 p.u. The applied exciter model is ESAC4A [96], whileGOV6 is selected for the governor model, and finally, the power system stabilizer modelis chosen PSS1. Both of the latter models can be found in the model manual of TSATtool [74]. The parameters of the models and their values are listed in Appendix A.We impose an abcg fault to the line with the length ` at three different locations: at0.25`, 0.5` , and 0.75` from bus 1. The fault initiates just after time instant kf = 6. Theamplitude of the residual current, e[k], at buses 1 and 2 for an abcg fault at 0.25` is shownin Fig. 4.3. Two different values of µ, i.e., 0.90, and 0.01, are chosen. As can be observedin Figs. 4.3a and 4.3b, the amplitude of the residual current phasor is zero in the pre-faultperiod for both buses 1 and 2. When the fault initiates, the residual current increasessharply, since the dynamics of the system and the observer do not match. However, aftera few time instants, the actual and the observer currents converge. The convergence ratedepends on the µ value. Comparing Figs. 4.3a and 4.3b, we can observe when a largevalue of µ is selected in the pre-defined range of −1 < µ < 1, i.e., µ = 0.9, the convergencerate is much slower than the case of small value µ = 0.01, in which the actual and theobserved currents converge to each other in a few time instants. This is because that thesmaller value of µ leads to a faster vanishing of the first term in (4.23). We also notethat choosing µ values from the range µ < 0 results in the same trend for convergence,while seeing fluctuations, i.e., overshoots and undershoots, in the observer current, sincethe first term in (4.23) is alternatively positive and negative in consecutive samples.861 2 3 4 5 6 7 8 9 10 11 12051015202530Amplitude (pu)(a) µ = 0.91 2 3 4 5 6 7 8 9 10 11 12051015202530Amplitude (pu)(b) µ = 0.01Figure 4.3: Amplitude of residual current (e[k]) for two different values of µ , i.e.,µ = 0.9, and µ = 0.01, for an abcg fault imposed at 0.25` of line (1,2) withthe length of ` in the two-bus system. The fault initiates immediately aftertime instant kf = 6 and is denoted by a bold dashed line.Finally, the obtained fault location and fault current results are shown in Tables 4.1and 4.2. The fault location x and the fault current changes at time instant k, i.e., ∆If [k]are obtained through (4.23) based on the actual voltage and current phasors and also,observer current phasors at earlier time instants. We note that the obtained values inthe pre-fault period are complex which do not make sense as distance quantities. Thesevalues are originally the result of dividing two very small residual current values obtainedat time instants in the pre-fault period and thus, should not be considered. But as soonthe fault initiates, the results change to real values which are correct for various locationson the transmission line. We also note in Table 4.2, that obtained fault current change inthe during-fault period increases sharply, as expected. These obtained quantities matchthe actual fault current values obtained in the system and can be used for more faultinspections in future work.87Table 4.1: Fault location results for different locations along the line (1, 2). Theresults obtained in the pre-fault period (k = 4–6) are complex which do notmake sense for a distance value, while those obtained from the during-faultperiod (k = 7–9) are real and correct.Location Obtained fault locations (x)of imposed Pre-fault period During-fault periodfault k = 4 k = 5 k = 6 k = 7 k = 8 k = 9x = 0.25 0.1729 + j0.0436 0.2965− j0.0108 −0.0478 + j0.2147 0.2500 0.2500 0.2500x = 0.50 0.2328− j0.0040 0.7861 + j0.1135 0.6031− j0.0176 0.5000 0.5000 0.5000x = 0.75 0.4843− j0.1059 0.9273 + j0.0410 0.8903 + j0.0944 0.7500 0.7500 0.7501Table 4.2: Fault current for different locations along the line (1, 2). The resultsobtained in the pre-fault period (k = 6) are almost zero, while those obtainedin the during-fault period (k = 7–9) are large values.Location Obtained fault currents (∆If [k])of imposed Pre-fault period During-fault periodfault k = 6 k = 7 k = 8 k = 9x = 0.25 1.3× 10−6 − j4.4× 10−6 −11.546 + j25.872 −1.8709− j3.4608 −1.4214− j2.2134x = 0.50 1.5× 10−6 − j1.6× 10−5 −10.099 + j22.918 −1.5863− j2.8004 −1.2562− j1.8572x = 0.75 −1.5× 10−6 − j1.6× 10−5 −8.9755 + j20.562 −1.3731− j2.3200 −1.1225− j1.58514.4.2 Three-bus SystemIn this section, we test the proposed method on a three-bus system shown in Fig. 4.4. Inthis system, faults can imposed to different lines in the system to check the capability ofthe proposed method for identification of faulted line, which is expected from a wide-areaFDI.The system is composed of two generators at buses 1 and 2, serving a load at bus 3,with P3 = 9.67 p.u. and Q2 = −1.06 p.u., using base power and voltage of 100 MVA and230 kV, similar to the previous two-bus system. Also, all the three lines in the systemconnecting the three buses of the system have the same impedance 0.001 + j0.01 p.u. Thepower flow results in power generations as P1 = 3.26 p.u., Q1 = 1.60 p.u., P2 = 6.50 p.u.,and Q2 = −1.67 p.u. The models of the two generators and their control systems aresimilar to the generator of the two-bus system in Fig. 4.2.We impose faults at different locations along the three different lines in the system.The obtained results are shown in Table 4.3 for the set of phasors obtained on differentlines2. Following the fault imposed on a line, the non-zero entries of matrix y correspondto the buses at either end of the faulted line. Similar to the two-bus system, the obtainedresults of fault location are correct after the fault is initiated, changing from complex toreal correct values. The obtained fault location currents are shown in Table 4.4 and matchthe actual current values with a good precision.2In this chapter, all the fault locations are considered from the smaller bus number. Also, the smallercomes first in our notations.88y12y13 y23V3P3+jQ3G1 G2V1 V2I3I1 I2Figure 4.4: Three-bus system used for testing the proposed method.Table 4.3: Fault location results for different lines and locations in the three-bussystem. The results obtained in the pre-fault period (k = 4–6) are complex,which do not make sense for a distance value, while those obtained from theduring-fault period (k = 7–9) are correct.Location Obtained fault locations (x)of imposed Pre-fault period During-fault periodfault k = 4 k = 5 k = 6 k = 7 k = 8 k = 9line (1, 2), x = 0.50 0.2770 + j0.0246 −0.5429 + j0.7471 0.4262 + j0.0059 0.5000 0.5000 0.5000line (1, 3), x = 0.75 1.2570 + j0.0876 0.9019 + j0.0131 0.8047 + j0.0124 0.7500 0.7500 0.7500line (2, 3), x = 0.25 −0.1799 + j0.2623 2.5831− j0.0061 0.3250− j0.0226 0.2500 0.2500 0.2500Table 4.4: Fault location results for different lines and locations in the three-bussystem. The results obtained in the pre-fault period (k = 6) are almost zero,while those obtained in the during-fault period (k = 7–9) are large values.Location Obtained fault currents (∆If [k])of imposed Pre-fault period During-fault periodfault k = 6 k = 7 k = 8 k = 9line (1, 2), x = 0.50 5.8× 10−7 + j3.4× 10−6 −11.120 + j46.599 −2.4991− j4.0238 −1.8530− j2.4569line (1, 3), x = 0.75 −5.6× 10−7 − j8.6× 10−7 −9.8988 + j38.584 −1.8850− j2.7972 −1.5213− j1.7864line (2, 3), x = 0.25 3.0× 10−6 + j5.9× 10−6 −11.557 + j43.540 −2.3184− j4.0772 −1.7709− j2.47914.4.3 Canonical Two-area SystemThe final test system under study is the canonical two-area test power system [5], whichis shown in Fig. 4.5. The only adjustment made here is we replace the multi-circuittransmission lines with an equivalent single line since the proposed method is not able todiscern faults in multi-circuit transmission lines. The system includes two loads at buses7 and 9 with P7 = 9.67 p.u., Q7 = −0.96 p.u., P9 = 17.67 p.u., and Q9 = −2.46 p.u.,using base power and voltage of 100 MVA and 230 kV. Again, all the lines in the systemare considered to have the same impedance 0.002 + j0.02 p.u.. Following the power flow,the generators generate power as P1 = 6.85 p.u., Q1 = 1.31 p.u., P2 = 7.00 p.u., Q2 =0.58 p.u., P3 = 7.19 p.u., Q3 = 1.26 p.u., P4 = 7.19 p.u., and Q4 = −0.46 p.u. The models89G1G2 G4G3125 67 8 910 11 34P7+jQ7 P8+jQ8Figure 4.5: Two-area system used for testing the proposed method.Table 4.5: Fault location results for different lines and locations in the canonicaltwo-area system. The results obtained in the pre-fault period (k = 4–6) arecomplex, while those obtained from the during-fault period (k = 7–9) arecorrect.Location Obtained fault locations (x)of imposed Pre-fault period During-fault periodfault k = 4 k = 5 k = 6 k = 7 k = 8 k = 9line (7, 8), x = 0.50 0.1493− j0.0025 0.6958− j0.0413 0.7206− j0.0413 0.5000 0.5000 0.5000line (8, 9), x = 0.25 0.6323 + j1.2599 0.4274 + j0.7154 0.2863 + j0.0707 0.2500 0.2500 0.2500line (9, 10), x = 0.75 0.8132− j0.3713 0.9358 + j0.0900 0.8745 + j0.0345 0.7500 0.7500 0.7500of all the generators and their control systems are similar to the generator of the two-bussystem in Fig. 4.2.Then, faults are imposed on three different locations along three lines in the system toassess whether or not the proposed method is capable of FDI in a multi-bus system. Theresults are again promising in the during-fault period as shown in Table 4.5. The obtainedchanges in fault current are shown in Table 4.6. The obtained values match the actualvalues, similar to the previous case studies.Table 4.6: Fault location results for different lines and locations in the canonicaltwo-area system. The results obtained in the pre-fault period (k = 6) arealmost zero, while those obtained in the during-fault period (k = 7–9) arelarge values.Location Obtained fault currents (∆If [k])of imposed Pre-fault period During-fault periodfault k = 6 k = 7 k = 8 k = 9line (7, 8), x = 0.50 −8.5× 10−6 + j8.5× 10−6 −7.5864 + j48.151 −2.8073− j2.2444 −2.4356− j1.5677line (8, 9), x = 0.25 2.3× 10−6 − j3.2× 10−5 −7.5244 + j47.956 −2.8436− j2.1953 −2.4716− j1.5462line (9, 10), x = 0.75 6.5× 10−6 + j5.8× 10−5 −9.6975 + j56.583 −3.6548− j3.4750 −2.9965− j2.3006904.5 SummaryIn this chapter, we proposed a wide-area FDI method based on estimation of fault currentthrough observers. This method is verified to locate the fault correctly in the case of havingfull observability on the system, i.e., having current- and voltage-phasor measurementsfrom all the system buses. As an immediate direction for future work, we plan to modifythe proposed method to accommodate cases in which measurements are not available fromall buses in the system.91Chapter 5Concluding RemarksIn this chapter, we provide a summary of the results presented in Chapters 2–4 of thethesis including the contributions. Then, we discuss the future work for each topic indetail.5.1 Thesis Summary and ContributionsA summary and contributions of each chapter are provided below.5.1.1 Real-time TSA using Synchrophasor-based WAMSIn Chapter 2, we focused on the practical implementation of classification and regressionmodels to assess transient stability in real time using synthesized PMU measurements ob-tained from the system. Through these models, we showed the possibility of predictionof stable and unstbale cases in the first few cycles of the post-contingency period with ahigh accuracy, which shows a superiority of WAMS over traditional SCADA systems. Weassumed that we have limited number of PMUs in the BCH power system at some spe-cific locations such as large power plants and interconnections with neighbouring utilites.This is in accordance with the real infrastrcuture available in the system. We also assumedthat the synthesized synchrophasors are available at the control centre with the rate of onesample per electrical cycle. Having these conditions in mind, first, we identified systemcharacteristics that are pertinent to TSA, including voltage angles, active power gener-ation of major power plants, deviations from COI angle and speed, and potential- andkinetic-energy related quantities. Next, these characteristic quantities, called “indices”or “features”, were evaluated using synthetic measurements obtained from time-domainsimulations of the full WECC system. These evaluated quantities were used to build classi-fication and regression models, which were then tested and verified on the WECC system,with emphasis on transient stability issues in the BCH power system that may emerge due92to high power exchange rates with the BPA system. The classification model was devel-oped using CART, which is a well-known DT algorithm [31]. DTs have been selected forclassification in this study because of their transparency, and simplicity for interpretationwhich can help for observing important features for classification. We trained five CARTmodels based on the measurements obtained sequentially in the first five cycles of the post-contingency period. This can be interpreted as an ensemble of CART models which forma voting system and makes a more reliable overall decision. Each CART model was ableto classify stable and unstable cases with 99% accuracy rate and showed the priority ofCOI- and energy function-based indices over the pure synthesized measurements obtainedfrom the PMUs. In addition to CART, which returns discrete-valued stable versus unstableclassifications, we proposed to use MARS [32] in conjunction with a continuous-valued TSIto assess the severity of the contingency and the system’s proximity to instability thatcan be useful to determine probable corrective actions. This regression method is alsotransparent. Similar to CART, we use an ensemble of MARS models over the time spanof the first five sequential cycles of the post-contingency period to decrease the error ofthe final predicted TSI values. The obtained MAE errors can decrease to around 2, whichshows the reliability of the trained models. We also showed the necessity and practical-ity of updating the obtained models after topological changes or in the presence of highpenetration of renewable resources, which is called “Adaptive Training”.5.1.2 Real-time Impedance-based FDI Using Synchrophasor-basedWAMSIn Chapter 3, we focused on another real-time application feasible based on WAMS, i.e.,FDI. We illustrated the usability of GoF metrics to quantify synchrophasor accuracy inreal-time applications. Specifically, for the applications that use synchrophasors obtainedin transient state, such as fault location along a transmission line, GoF metrics indeed yieldan assessment for the confidence level of results. The GoF metric is calculated over thetime frame of a measurement window in a PMU and thus, can show the match betweenthe original signal and the one reconstructed from the resultant synchrophsor values. TheGoF introduced in [82] and its variant introduced in this chapter (GoF) are two potentialmetrics for synchrophasor data that can help synchrophasor applications provide credibil-ity indications for their end results. Through numerous examples, we demonstrated theuse of the GoF and (GoF) metrics to quantify the accuracy of synchrophasors in differenttransient periods, i.e., for various fault initiation times, and also in presence of CT satu-ration. Next, we generated a big data set composed of simulation cases for various faulttypes, fault initiation times, fault resistances, and fault locations, and showed that there isa correlation between high values of GoF metrics and the fault location error for a specificsynchrophasor. Actually, smaller fault location errors are obtained from synchrophsors93with larger GoF values. Therefore, based on GoF metrics, the results of synchrophasor-based fault location applications can be evaluated credibly. Finally, we showed that theproposed metrics can be used for a practical fault location application at a power utility.For example, specific GoF thresholds and limits can be set based on the accuracy andconfidence level desired for the fault location results in a specific system.5.1.3 Real-time Wide-area Observer-based FDIFinally, we proposed a wide-area fault location method using a linear observer in an effortto locate the fault in a power system without having full observability through PMUs.Wide-area FDI is a real-time application of WAMS that has received attention recently. Inthis method, the faulted line is detected prior to the fault location, independent of relays’operation, which can be a back up for local protection systems. Observers have beenpreviously used for state estimation. Here, we propose to use them for FDI. Therefore,we modeled the problem of fault location in state space to be able to use state estimationtechniques for fault location. The input states are voltage and current phasors and assumedto be obtained once per cycle from WAMS. We calculated the net injected actual andestimated current phasors at each bus through the admittance matrix and voltage andcurrent phasors from the previous step and then, obtained their difference called residualerror. This residual error can be used to detect the faulted line and also, calculate the faultlocation and the fault current. Therefore, the proposed method is capable of detectingthe faulted line in addition to locating the fault in real time. We have shown this throughthree different examples on a two-bus, a three-bus, and a canonical two-area power system.As soon as the fault initiates in the system, the residual error for the buses at the ends ofthe faulted line increase. These error values are used for fault location along the faultedline. The obtained results are promising.5.2 Future WorkInspired by the results and contributions of Chapters 2–4, here we discuss suggested ideasfor future work in the following sections.5.2.1 Robustness of the Developed Models for TSA againstMeasurement Errors, Bad Data and Missing ValuesMost studies have assumed that WAMS provide reliable and accurate synchrophasors forreal-time TSA. However, in real applications, synchrophasors may become unavailabledue to unexpected failures in PMUs, communication infrastructure, and PDCs. Also, intransient states, many measurements may not be as accurate as expected. Further, thedelivery of synchrophasors from remote locations of power systems to control centres can94experience high latency when communication networks are heavily congested, which alsomay result in the unavailability of synchrophasors. Therefore, it is necessary to developclassification and regression models which are robust against bad data and missing valuesof synchrophasors. On the other hand, if the post-fault synchrophasors are applied toclassification or regression models, the problem of corrective or emergency control canbe tackled [72, 97, 98]. Based on the classification and especially regression results, anappropriate emergency control action can be initiated. Therefore, another future work forthis part is investigating the possible corrective control actions considering the selectedfeatures in the classification or regression stage, and also the severity of instability issue,with the purpose of mitigating unstable or marginally stable cases. The predicted TSIvalues can be very useful for this purpose. Different approaches can be used for correctivecontrol such, as load shedding and generation shedding.5.2.2 Incorporating GoF Metrics to Improve Synchrophasor-basedApplicationsGoF can also be important when an undetected high-impedance fault results in distortedwaveforms and inaccurate synchrophasor measurements. The presence of noise can addi-tionally affect measurement signals. Therefore, we see that there are even more examplesin which the importance of GoF is clear. Furthermore, the same approach that we used forreal-time fault location applications can be extended to a number of other synchropha-sor applications. For example, GoF can be used to qualify synchrophasor data that maybe used in a dynamic state estimation (or a linear state estimation), special protectionschemes, transient stability assessment, and dynamic state estimation. Specifically, a sen-sitivity study can be implemented on the CART and MARS models developed for real-timeTSA.5.2.3 Wide-area Fault Location with Limited Number of PMUsThe wide-area FDI method proposed in Chapter 4 needs to be modified to solve specificchallenges. First, in real power systems, only some of the buses are equipped with PMUs.Therefore, this method should be modified in order to operate correctly without full ob-servability on the system via measurements. One of the well-known techniques used for“network-reduced” models of a power system is the Kron reduction [99]. The second chal-lenge of the proposed method is its capability to identify different types of short-circuitfaults. The results shown in Chapter 4 were obtained based on positive-sequence sim-ulations. One solution for this challenge is to apply the proposed method for all threephases separately to investigate if the fault is imposed to a specific phase. Although thiswould introduce a larger computational burden, the process ensures correct identification95of the fault type. Finally, a challenge in FDI arises in multi-circuit transmission lines whenimpedance or an admittance matrix is used, since the admittance matrix considers onesingle equivalent value for the admittance of parallel lines. To overcome this challenge, anode breaker model can be used. In common system impedance or admittance matrices,bus-branch information is employed, in which it is assumed that all lines at one substationare connected to each other, ignoring the breakers. Considering a node-breaker model, asimple bus can be divided into multiple elements (buses), connected via zero impedancebranches. Through this technique, the topology of the system can be reflected in theimpedance or admittance matrix [100, 101].96Bibliography[1] North American Electric Reliability Corporation (NERC). (2015, Nov.) NERCoperation time horizons. [Online]. Available: http://www.nerc.com/files/TimeHorizons.pdf → pages iii[2] F. Rahmatian, D. Atanakovic, M. Rahmatian, G. Stanciulescu, and M. Nagpal,“BC Hydro synchrophasor system for wide area monitoring, protection and controlfunctional requirements and system architecture considerations,” in Proceeding ofCIGRE´ Canada Conference, 16–18 Oct., 2016. → pages xiv, 2, 3, 4, 5, 6, 8[3] M. Pavella, D. Ernst, and D. Ruiz-Vega, Transient stability of power systems: aunified approach to assessment and control. Springer Science & Business Media,2012. → pages xiv, 16, 17, 18[4] Western Electricity Coordinating Council (WECC). (2013, Dec.) WECC 2014operation case map. [Online]. Available: https://www.wecc.biz/ → pages xiv, 35[5] P. Kundur, N. J. Balu, and M. G. Lauby, Power system stability and control.McGraw-hill New York, 1994. → pages xv, 16, 25, 54, 89[6] E. Vaahedi, Practical Power System Operation. John Wiley & Sons, 2014. → pages1[7] B. Liscouski and W. Elliot, “Final report on the august 14, 2003 blackout in theunited states and canada: Causes and recommendations,” A report to US Depart-ment of Energy, vol. 40, no. 4, 2004. → pages 2[8] A. G. Phadke and J. S. Thorp, Synchronized phasor measurements and their appli-cations. Springer Science & Business Media, 2008. → pages 2, 6[9] “IEEE Standard for Synchrophasor Measurements for Power Systems,” IEEE StdC37.118.1-2011 (Revision of IEEE Std C37.118-2005), Dec. 2011. → pages 3, 5, 6,14, 48, 50, 51, 59, 7897[10] “IEEE Standard for Synchrophasor Measurements for Power Systems – Amendment1: Modification of Selected Performance Requirements,” IEEE Std C37.118.1a-2014(Amendment to IEEE Std C37.118.1-2011), Apr. 2014. → pages 3, 6, 36, 48, 50, 51[11] “IEEE Standard for Synchrophasor Data Transfer for Power Systems,” IEEE StdC37.118.2-2011 (Revision of IEEE Std C37.118-2005), Dec. 2011. → pages 3, 5, 6,48[12] “IEEE Guide for Phasor Data Concentrator Requirements for Power System Pro-tection, Control, and Monitoring,” IEEE Std C37.244-2013, May 2013. → pages 3,7[13] M. Kezunovic, S. Meliopoulos, V. Venkatasubramanian, and V. Vittal, Applica-tion of Time-synchronized Measurements in Power System Transmission Networks.Springer, 2014. → pages 3, 5, 9, 10[14] B. Naduvathuparambil, M. C. Valenti, and A. Feliachi, “Communication delays inwide area measurement systems,” in Proceedings of the Thirty-Fourth SoutheasternSymposium on System Theory (Cat. No.02EX540), 2002, pp. 118–122. → pages 6[15] M. Rahmatian, Y. C. Chen, W. Dunford, and A. Moshref, “Transient stability as-sessment of bc hydro power system using decision trees,” in Proceeding of CIGRE´Canada Conference, 16–18 Oct., 2016. → pages 8, 10[16] V. Madani, J. Giri, D. Kosterev, D. Novosel, and D. Brancaccio, “Challenging chang-ing landscapes: Implementing synchrophasor technology in grid operations in thewecc region,” IEEE Power and Energy Magazine, vol. 13, no. 5, pp. 18–28, 2015. →pages 8, 9, 10, 11, 13[17] V. Madani, D. Novosel, S. Horowitz, M. Adamiak, J. Amantegui, D. Karlsson,S. Imai, and A. Apostolov, “IEEE PSRC report on global industry experiences withsystem integrity protection schemes (sips),” IEEE Transactions on Power Delivery,vol. 25, no. 4, pp. 2143–2155, 2010. → pages 8, 10, 11[18] Z. Huang, P. Du, D. Kosterev, and S. Yang, “Generator dynamic model validationand parameter calibration using phasor measurements at the point of connection,”IEEE Transactions on Power Systems, vol. 28, no. 2, pp. 1939–1949, 2013. → pages9, 11[19] M. Rahmatian, Y. C. Chen, W. G. Dunford, and F. Rahmatian, “Incorporatinggoodness-of-fit metrics to improve synchrophasor-based fault location,” IEEE Trans-actions on Power Delivery, vol. PP, no. 99, pp. 1–1, Jan. 2018. → pages 1198[20] S. Zadkhast, J. Jatskevich, and E. Vaahedi, “A multi-decomposition approach foraccelerated time-domain simulation of transient stability problems,” IEEE Transac-tions on Power Systems, vol. 30, no. 5, pp. 2301–2311, Sept 2015. → pages 11[21] H. Chiang, Direct methods for stability analysis of electric power systems: theoreticalfoundation, BCU methodologies, and applications. John Wiley & Sons, 2011. →pages 12, 18, 19, 23[22] R. Diao, V. Vittal, and N. Logic, “Design of a real-time security assessment tool forsituational awareness enhancement in modern power systems,” IEEE Transactionson Power Systems, vol. 25, no. 2, pp. 957–965, May 2010. → pages 12, 21, 23[23] M. He, J. Zhang, and V. Vittal, “A data mining framework for online dynamicsecurity assessment: Decision trees, boosting, and complexity analysis,” in 2012IEEE PES Innovative Smart Grid Technologies, Jan. 2012, pp. 1–8. → pages 12, 21[24] Y. Xu, Z. Y. Dong, J. H. Zhao, P. Zhang, and K. P. Wong, “A reliable intelligentsystem for real-time dynamic security assessment of power systems,” IEEE Trans-actions on Power Systems, vol. 27, no. 3, pp. 1253–1263, Aug. 2012. → pages 12,21, 23[25] I. Kamwa, S. Samantaray, and G. Joos, “Development of rule-based classifiers forrapid stability assessment of wide-area post-disturbance records,” IEEE Transac-tions on Power Systems, vol. 24, no. 1, pp. 258–270, Feb. 2009. → pages 12, 21, 23,24, 40[26] A. Kaci, I. Kamwa, L. Dessaint, S. Guillon et al., “Synchrophasor data baseliningand mining for online monitoring of dynamic security limits,” IEEE Transactionson Power Systems, vol. 29, no. 6, pp. 2681–2695, Nov 2014. → pages 12[27] J. Geeganage, U. Annakkage, T. Weekes, B. Archer et al., “Application of energy-based power system features for dynamic security assessment,” IEEE Transactionson Power Systems, vol. 30, no. 4, pp. 1957–1965, Jul. 2015. → pages 12[28] M. He, V. Vittal, and J. Zhang, “Online dynamic security assessment with miss-ing PMU measurements: A data mining approach,” IEEE Transactions on PowerSystems, vol. 28, no. 2, pp. 1969–1977, May 2013. → pages 12, 40[29] I. Kamwa, S. Samantaray, and G. Joos, “Catastrophe predictors from ensembledecision-tree learning of wide-area severity indices,” IEEE Transactions on SmartGrid, vol. 1, no. 2, pp. 144–158, Sept 2010. → pages 12, 21, 23, 3599[30] Y. Xu, Z. Y. Dong, K. Meng, R. Zhang, and K. P. Wong, “Real-time transientstability assessment model using extreme learning machine,” IET generation, trans-mission & distribution, vol. 5, no. 3, pp. 314–322, Mar. 2011. → pages 12, 21[31] L. Breiman, J. Friedman, C. J. Stone, and R. Olshen, Classification and RegressionTrees. Wadsworth, 1984. → pages 12, 27, 28, 93[32] J. H. Friedman, “Multivariate adaptive regression splines,” The annals of statistics,pp. 1–67, 1991. → pages 12, 27, 31, 32, 34, 93[33] M. He, J. Zhang, and V. Vittal, “Robust online dynamic security assessment usingadaptive ensemble decision-tree learning,” IEEE Transactions on Power Systems,vol. 28, no. 4, pp. 4089–4098, 2013. → pages 13, 40, 41[34] S. Samantaray, I. Kamwa, and G. Joos, “Ensemble decision trees for phasor mea-surement unit-based wide-area security assessment in the operations time frame,”IET Generation, Transmission & Distribution, vol. 4, no. 12, pp. 1334–1348, 2010.→ pages 13, 21, 23, 40[35] T. Takagi, Y. Yamakoshi, M. Yamaura, R. Kondow, and T. Matsushima, “Develop-ment of a new type fault locator using the one-terminal voltage and current data,”IEEE Transactions on Power apparatus and systems, no. 8, pp. 2892–2898, Aug.1982. → pages 13, 50[36] C. E. de Morais Pereira and L. C. Zanetta, “Fault location in transmission linesusing one-terminal postfault voltage data,” IEEE Transactions on Power Delivery,vol. 19, no. 2, pp. 570–575, Mar. 2004. → pages 13, 50[37] A. T. Johns and S. Jamali, “Accurate fault location technique for power transmissionlines,” IEE Proceedings C - Generation, Transmission and Distribution, vol. 137,no. 6, pp. 395–402, Nov. 1990. → pages 13, 50[38] J. Izykowski, E. Rosolowski, M. M. Saha, M. Fulczyk, and P. Balcerek, “A fault-location method for application with current differential relays of three-terminallines,” IEEE Transactions on Power Delivery, vol. 22, no. 4, pp. 2099–2107, Oct.2007. → pages 13, 14, 50[39] D. Novosel, D. G. Hart, E. Udren, and J. Garitty, “Unsynchronized two-terminalfault location estimation,” IEEE Transactions on Power Delivery, vol. 11, no. 1, pp.130–138, Jan. 1996. → pages 13, 50100[40] J. Izykowski, E. Rosolowski, P. Balcerek, M. Fulczyk, and M. M. Saha, “Accu-rate noniterative fault-location algorithm utilizing two-end unsynchronized mea-surements,” IEEE Transactions on Power Delivery, vol. 26, no. 2, pp. 547–555,Apr. 2011. → pages 13, 14, 50[41] D. Novosel, V. Madani, B. Bhargave, K. Vu, and J. Cole, “Dawn of the grid synchro-nization,” IEEE Power and Energy Magazine, vol. 6, no. 1, pp. 49–60, Jan. 2008.→ pages 13, 16, 50[42] A. S. Dobakhshari, “Noniterative parameter-free fault location on untransposedsingle-circuit transmission lines,” IEEE Transactions on Power Delivery, vol. 32,no. 3, pp. 1636–1644, Jun. 2017. → pages 14, 48, 50, 51, 52, 78[43] M. Kezunovic, “Smart fault location for smart grids,” IEEE Transactions on SmartGrid, vol. 2, no. 1, pp. 11–22, Mar. 2011. → pages 48[44] G. Feng and A. Abur, “Fault location using wide-area measurements and sparseestimation,” IEEE Transactions on Power Systems, vol. 31, no. 4, pp. 2938–2945,Jul. 2016. → pages 14, 50, 79, 80[45] T. Becejac, P. Dehghanian, and M. Kezunovic, “Impact of the errors in the PMU re-sponse on synchrophasor-based fault location algorithms,” in North American PowerSymposium (NAPS), 2016. IEEE, 2016, pp. 1–6. → pages 14, 78[46] M. Young and A. Silverstein, “Factors affection PMU installation cost,” AmericanRecovery and Reinvestment Act of 2009, 2014. → pages 14, 78[47] M. M. Saha, J. J. Izykowski, and E. Rosolowski, Fault location on power networks.Springer Science & Business Media, 2009. → pages 14, 47, 48, 49, 50, 79[48] G. Feng and A. Abur, “Identification of faults using sparse optimization,” in 201452nd Annual Allerton Conference on Communication, Control, and Computing(Allerton), Sept 2014, pp. 1040–1045. → pages 14, 79[49] Y. Liao, “Fault location for single-circuit line based on bus-impedance matrix uti-lizing voltage measurements,” IEEE Transactions on Power Delivery, vol. 23, no. 2,pp. 609–617, Apr. 2008. → pages 80[50] Q. Jiang, X. Li, B. Wang, and H. Wang, “PMU-based fault location using volt-age measurements in large transmission networks,” IEEE Transactions on PowerDelivery, vol. 27, no. 3, pp. 1644–1652, Jul. 2012. → pages 80101[51] S. Azizi and M. Sanaye-Pasand, “A straightforward method for wide-area fault lo-cation on transmission networks,” IEEE Transactions on Power Delivery, vol. 30,no. 1, pp. 264–272, Feb. 2015. → pages 80[52] A. Salehi-Dobakhshari and A. M. Ranjbar, “Application of synchronised phasor mea-surements to wide-area fault diagnosis and location,” IET Generation, TransmissionDistribution, vol. 8, no. 4, pp. 716–729, Apr. 2014. → pages 14, 79, 80[53] P. Kundur, J. Paserba, V. Ajjarapu, G. Andersson, A. Bose, C. Canizares,N. Hatziargyriou, D. Hill, A. Stankovic, C. Taylor, T. V. Cutsem, and V. Vittal,“Definition and classification of power system stability ieee/cigre joint task forceon stability terms and definitions,” IEEE Transactions on Power Systems, vol. 19,no. 3, pp. 1387–1401, Aug. 2004. → pages 17[54] S. Zadehkhost, “Efficient algorithms to expedite transient stability analysis of powersystems,” Ph.D. dissertation, The Univeristy of British Columbia, 2015. → pages18[55] A. Pai, Energy function analysis for power system stability. Springer Science &Business Media, 1989. → pages 18[56] S. Zhao, H. Jia, D. Fang, Y. Jiang, and X. Kong, “Criterion to evaluate power systemonline transient stability based on adjoint system energy function,” IET Generation,Transmission & Distribution, vol. 9, no. 1, pp. 104–112, 2014. → pages 20[57] I. Kamwa, S. Samantaray, and G. Joo´s, “On the accuracy versus transparency trade-off of data-mining models for fast-response pmu-based catastrophe predictors,” IEEETransactions on Smart Grid, vol. 3, no. 1, pp. 152–161, 2012. → pages 21, 22, 35,40, 45[58] W. P. Servcices. (2013, Dec.) WECC planing services. [Online]. Available:https://www.wecc.biz/PlanningServices/Pages/Default.aspx → pages 22, 34[59] M. Rahmatian, W. G. Dunford, and A. Moshref, “PMU based system protectionscheme,” in Proceeding of the 14th Electrical Power and Energy Conference (EPEC),12–14 Nov, 2014. → pages 22, 25[60] M. Rahmatian, W. G. Dunford, A. Palizban, and A. Moshref, “Transient stabilityassessment of power systems through wide-area monitoring system,” in Proceedingof 2015 IEEE PES General Meeting, 26–30 Jul., 2015. → pages 22[61] A. N. Michel, A. Fouad, and V. Vittal, “Power system transient stability usingindividual machine energy functions,” IEEE Transactions on Circuits and Systems,vol. 30, no. 5, pp. 266–276, 1983. → pages 23102[62] A. Fouad and S. Stanton, “Transient stability of a multi-machine power system parti: investigation of system trajectories,” IEEE Transactions on Power Apparatus andSystems, no. 7, pp. 3408–3416, 1981. → pages 23[63] J. C. Cepeda, J. L. Rueda, D. G. Colome´, and D. E. Echeverr´ıa, “Real-time transientstability assessment based on centre-of-inertia estimation from phasor measurementunit records,” IET Generation, Transmission & Distribution, vol. 8, no. 8, pp. 1363–1376, 2014. → pages 23, 36[64] North American SynchroPhasor Initiative (NASPI) RITT Report. (2011, Jun.)Guidelines for siting phasor measurement units. [Online]. Available:https://www.naspi.org/documents → pages 23, 35[65] P. W. Sauer, Power system dynamics and stability. Prentice Hall, 1998. → pages24, 25, 26[66] C. Fu and A. Bose, “Contingency ranking based on severity indices in dynamicsecurity analysis,” IEEE Transactions on Power Systems, vol. 14, no. 3, pp. 980–985, Aug. 1999. → pages 24[67] S. M. Rovnyak, “Integral square generator angle index for stability assessment,” inProceedings of the IEEE Power Engineering Society Winter Meeting, 28 Jan.–1 Feb.,2001. → pages 25[68] I. Kamwa, A. K. Pradhan, G. Joos, and S. Samantaray, “Fuzzy partitioning of a realpower system for dynamic vulnerability assessment,” IEEE Transactions on PowerSystems, vol. 24, no. 3, pp. 1356–1365, Aug. 2009. → pages 25[69] L. Wehenkel, V. Cutsem, and M. Ribbens-Pavella, “An artificial intelligence frame-work for online transient stability assessment of power systems,” IEEE Transactionson Power Systems, vol. 4, no. 2, pp. 789–800, 1989. → pages 26[70] L. Wehenkel and M. Pavella, “Decision trees and transient stability of electric powersystems,” Automatica, vol. 27, no. 1, pp. 115–134, 1991. → pages 26[71] T. Van Cutsem, L. Wehenkel, M. Pavella, B. Heilbronn, and M. Goubin, “Decisiontree approaches to voltage security assessment,” IEE Proceedings C (Generation,Transmission and Distribution), vol. 140, no. 3, pp. 189–198, 1993. → pages 26[72] N. Senroy, G. T. Heydt, and V. Vittal, “Decision tree assisted controlled islanding,”IEEE Transactions on Power Systems, vol. 21, no. 4, pp. 1790–1797, 2006. → pages26, 95103[73] T. Hastie, R. Tibshirani, and J. Friedman, The elements of statistical learning.Springer, 2009, vol. 2, no. 1. → pages 28, 29, 30, 31, 32, 33, 34[74] DSATools. (2014, January) TSAT, Transient Stability Assessment Tool. [Online].Available: http://www.dsatools.com → pages 31, 35, 85, 86[75] N. Amjady and S. Banihashemi, “Transient stability prediction of power systems bya new synchronism status index and hybrid classifier,” IET generation, transmission& distribution, vol. 4, no. 4, pp. 509–518, 2010. → pages 35[76] A. D. Rajapakse, F. Gomez, K. Nanayakkara, P. A. Crossley, and V. V. Terzija, “Ro-tor angle instability prediction using post-disturbance voltage trajectories,” IEEETransactions on Power Systems, vol. 25, no. 2, pp. 947–956, 2010. → pages 35[77] J. Cepeda, D. Colome´, and N. Castrillo´n, “Dynamic vulnerability assessment due totransient instability based on data mining analysis for smart grid applications,” inProceedings of the IEEE PES Conference on Innovative Smart Grid Technologies,19–21 Oct., 2011. → pages 36[78] SPM. (2015, February) CART, Classification And Regression Tree. [Online].Available: http://www.salford-systems.com → pages 36[79] ——. (2015, February) MARS, Multivariate Adaptive Regression Splines. [Online].Available: http://www.salford-systems.com → pages 36[80] M. Rahmatian, Y. C. Chen, A. Palizban, A. Moshref, and W. G. Dunford, “Transientstability assessment via decision trees and multivariate adaptive regression splines,”Electric Power Systems Research, vol. 142, pp. 320–328, 2017. → pages 36, 57[81] K. Sun, S. Likhate, V. Vittal, V. S. Kolluri, and S. Mandal, “An online dynamicsecurity assessment scheme using phasor measurements and decision trees,” IEEETransactions on Power Systems, vol. 22, no. 4, pp. 1935–1943, 2007. → pages 40[82] A. Riepnieks and H. Kirkham, “An introduction to goodness of fit for PMU parame-ter estimation,” IEEE Transactions on Power Delivery, vol. 32, no. 5, pp. 2238–2245,Oct. 2017. → pages 48, 50, 51, 59, 61, 93[83] A. J. Roscoe, I. F. Abdulhadi, and G. M. Burt, “P and M class phasor measure-ment unit algorithms using adaptive cascaded filters,” IEEE Transactions on PowerDelivery, vol. 28, no. 3, pp. 1447–1459, Jul. 2013. → pages 51, 59[84] P. Castello, J. Liu, C. Muscas, P. A. Pegoraro, F. Ponci, and A. Monti, “A fastand accurate PMU algorithm for P+M class measurement of synchrophasor and fre-104quency,” IEEE Transactions on Instrumentation and Measurement, vol. 63, no. 12,pp. 2837–2845, Dec 2014. → pages 51[85] I. Ivanov and A. Murzin, “Synchrophasor-based transmission line parameter es-timation algorithm taking into account measurement errors,” in 2016 IEEE PESInnovative Smart Grid Technologies Conference Europe, Oct. 2016. → pages 52[86] J. J. Grainger and W. D. Stevenson, Power system analysis. McGraw-Hill, 1994.→ pages 53[87] H. Kirkham and A. Riepnieks, “Dealing with non-stationary signals: Definitions,considerations and practical implications,” in Proc. Power and Energy Society Gen-eral Meeting, Aug. 2016. → pages 59[88] B. M. Schettino, C. A. Duque, P. M. Silveira, P. F. Ribeiro, and A. S. Cerqueira, “Anew method of current-transformer saturation detection in the presence of noise,”IEEE Transactions on Power Delivery, vol. 29, no. 4, pp. 1760–1767, Aug. 2014. →pages 66[89] E. Hajipour, M. Vakilian, and M. Sanaye-Pasand, “Current-transformer saturationprevention using a controlled voltage-source compensator,” IEEE Transactions onPower Delivery, vol. 32, no. 2, pp. 1039–1048, Apr. 2017. → pages 66[90] M. Ramanathan, C. Chen, W. G. Dunford, and F. Rahmatian, “Synchrophasor dataqualification measures for reliable fault location,” in Protection, Automation, andControl (PAC) World Conference, Aug. 2017. → pages 76[91] A. S. Dobakhshari and A. M. Ranjbar, “A novel method for fault location of trans-mission lines by wide-area voltage measurements considering measurement errors,”IEEE Transactions on Smart Grid, vol. 6, no. 2, pp. 874–884, Mar. 2015. → pages79[92] J. Chang, G. N. Taranto, and J. H. Chow, “Dynamic state estimation in powersystem using a gain-scheduled nonlinear observer,” in Proceedings of the 4th IEEEConference on Control Applications, 1995. → pages 80[93] E. Scholtz, P. Sonthikorn, G. Verghese, and B. Lesieutre, “Observers for swing-stateestimation of power systems,” in North American Power Symposium, 2002. → pages80[94] I. Stolz, “Observers and graphic displays for the swing motions of a power system,”Ph.D. dissertation, Masters thesis, Technical University of Berlin (research done atMIT), 1994. → pages 81105[95] H. Modir and R. A. Schlueter, “A dynamic state estimator for dynamic securityassessment,” IEEE Transactions on Power Apparatus and Systems, no. 11, pp. 4644–4652, 1981. → pages 81[96] “IEEE recommended practice for excitation system models for power system sta-bility studies - redline,” IEEE Std 421.5-2016 (Revision of IEEE Std 421.5-2005) -Redline, pp. 1–453, Aug. 2016. → pages 86[97] S. Rovnyak, S. Kretsinger, J. Thorp, and D. Brown, “Decision trees for real-timetransient stability prediction,” IEEE Transactions on Power Systems, vol. 9, no. 3,pp. 1417–1426, 1994. → pages 95[98] K. Mei and S. M. Rovnyak, “Response-based decision trees to trigger one-shot sta-bilizing control,” IEEE Transactions on Power Systems, vol. 19, no. 1, pp. 531–537,2004. → pages 95[99] F. Dorfler and F. Bullo, “Kron reduction of graphs with applications to electricalnetworks,” IEEE Transactions on Circuits and Systems I: Regular Papers, vol. 60,no. 1, pp. 150–163, Jan. 2013. → pages 95[100] R. Ramanathan and B. Tuck, “BPA’s experience of implementing node breakermodel for power system operations studies,” in 48th International Universities PowerEngineering Conference (UPEC), Sep. 2013. → pages 96[101] M. Heidarifar and H. Ghasemi, “A network topology optimization model basedon substation and node-breaker modeling,” IEEE Transactions on Power Systems,vol. 31, no. 1, pp. 247–255, Jan. 2016. → pages 96106Appendix AParameters of Applied ModelsThe specs of the exciter, governor and power systems stabilizer applied in Chapter 4 areshown Tables A.1–A.3.Table A.1: Parameters of the exciter model ESAC4A.Parameter ValueTR 0.0000VI,max 0.5000VI,min −0.5000TC 0.0000TB 0.0000KA 100.0TA 0.0200VR,max 7.0000VR,min −7.0000KC 0.0000Table A.2: Parameters of the governor model GOV6.Parameter ValueR 50.00Vmax 1.000Vmin 0.000T1 0.500T2 3.000T3 9.000Pmax 1000107Table A.3: Parameters of the power system stabilizer model PSS1.Parameter ValueT1 0.08T2 0.08T3 0.04T4 0.04T5 2.00K5 6.00Hlim 0.10108