A Tool for On-line Stability Determination and Control for Coordinated Operations between Regional E

A Tool for On-line Stability Determination and Control for Coordinated Operations between Regional E

Info iconThis preview shows page 1. Sign up to view the full content.

View Full Document Right Arrow Icon
This is the end of the preview. Sign up to access the rest of the document.

Unformatted text preview: PSERC A Tool for On-line Stability Determination and Control for Coordinated Operations between Regional Entities Using PMUs Final Project Report Power Systems Engineering Research Center A National Science Foundation Industry/University Cooperative Research Center since 1996 Power Systems Engineering Research Center A Tool for On-line Stability Determination and Control for Coordinated Operations between Regional Entities Using PMUs Final Project Report Project Team Vijay Vittal, Project Leader, Arizona State University Gerald T. Heydt, Arizona State University A.P. Sakis Meliopoulos, Georgia Institute of Technology PSERC Publication 08-01 January 2008 Information about this project For information about this project contact: Vijay Vittal Ira A. Fulton Chair Professor Department of Electrical Engineering Arizona State University PO Box 875706 Tempe, AZ 85287-5706 Tel: 480-965-1879 Fax: 480-965-0745 Email: [email protected] Power Systems Engineering Research Center This is a project report from the Power Systems Engineering Research Center (PSERC). PSERC is a multi-university Center conducting research on challenges facing a restructuring electric power industry and educating the next generation of power engineers. More information about PSERC can be found at the Center’s website: http://www.pserc.org. For additional information, contact: Power Systems Engineering Research Center Arizona State University 577 Engineering Research Center Box 878606 Tempe, AZ 85287-8606 Phone: 480-965-1643 FAX: 480-965-0745 Notice Concerning Copyright Material PSERC members are given permission to copy without fee all or part of this publication for internal use if appropriate attribution is given to this document as the source material. This report is available for downloading from the PSERC website. © 2007 Arizona State University and Georgia Institute of Technology. All rights reserved. Acknowledgements This is the final report for the Power Systems Engineering Research Center (PSERC) research project S-27 titled “A Tool for On-line Stability Determination and Control for Coordinated Operating between Regional Entities Using PMUs.” We express our appreciation for the support provided by PSERC’s industrial members, and by the National Science Foundation under grants NSF EEC-0001880 at Arizona State University and NSF EEC-0080012 at Georgia Institute of Technology received under the Industry/ University Cooperative Research Center program. The authors wish to recognize their postdoctoral researchers and graduate students that contributed to the research and creation of the reports: Dr. Kai Sun (Arizona State University), Siddharth Likhate (Arizona State University), George Stefopoulos, (graduate research assistant at Georgia Institute of Technology), and Dr. Salman Mohagheghi, (Postdoctoral at Georgia Institute of Technology). The authors thank all PSERC members for their technical advice on the project, especially the industry advisors for the project: • • • • • Sharma Kolluri – Entergy Corporation Sujit Mandal – Entergy Corporation Navin Bhatt – American Electric Power Co. Lisa Beard – Tennessee Valley Authority Kip Morison – Powertech Labs. . i Executive Summary This project examines the application of phasor measurement units (PMUs) for on-line transient stability assessment in an operational decision-making environment. The project involved the participation of three PSERC member companies, American Electric Power, Entergy and Tennessee Valley Authority, who are active participants in the Eastern Interconnection Phasor Project (EIPP) sponsored by the U.S. Department of Energy’s Consortium for Electric Reliability Technology Solutions. The primary motivation for the PSERC PMU project arises from the critical need to develop analytical tools to enhance power system reliability, from the goal to provide superior operational situation awareness using accurate real-time measurements, and from the search for innovative applications for measurements made using PMUs. A survey of the three companies was conducted to determine key features of their use of PMUs. The goal of the survey was to determine the type of data that can be obtained from the PMU measurements, the frequency with which the data is obtained, and the feasibility of using this data for real-time transient stability assessment. The survey established the feasibility and advantages of using the PMU data in real-time transient stability assessment. The results of the surveys, and a literature review, also showed that nomograms are the main method used for insuring system stability during system operations (e.g., line switching, generation level changes). A nomogram is a tool is created using off-line analyses, but used in real-time decision-making by system operators. The nomogram approach provides an approximation of actual transmission loading based on power flow and stability studies done off-line. The results are used in real time by matching the operating conditions to the “stored” off-line simulations. With a nomogram, several key system parameters are plotted on a graph-like sheet. By joining points on the graph, the transmission loading is estimated. While the word nomogram implies a physical graph (e.g., drawn on a sheet of paper), the present practice is to automate the process via a table look-up. In effect, the given system conditions for the stability assessment (i.e., present line loadings, system loads and key generation levels) are entered and maximum individual circuit loads are then read in the table look-up. The main thrust of this project was to develop methodologies that use real-time data to predict instabilities for the present operating conditions and to guide actions that can prevent in the instabilities from occurring under a given system operating conditions. Since the nomogram approach does not capture the value of accurate real-time PMU data, different approaches were developed for on-line transient stability assessment using PMU measurements. This approach requires substantial computational power in real time as compared to the use of nomograms, but provides the advantage of providing more accurate information needed for operating a power system reliably, securely and efficiently. ii Volume I Volume I of this report presents an approach to on-line dynamic security assessment using PMU measurements in conjunction with off-line transient stability analysis on a near real-time operational model of the power system. The technique is developed using the operational model of the Entergy power system with 2100 buses, 2600 transmission lines, and 240 generators. Based on the load forecast and information regarding component availability, Entergy has the capability to create an accurate representation of the on-line system for a 24-hour horizon. This system representation includes the network power flow database, and the associated dynamic data and modeling specification needed to conduct detailed time-domain (T-D) simulations. The proposed on-line transient stability assessment approach is developed in three stages: i. For the 24-hour horizon, on-line system data, power flow analysis is used to model a series of operating conditions representing the projected variation in daily load with the unit-commitment-based generation pattern. Exhaustive T-D simulations for “n-1” contingencies and probable “n-k” contingencies are conducted on the generated cases and stored in a database. ii. A decision tree is then trained using the database obtained. The decision tree is used to identify critical attributes from system parameters that characterize the system dynamic performance and to determine the critical attribute thresholds that result in system insecurity. iii. These critical attributes serve as the measurements to be made using PMUs. For contingencies analyzed, the real-time measurements are compared to their thresholds stored in the decision tree to determine related paths and terminal nodes. An insecurity score is calculated for each path. If any score exceeds a preset limit and the associated contingencies have high probabilities of occurrence, then appropriate preventive control actions can be determined and prepared as necessary. The results obtained using the approach developed in Volume I show that the proposed approach can build decision trees with high prediction accuracy to reliably identify stability problems for prospective operating conditions and provide guidelines for preventive control. The approach also identifies critical security indicators which could indicate candidates for new PMU locations. A new path-based decision tree classification method is also proposed and compared with the traditional terminal nodes-based method to regulate the decision tree’s prediction reliability for changes in operating conditions. This new method focuses on the development of self-adaptive rules to make better use of available decision trees and introduces a shift from “Black-Box” decision trees to “White-Box” decision trees, or rule-based decision tree classification systems. The decision trees built in this project can reach prediction accuracies of about 97% for insecure cases. This research is a first step towards the goal of more practical versions of the proposed scheme that can achieve high accuracy (e.g. 99% or even 99.9%). To realize that goal, several issues will have to be carefully analyzed. One aspect that we are currently investigating is the impact of the use of a larger number of operating conditions to increase the number of insecure cases used in decision tree training. iii Volume II Volume II describes an approach to determining the dynamic real-time model of the system based on a dynamic state estimator that uses PMU measurements in addition to other available measurements. The methodology is based on energy methods. Specifically, the dynamic real time model is used to compute the energy function of the system (Lyapunov function) and subsequently the stability margin of the system and whether instability is eminent. This approach is complimentary to the approaches described in Volume I. A real-time dynamic model is used to determine the stability margin of the system and, therefore, how close the system is to instability of the present operating conditions as well as a number of plausible contingencies. To achieve a reliable model requires accurate information about current operating conditions and about measurement errors from the instrumentation channels normally used in substations. Conceptually, PMU data provide the needed accuracy because they are time-tagged with accuracy better than 1 microsecond and magnitude accuracy that is better than 0.1%. In addition, PMU data include frequency and rate of change of frequency. An analysis of the observability of the real-time dynamic model was conducted using a specialized application of the SuperCalibrator, a dynamic state estimator that operates locally in a substation using PMU data, and other data available from relays, RTUs, etc. The analysis showed that the dynamic state of the system can be obtained from standard substation data and GPS synchronized devices. The conclusion is that present day PMUs have the accuracy need to extract the real-time dynamic model. A real-time dynamic model was built for a case study with a test system with two generating units and interconnecting transmissions lines. Several scenarios of stable and unstable swings were assessed. The performance of the proposed method in predicting the outof-step condition showed that out-of-step conditions can be assessed much earlier than conventional methods (such as standard out-of-step relaying). The results obtained in this project indicate that, with further development, the tools developed can be used in an operational EMS setting for evaluating stability limits and developing preventive and control options to steer the system away from vulnerable situations. Research grade software implements the entire approach based on data provided by the DSATools suite of software developed by Powertech Labs. The approach also uses commercial grade decision tree software called CART developed by Salford Systems. The developed technology in this project is commercializable with more development work. Specifically, the proposed method can lead to a new predictive on-line dynamic security assessment scheme. iv A Tool for On-line Stability Determination and Control for Coordinated Operations between Regional Entities Using PMUs Volume 1 ` Project Team Vijay Vittal, Project Leader, Arizona State University Gerald T. Heydt, Arizona State University A.P. Sakis Meliopoulos, Georgia Institute of Technology Information about this volume For information about this volume contact: Vijay Vittal Ira A. Fulton Chair Professor Department of Electrical Engineering Arizona State University PO Box 875706 Tempe, AZ 85287-5706 Tel: 480-965-1879 Fax: 480-965-0745 Email: [email protected] Power Systems Engineering Research Center This is a project report from the Power Systems Engineering Research Center (PSERC). PSERC is a multi-university Center conducting research on challenges facing a restructuring electric power industry and educating the next generation of power engineers. More information about PSERC can be found at the Center’s website: http://www.pserc.org. For additional information, contact: Power Systems Engineering Research Center Arizona State University 577 Engineering Research Center Box 878606 Tempe, AZ 85287-8606 Phone: 480-965-1643 FAX: 480-965-0745 Notice concerning copyright material PSERC members are given permission to copy without fee all or part of this publication for internal use if appropriate attribution is given to this document as the source material. This report is available for downloading from the PSERC website. © 2007 Arizona State University and Georgia Institute of Technology. All rights reserved. Table of Contents 1. Introduction................................................................................................................. 1 1.1 Background......................................................................................................... 1 1.2 Off-line Approaches ........................................................................................... 2 1.3 Phasor Measurement Units: a Review............................................................. 2 1.4 Phasor Measurement Units: The IEEE Standard for Synchrophasors for Power Systems .................................................................................................... 5 1.5 The Migration from Wide Area Measurements to Wide Area Control........ 7 1.6 Voltage Stability Limits ..................................................................................... 8 1.7 Angular Stability Limits .................................................................................. 10 1.8 Additional Considerations ............................................................................... 11 1.9 On-line Transient Stability Assessment ......................................................... 13 1.10 Overview of Analysis Techniques ................................................................... 13 1.10.1 Traditional time domain simulation ............................................................ 14 1.10.2 Direct methods ............................................................................................ 14 1.10.3 Artificial learning ........................................................................................ 15 1.11 Decision Trees for Transient Stability Assessment ....................................... 16 1.12 Report Organization ........................................................................................ 17 2. Theoretical Background........................................................................................... 19 2.1 Definition........................................................................................................... 19 2.2 CART Methodology ......................................................................................... 19 2.2.1 Splitting criterion ........................................................................................... 19 2.2.2 Stop-splitting criterion and DT pruning...................................................... 21 2.3 Decision Tree Building (preliminary tasks) ................................................... 22 2.3.1 Generation of operating scenarios .............................................................. 22 2.3.2 Stability assessment .................................................................................... 22 2.3.3 Selection of attributes .................................................................................. 22 2.3.4 Preparation of database .............................................................................. 23 2.3.5 Multi-class and multi-contingency DT considerations ............................... 23 2.4 Use of Phasor Measurement Units for Transient Stability Assessment ...... 24 3. Proposed Scheme ...................................................................................................... 25 3.1 DT Based Online Security Assessment Scheme............................................. 25 3.1.1 Offline DT building ..................................................................................... 26 3.1.2 Periodic OC prediction and DT updating ................................................... 27 3.1.3 Online security assessment.......................................................................... 27 3.2 Reliability of the DT ......................................................................................... 28 3.3 Preventive Control ........................................................................................... 31 3.4 Building a Better DT by Adjusting Penalties for Predictors........................ 32 4. Entergy Case Study................................................................................................... 34 4.1 System Description ........................................................................................... 34 4.2 Database Generation........................................................................................ 35 4.2.1 Operating conditions ................................................................................... 35 4.2.2 Contingencies .............................................................................................. 35 4.2.3 Candidate attributes .................................................................................... 36 Table of Contents (continued) 4.3 Results................................................................................................................ 37 4.3.1 Decision trees for transient stability (DT-U) ................................................. 38 4.3.2 Decision trees for damping (DT-D) ............................................................ 41 4.4 Evaluation of DT reliability............................................................................. 47 4.5 Conclusions ....................................................................................................... 50 References........................................................................................................................ 51 Project Publications ........................................................................................................ 55 ii List of Figures Fig. 1.1 Conceptual diagram of a synchronized phasor measuring system redrawn from [4] ............................................................................................................. 3 Fig. 1.2 Pictorial of the transition of WAMS to WACS..................................................... 8 Fig 1.3 Classification of on-line transient stability assessment methods taken directly from [34] ......................................................................................................... 14 Fig 1.4 Illustration of the Transient Energy Function (TEF) method............................... 15 Fig. 2.1 Example of a decision tree................................................................................... 20 Fig. 3.1 Flowchart describing the proposed scheme......................................................... 25 Fig. 3.2 Sample DT illustrating new path based approach (Make all Italics in figure) .... 29 Fig.4.1 Entergy transmission system map ........................................................................ 34 Fig.4.2 Amite South area (showing PMU locations)........................................................ 36 Fig. 4.3 Details of predictors used at each stage of the DT-U1 ........................................ 39 Fig. 4.4 Details of predictors used at each stage of DT-U2.............................................. 40 Fig. 4.5 Details of predictors used at each stage of DT-D1.............................................. 42 Fig. 4.6 Details of predictors used at each stage of DT-D2.............................................. 43 Fig. 4.7 Key transmission path ......................................................................................... 44 Fig. 4.8 Operational nomogram built using DT-D2 ......................................................... 46 iii List of Tables Table 1.1 Main contents of IEEE C37.118 ......................................................................... 6 Table 1.2 Influence quantities and allowable error limits for compliance levels 0-1......... 7 Table 1.3 Status and levels monitored for the calculation of voltage stability limit calculations .......................................................................................................... 9 Table 2.1 Sample database................................................................................................ 23 Table 3.1 Insecurity scores of paths for sample tree......................................................... 29 Table 3.2 Prediction accuracy for traditional method....................................................... 31 Table 3.3 Prediction accuracy for path based method ...................................................... 31 Table 4.1 Prediction accuracy for Learning Set (DT-U1) ................................................ 38 Table 4.2 Prediction accuracy for Test Set (DT-U1)........................................................ 38 Table 4.3 Prediction accuracy for Learning Set (DT-U2) ................................................ 38 Table 4.4 Prediction accuracy for Test Set (DT-U2)........................................................ 38 Table 4.5 Prediction accuracy for Learning Set (DT-D1) ................................................ 41 Table 4.6 Prediction accuracy for Test Set (DT-D1)........................................................ 41 Table 4.7 Prediction accuracy for Learning Set (DT-D2) ................................................ 41 Table 4.8 Prediction accuracy for Learning Set (DT-D2) ................................................ 41 Table 4.9 Prediction accuracy for July 26 cases (DT-D2)................................................ 47 Table 4.10 Insecurity scores of paths for DT-D2 ............................................................. 48 Table 4.11 Terminal node details for DT-D2 ................................................................... 49 Table 4.12 Prediction accuracy for July 26 cases using path based method (DT-D2) ..... 49 iv NOMENCLATURE A_i_j Phasor angle difference between buses i and j c(i | j ) FB I ( n) N ts t N ijs Cost of misclassifying a class j case as class i Fault bus Impurity score for node n Total number of cases in the test set Total number of class j cases misclassified as class i Probability that a case in node n is secure Probability that a case in node n is insecure Percentage of insecure cases in node i Power flow of line i-j (‘H’ indicates a 345 kV line) MW output of generator i Misclassification cost Insecurity score Weight of node i p ( S | n) p ( I | n) pi P_i_j PG_i R ts S ωi v 1. 1.1 Introduction Background The primary objective of this project was to develop applications using synchronized phasor data measurements that are being implemented at several electric utility companies in the US. The application that was specifically addressed in this project dealt with on-line transient stability assessment using the measurements from phasor measurement units (PMUs). In order to develop a clearer understanding of the capabilities of the PMUs installed in the field and the nature and quality of the measurements obtained by these units a questionnaire was developed and distributed to the three industry advisors on this project. The questionnaire is shown in the Appendix. The analysis of the responses to the questionnaire showed that at present many companies in the US have only a limited number of PMUs installed at specific locations in the bulk power network. In addition no specific criteria were used to actually choose the locations for placements of these PMUs and as a result the existing PMU locations may not provide the best measurements to characterize the transient stability behavior of the system. No PMUs currently exist at generator buses hence it is not possible to make a direct measurement of stability behavior primarily associated with generator response to disturbances. Having obtained this valuable insight from the results of the questionnaire, it was decided that the best option for using the existing PMU measurements for transient stability analysis consists of using an operational model of the system under consideration as close to real time as possible. This would then account for changes in system configuration and component outages for maintenance. Off-line analysis would then be conducted for a large number of probable contingencies over a range of operating conditions on this operational system model using conventionally transient stability analysis tools and transient stability limits would be obtained in terms of pre-disturbance system attributes corresponding to the measurements made by the existing PMUs. Since these attributes would then define a nomogram in the space of critical attributes, the synchronized PMU measurements could be monitored at defined time intervals to see if the operating condition would cross the boundary of an evaluated nomogram. If the boundary was approached then appropriated preventing or corrective control could be armed in anticipation of the contingency which would cause a violation of the threshold defined by the nomogram. If the contingency did not occur then the operation sequence would continue with any alteration to the normal operations procedures. If the contingency does occur then the armed control actions would be implemented to steer the system away from instability. 1 Given this setting the industry advisors on this project were then visited by each of the three PIs and a careful understanding of the existing off-line transient stability analysis policies and procedures was obtained. A summary of the existing practices is provided in what follows. A brief overview of the basis behind the synchrophasor measurements is also provided. 1.2 Off-line Approaches In this section several traditional concepts are revisited relating to stability assessment. For purposes of this section, off-line methods are discussed. The reason for discussing off-line methods is so that a base line can be obtained for assessing on-line stability assessment methods. Main inputs for this chapter were obtained from visits to electric utility companies that use nomogram type approaches to stability assessment. Stability, as viewed in this section, is divided into two areas: voltage stability and angular stability. In both cases, the problem is viewed as the calculation of operating limits. Because phasor measurement units are considered central in stability assessment (either on-line or off-line), a section of this chapter is relegated to this subject. 1.3 Phasor Measurement Units: a Review Phasor measurement units (PMUs) are instruments that take measurements of voltages and currents and time-stamp these measurements with high precision. PMUs are equipped with Global Positioning System receivers. The GPS receivers allow for the synchronization of the several readings taken at distant points [1]. To accomplish synchronization of measurements taken at distant points, several measurements are taken, and the measurements are time stamped. Interpolation is used to obtain estimates of measurements at a given time within the time horizon of the measurements. PMUs were developed from the invention of the symmetrical component distance relay (SCDR). The SCDR development outcome was a recursive algorithm for calculating symmetrical components of voltage and current [2]. Synchronization is made possible with the advent of the GPS satellite system [3]. The GPS system is a system of 36 satellites (of which 24 are used at one time) to produce time signals at the Earth’s surface. GPS receivers can resolve these signals into (x, y, z, t) coordinates. The t coordinate is time. This is accomplished by solving the distance= (rate) (time) in three dimensions using satellite signals. The PMU records the sequence currents and voltages and time stamps the reading with time obtained by the GPS receiver. It is possible to achieve accuracy of synchronization of 1 microsecond or 0.021° for 60 hertz signal. This is well in the suitable range of measuring power frequency voltages and currents [2]. The basic distance-rate-time formulation of this problem is “solved” using state estimation. That is, a least squares problem 2 formulation is used to find {x, y, z, t} which makes the distance-rate-time equations agree in the least squares sense. A minimum of four satellite readings are needed to obtain an observable problem to calculate {x, y, z, t}. Most PMUs have the ability to receive signals from at least ten satellites (i.e., ten channel receivers). Based upon the research done at Virginia Tech, the Macrodyne Company was able to begin production of PMU devices which led to the IEEE Standard 1344 “Sychrophasor” which defines the output data format of a PMU [2]. Figure 1.1 is a pictorial of PMU measurement system from [4]. Figure 1.1 Conceptual diagram of a synchronized phasor measuring system redrawn from [4] The calculation of the phasor measurement can be done using discrete Fourier transforms. A sinusoidal quantity representing voltage v(t ) = Vm cos(ωt + δ ) has a phasor representation 3 V jδ V= me . 2 By sampling v(t) every τ seconds, the total duration of the ensemble of samples is T=kτ. The discrete Fourier transform of v(t) is V (kΔΩ ) = 12 2N ⎛⎛ ⎜⎜ ⎜ ⎝⎝ ∑ i =1 N ⎞⎛ v(iτ ) cos(ikτΔΩ) ⎟ − j ⎜ ⎠⎝ ∑ v(iτ )sin(ikτΔΩ) ⎟ ⎟ ⎟ ⎠⎠ i =1 N ⎞⎞ where ΔΩ=2π/Nτ, and N is the number of samples in one period of the nominal power system frequency [5]. PMUs measure voltage and current with high accuracy at a rate of 2.88 kHz. It can calculate watts, vars, frequency, and phase angle 12 times per 60 hertz cycle. The actual sampling rate used to achieve this output is 1.4 MHz [6]. Some examples of uses of PMUs are fault recording, disturbance recording, and transmission and generation modeling [6]. PMUs are able to measure what was once virtually immeasurable: phase difference at different substations. With PMUs, the utilities are able to measure voltage phase angle The integration of PMUs into state estimation has been discussed in the literature [7-16]. There is a school of thought that the measurements from the PMU are far superior than the SCADA data used in traditional state estimation and should be collected and used separate from this data [17]. Others admit there is difference in the information and it is viable to use PMU measurements in with SCADA data [18]. Hydro-Quebec believes that the PMUs are accurate enough to not need correlation between PMU measurements. Their algorithm is to place the PMUs based on the buses which minimize the correlation between measurements between proposed PMU locations and present measurement locations [19]. An dramatic improvement in the state estimate has been seen by using a threephase model and the use of GPS synchronized measurements [20]. With increased need for multi-area of state estimation, there has been noted the possibility of the increased error in the state estimate as the size of systems grows [21]. PMUs are being investigated as a solution to this problem. The electric utilities in the regional transmission organization (RTO) would still do their state estimation. The RTO receives the results from the various state estimates of the areas under its control and PMU measurements from boundaries between electric utilities. The individual state estimators do not interact or exchange data with other state estimators. This allows for each estimator to have its own unique algorithm and not effect performance of other area estimators [12]. The Eastern Interconnect Phasor Project (EIPP) was created to encourage the use of advanced metering technology in the eastern interconnect. The primary technology focused on is PMUs. The Eastern Interconnect Phasor Project (EIPP) is divided into two stages. The near term goal is to use the expertise and equipment developed with sponsor- 4 ship from the U.S. Department of Energy to deliver immediate value to project participants in the eastern interconnect. Most of the existing expertise involves off-line analysis and is supportive of planning activities. The long-term goal is to add value to the interregional information system and measurement system using PMUs [22]. 1.4 Phasor Measurement Units: The IEEE Standard for Synchrophasors for Power Systems The IEEE has recognized the need for standards for PMUs, also known as synchrophasors. The first standard for PMUs, IEEE 1344 [23], was written in 1995. The drafting of a standard for PMUs is, perhaps, documentation that PMUs are expected to occupy a significant role in power systems instrumentation. A working group was formed in January 2001 to create a new standard for PMUs, IEEE C37.118 [24]. The updated standard provides clarification for phasor and synchronized phasor definitions. The standard defines synchronized phasor measurements in substations so that the measurement equipment can be readily interfaced with associated systems. Table 1.1 lists the major contents of the updated standard. To allow for integration of PMUs with other equipment, the standard provides common data format for exchanging information with PMUs. The data for time measurement shall consist of second-of-century (SOC) counts, fraction of second count, and a time status value. The SOC count is the number of seconds since the calendar time from midnight January 1, 1970. The accuracy of the time stamp required is ±1 μs. The maximum phase time error allowable is 26 microseconds. Table 1.2 shows all the limitations imposed on the PMU. 5 Table 1.1 Main contents of IEEE C37.118 (taken directly from [24]) Body of Standard Synchrophasor measurement • o Definition of phasor and synchrophasor o Measurement time tag • o System time synchronization Synchrophasor measurement requirements and • compliance verification • o Synchrophasor estimation o Accuracy limits • o Compliance verification Message format • o Message application • o Message framework o Data frame • o Configuration frame o Header frame o Command frame Appendices Cyclic redundancy check codes Time tagging and transient response Message examples Sources of synchronization Time and synchronization communication Benchmark tests TVE evaluation and PMU testing Synchrophasor message mapping into communications • • • 6 Table 1.2 Influence quantities and allowable error limits for compliance levels 0-1 (taken directly from [24]) Range of influence quantity change with respect to reference and maximum allowable TVE in percent Influence Reference (%) for each compliance level quantity condition Level 0 Level 1 Range TVE (%) Range TVE (%) Signal frequency fnominal ± 0.5 Hz 1 ± 5 Hz 1 Signal magnitude 100 % 80 – 120% 1 10 – 120% 1 rated rated rated Phase angle 0 radian ± π radians 1 ± π radians 1 Harmonic distor< 0.2 % 1% any 1 10% any 1 tion THD harmonic harmonic up up to 50th to 50th Our band of in< 0.2 of 1% of input 1 10% of in1 terfering signal, input sigsignal magput signal at frequency fi, nal magninitude magnitude where tude |fi – f0|>fs/2, fs= phasor reporting rate, f0=fnominal 1.5 The Migration from Wide Area Measurements to Wide Area Control Wide area measurement systems (WAMS) are instrumentation infrastructures that span a wide geographic area (e.g. typically several control areas and potentially several operating companies). Under WAMS, the time required to transmit the sensory information (latency) back to the central control center is significant compared to the dynamics of the measurement. The latency for a measurement sent via satellite can be as high as 250 ms, and up to 100 ms for a signal sent using a 4800 bits/s modem [25]. Because of the latency issue under WAMS as well as in wider are control systems (WACS), PMUs offer time stamp measurements. PMUs allow several different state estimations to be integrated into a complete set of state estimates of an area [26]. In the deregulated market, the system operational conditions may change quickly and dynamic power flow patterns appear to the system operator [27]. Reference [28] discusses the ability to capture system dynamics using WAMS in the case of the North East blackout of August 14, 2003. A conclusion is that WAMS offer better characterization than digital fault recorders. The real need for capability to capture dynamic system data relates to system control (i.e. WACS). Control signals need to be sent within one cycle of the disturbance to effectuate system control [29]. 7 Figure 1.2 is intended to depict the roles of measurements versus control. With WAMS becoming increasingly used, researchers have begun examining the concept of (WACS). Some potential elements of WACS are depicted in Figure 1.2. Such studies include the use of WAMS to control power system stabilizers [30,31]. Other system stability controls are being researched also as seen in [5, 27,32]. Figure 1.2 Pictorial of the transition of WAMS to WACS 1.6 Voltage Stability Limits Many large electric utility companies use nomograms for line loading based on voltage stability. As an example, American Electric Power (AEP) focuses on key 345 kV circuits, e.g., Kanawa River – Matt Funk. The voltage stability assessment method develops loadability limits on this line based on the status / levels indicated in Table 1.3. 8 Table 1.3 Status and levels monitored for the calculation of voltage stability limit calculations Status or level monitored Shunt reactors Series capacitors Outages in effect Generator status Load levels Season Typical number of cases considered, based on AEP practices 3 possible settings are precalculated 7 settings about 19 possibilities 3 settings 2 settings summer, winter With reference to Table 1.1, the precalculation gives in the order of 24,000 cases. If the number of cases in Table 1.1 is accepted as typical, then the 24,000 cases is simply obtained as the product 3 * 7 * 19 * 3 * 2 * 2. The table lookup is somewhat automatic coming from SE outputs. The operator sees the loadability limit. The |V| limit is generally 0.920 per unit, although there are more stringent limits (e.g., 0.940) at some generation busses. The base case is established by the MMWG working group, and this is updated seasonally or after a topological change. The method is applied in detail to this 345 kV circuit, and some form of the method may be applied elsewhere in the system. The operator does not exceed the loadability limits. If an event occurs which causes 85% of the limit on the line loading of the Kanawa River – Matt Funk line, alarms and transmission loading relief (TLR) measures go into effect. The cited TLR measures begin at a simple alarm if loading on the key circuit exceeds 85%. Also, TLR STEP 1 goes into effect and series capacitors may be switched. STEP refers to a special transmission emergency relief measure. At 95%, TLR 2 through TLR 3a may be used. This is a PJM (AEP’s ISO) curtailment measure. At 100% of loading, TLR 3b goes into effect. Above 100% loading, load shedding and higher TLRs may be used – including rolling blackouts above 110% line load. At the higher levels of TLR, the operator has full responsibility for these measures, including release of firm power sales. The nomogram tables are rechecked as needed, and the PTI software called PSSE is used. AEP also has TSAT but they do not use this often. The modeling is basically steady state, although transient stability checks are used on some occasions as a check. The |V| limits are more operationally restrictive than the angle stability limits. This seems to be basically true throughout the AEP system with at least one noteworthy exception. 9 The concept used considers the N-1 criterion. Operating conditions far (e.g., more than 2 busses away) from the line under study do not seem to effect the loadability of the line, and therefore this is not considered. 1.7 Angular Stability Limits Angular stability refers to the stability of bus voltage phase angles. That is, the phase angles do not loose synchronism after a disturbance. The calculation of maximum phase angle swings is central to the concept of angular stability. The approach taken in an off-line pre-calculation of stability entails the establishment of limits on operating conditions. These limits are collectively referred to as ‘angular stability limits’ for this purpose. The subsequent discussion relates to the calculation and use of angular stability limits at a major electric utility company. The majority of these remarks come from a visit to AEP. Some AEP generating plants are near the system’s edge, and there is a stability limit issue in these cases. At least one such case is handled using a nomogram: the Rockport generating station in Southern Indiana. This plant basically has two 345 kV transmission circuits leaving, and the lines are long. This results is angular stability limits. AEP considers this stability to be plant stability since it basically refers to the generation level at a single generating station. For the cited case, a number of prefault operating conditions are listed, and a table is formed using these conditions as rows, and fault type as columns. Again, the calculation is done off line and a nomogram is generated in tabular form. The prefault conditions can be accessed using SE and SCADA information, and the table gives the maximum curtailed generation limit. A quick inspection of typical limits shows that some curtailments can be in the 25% range – that is, the plant is operated only to 75% of capacity. Of course, these curtailed limits go into effect only for N hours per year (N may be in the range of 100). Studies of the number of megawatt hours curtailed have been done. With the exception of unusual unmodeled factors, the method has never resulted in an unstable operating condition for, perhaps, 20 years. A very quick assessment seems to imply that the method used is conservative, and may result in excess curtailed megawatt hours – but the conclusion is uncertain. Like the voltage stability nomogram, the angular stability nomogram is recalculated and checked on a bi-annual basis, and after any change in topology. In addition to the nomogram for angular stability and generation curtailment, there is a logic tree that is implemented to arm relays in case of emergencies. Fast valving is armed at all times. Other measures (relays) are armed in what appears to be special protection schemes based on the operating condition versus the nomogram dictated limits. 10 Some other generating plants have similar curtailment limits, but the detail of the calculation is simpler than in the case of Rockport. There was no significant loss of AEP generation or line outages in the August 2003 blackout. 1.8 Additional Considerations The general subject of power system stability, whether assessed on-line or offline, entails a large number of variables. Necessarily, in an actual operating environment, the electric utility companies must solve the practical problem of pre-calculating operating limits, and then assessing what operations would cause compliance or noncompliance with these limits. For purposes of highlighting the practical aspect of power system stability assessment and presently practiced, a series of questions were posed to electric utility companies. The following is a collection of these questions along with the responses. Relating to setting operating limits a. How operating limits are set in each system for transient stability? In particular, what kind of off line calculations are done, how many studies are done and how are they generated, what software is used, how is system reconfiguration and topology changes accommodated. The response is mainly that of AEP. AEP uses the angular stability nomogram for the Rockport plant as described above, and they use simpler but similarly generated curtailment levels for other stations. The calculations of the tabular nomogram are done off line using PSSE. This is required by FERC 715. System reconfiguration is accommodated as needed. Because the main application is at Rockport, with only two 345 kV lines as the main exit points for power, there is a finite and small number of cases of expected topology. b. What is the basis for setting operating limits for transient stability? In particular, are nomograms used, are these hard copy or in the computer. There seems to be no firmly consistent method cited for setting operating limits for transient stability nationwide. At one electric utility company, the basis was said to be that the operating power level of the generating station under consideration, and the accommodation of several different topologies and fault combinations are used in off-line studies. This is a ‘brute force’ approach that is suitable for off-line calculation. The tabular nomogram is a numerical table, precalculated, and made available to the operator through the EMS. 11 Relating to the utilization of base cases c. How does the operations planning base case relate to the system available in the EMS? At AEP, the MMWG base case is used throughout. Relating to contingencies d. What are the critical contingencies normally considered? ranked? How many are considered? How are these The general practice in North America appears to be the use of N-1 contingency studies, with contingency ranking used. Commercial software packages are used to obtain the contingency ranking. In the case of known critical circuits, additional contingency studies are used. As an example of the latter, at AEP, for the tabular nomograms developed for Rockport, all the contingency cases are studied. However, in general, N-1 contingencies are ranked (variable number – depending on the severity of the contingency) for studies. The role of PMUs e. What is the feeling of the engineers on the use of PMUs, how many are in the system now, how are the measurements brought back, how is the information used, is any three phase detail brought back (or is it all positive sequence)? Phasor measurement units are rapidly being deployed in the Eastern Interconnection in North America. These devices have a remarkable range of outputs including phase angle measurement (in three phase detail). Several eastern and western U. S. utility companies indicate that they have several PMUs and only positive sequence data are brought back. Perhaps the experience of AEP is typical in this regard: AEP has a long history (back to at least 1980) of PMU use – directly back to the inception of the device by Dr. Arun Phadke. At the moment, five PMUs are in place, and the data are brought to a data concentrator. This is done through a WAN at the rate of 30 samples per second. These data are made available to the operator through the EMS. One engineer at AEP compares PMU data with SE data – presumably this is the only actual use at the moment. The data are retained for off line analysis in a circular buffer. Only positive sequence data are brought back. The general feeling among the engineers is that PMUs are reliable and accurate. There is a fuzzy general plan to increase their use and bring back data, after some local processing, via the Ethernet. That is, the SCADA system shall be used. They are seeking advice on the suitability of the present sampling rate, namely 30 samples / second. 12 1.9 On-line Transient Stability Assessment With the advent of deregulation in the power industry and with the lack of transmission investment, today’s power systems are operated much closer to their limits. The stressed system conditions mean that operators are faced with scenarios that never occurred during the past. Various critical decisions need to be taken in real time and this requires sophisticated software tools. One of the main issues that needs to be tackled deals with online dynamic security assessment and control. The objective of dynamic security assessment and control is to ensure that the system can withstand unforeseen contingencies and return to an acceptable steady state condition [33] without transient instability or voltage instability problems. In this report the development of a software tool that uses artificial intelligence (decision trees) has been discussed. The proposed tool is developed by training a set of trees based on simulations conducted offline. With advances in computer technology it is now possible to create and store large databases that can be used to train agents. In this case the agents are the decision trees and the decision tree building procedure identifies critical attributes and parameter thresholds. In order to evaluate stability limits in terms of critical interface flows and plant generation limits, the use of synchronized measurements from PMUs has been proposed. 1.10 Overview of Analysis Techniques The problem of online transient stability assessment has been studied extensively over the last few decades. Various methods have been suggested in order to help in maintaining system security in real time. These can be broadly classified into efficient time domain simulation algorithms, direct methods and artificial learning approaches. [34]. Fig. 1.3 illustrates this classification. 13 Figure 1.3 Classification of on-line transient stability assessment methods ( taken directly from [34]) 1.10.1 Traditional time domain simulation Time domain simulation is the most accurate means of determining whether a system will remain stable following a particular contingency. The differential equations of the system are solved step by step using numerical techniques in order to get the actual values of the state variables. These state variable values then yield important information regarding transient stability. Typically the machine swing curves are computed in order to get an estimate of the maximum deviation of machine relative rotor angles. The system dynamics are simulated for the faulted and post fault period. The post fault period is generally of the order of a few seconds during which the system behavior is observed. However, this method is very time consuming and does not yield any information regarding stability margins [33]. 1.10.2 Direct methods The direct methods consist of using Lyapunov’s functions in order to assess stability [34]. The Transient Energy Function (TEF) is one such Lyapunov function that can be used to assess stability. It can be thought of as a multi-machine equivalent of the more simplified equal area criterion. The basic idea consists of evaluating the kinetic energy 14 gained by the machines during the fault period and comparing it with the maximum gain in potential energy that the system can withstand (without becoming unstable) after the fault is cleared as shown in Fig. 1.4. Figure 1.4 Illustration of the Transient Energy Function (TEF) method The main advantage of these direct methods is the significant reduction in computing time. Also these methods give an idea of the stability margin of the system and hence help in evaluating the proximity of the system to instability. 1.10.3 Artificial learning With advances in computer technology coupled with the decreasing cost of mass storage devices, it has become possible to create large databases that can be used to train agents. The underlying principle in this method can be described as follows [35]: Given a set of preclassified examples (learning set) along with their characteristic attributes, derive a general rule (using thresholds) that can explain the input-output relationship of the preclassified cases and correctly classify new or unseen cases.’ In the context of transient stability assessment, the database would consist of numerous cases that have been assessed using time domain techniques. The attributes would consist of parameters such as generator outputs, critical line flows and relative phase angles. The output could be a stability margin or a simple classification (secure/insecure). Artificial neural networks (ANN’s) The basic building block in a neural network is a neuron or a perceptron [35]. The neuron is a linear threshold unit that can be used to model simple linear constraints. A neural network consists of many such neurons arranged in different layers. The learning process consists of identifying the structure of the network, selecting various parameters and defining thresholds. 15 Statistical pattern recognition (K-NN) This artificial learning technique consists of matching an unseen situation with the examples stored in the database. The unseen state is assigned the class of that preclassified example whose attribute values most closely match its own. The K-NN method consists of classifying a state into the majority class among its k nearest neighbors in the stored database. Hybrid techniques that involve coupling of K-NN with decision trees for power system transient stability assessment are dealt with in [36]. Machine learning (Decision trees) Machine learning is a field of automatic learning that consists of building if-then rules similar to that used by humans. This report focuses on the use of decision trees (DTs) for online transient stability assessment. Decision trees belong to the class of techniques which deal with machine learning from examples. One of the first attempts at a detailed discussion of applying decision trees to assess the transient stability of a power system can be found in [37]. The paper lays the foundation of the inductive inference method that can automatically build decision trees based on a preclassified set of examples. The rules thus obtained relate the system stability to pre fault static parameters of the system. A deductive inference method is then used to classify unseen cases in real time. 1.11 Decision Trees for Transient Stability Assessment Transient stability assessment deals with analysis and preventive control. Since the DT based method uses pre fault parameters in order to assess stability, it can be used as an effective approach for preventive control. The use of controllable parameters in decision rules often becomes an important aspect while taking preventive action. In this regard, stability constrained optimal rescheduling of generation is an area that has received considerable attention. Reference [38] deals with the use of the transient energy function (TEF) method of transient stability assessment in order to carry out rescheduling of generation and critical line flows for a given initial operating condition and specified contingency. Reference[39] uses DTs as a tool for security constrained generation redispatch. The paper takes an interesting approach and uses a two step procedure. The first step consists of running an economic dispatch on all generators without considering security constraints. PT = P1 + P2 + P3 + .................. + Pi (1.1) Following this, the generator set points are passed through a decision tree which checks the various generator outputs. In case there are no violations the results of the economic dispatch are displayed to the operator. However, in case one of the generators violates a decision rule a generation redispatch is carried out as follows: the decision rule 16 used in the last splitting node of the DT is used as an equality constraint in order to fix the generation of the corresponding unit. Equation (1.1) now becomes PT − PU = P1 + P2 + P3 + .................. + Pi −1 (1.2) Thus the total generation now available for dispatch changes to PT - PU and this has to be distributed economically amongst i-1 generators. The above step is repeated iteratively till a secure solution is obtained. The essential features of the decision tree transient stability method are documented and discussed in detail in [40]. The paper deals with the selection of splitting and stop-splitting criteria that are used when building a DT. The building of the learning set and its effects on the accuracy of the method is analyzed. Questions related to optimal size of the learning set, selection of candidate attributes and building of multi-class DTs are dealt with. Tradeoffs related to complexity versus number of classes and computational issues relating to the building and deployment of DTs in online environments are discussed. A comprehensive discussion of the method and its critical aspects is provided. The answers provided to different basic questions about the method make it an important piece of literature in this field. The decision tree transient stability method is further developed and deals with the practical aspects regarding deployment. Different methods for enhancing the reliability of the DTs and for generalizing them to multi-contingency stability assessment are presented. Unlike [39-43], [44-46] tackle the problem of emergency control using decision trees. Instead of using pre fault static parameters of the system as candidate attributes for building the DT, post fault measurements from a fault that is currently in progress are used. Using this information, the decision tree is used to predict whether the system is moving towards transient instability. Based on the prediction appropriate emergency control action can be initiated. 1.12 Report Organization This report focuses on using decision trees for transient stability assessment by using pre-fault static parameters of the system. The critical thresholds obtained can then be used for preventive control. Section 2 provides an overview of the decision tree methodology. The relevant theoretical aspects of the method are discussed in detail. The steps involved in building a decision tree for transient stability prediction are enumerated and described. Since this 17 work uses the Classification and Regression Tree (CART) algorithm [47] for DT building, an overview of the salient features of this algorithm is also provided. Section 3 discusses the three stages of the proposed DT building scheme. The process of offline DT building, periodic DT updating and the use of DTs in an online environment is dealt with in detail. Certain issues related to DT reliability and preventive controls are discussed. Section 4 deals with the case study carried out on the Entergy system. Based on the DT generated for prediction of damping problems a nomogram identifying insecure regions in the attribute space is generated. A new path based classification approach is proposed and tested. 18 2. Theoretical Background 2.1 Definition A decision tree can be thought of as a flowchart representing a classification system [website need to fill this]. It consists of a sequence of simple questions regarding critical attributes (CAs). The answers to these questions trace a path down the tree. The classification made by the decision tree depends on the class of the terminal node. In other words, a decision tree can be thought of as a series of tests that are conducted on attributes of an object in order to ascertain its class. In Fig. 2.1, A1, A2 and A3 are the critical attributes. For numerical attributes A1 and A2, the question compares the value with a threshold. For categorical attribute A3 the question checks if it belongs to a particular set. With respect to power systems, A1 and A2 could be thought of as generator outputs (numerical attributes) and A3 would be the fault location or fault bus number (categorical attribute). Finally, a class is assigned to each object or case at the terminal node. In this report the Classification and Regression Trees (CART) software is used to build decision trees. CART is a decision tree procedure introduced in 1984 by Stanford statisticians Leo Breiman, Jerome Friedman, Richard Olshen and Charles Stone [47]. A short description of the CART methodology is given below. 2.2 CART Methodology The procedure used by CART is called binary recursive partitioning. The term binary refers to the fact that CART splits each node into two child nodes. This splitting procedure is then applied to the child nodes and so on. The main aspects of tree building that CART deals with are: 1. 2. Splitting criteria Deciding when to stop growing the tree 2.2.1 Splitting criterion In order to split a node, CART asks questions that have a yes/no answer. Based on the answer the cases are divided into two parts. There are several questions that can be asked regarding the value of a particular attribute. Consider a simple example having 100 cases with 5 attributes. Assuming that every attribute has a unique value for each case, CART will consider 100 candidate splits for each attribute. Therefore, at each step of tree building CART evaluates 500 possible splits and chooses the best split based on certain criteria. 19 Root node Numerical attribute Categorical attribute Terminal node Figure 2.1 Example of a decision tree CART ranks candidate splits using a measure called node impurity [47]. The target variable in this case is the security of the contingency and can take on two values (‘Secure’ or ‘Insecure’). Each node of the tree is assigned a unique number n and will have an estimated probability ( p( S | n) or p( I | n) ) for secure and insecure cases. The Gini impurity criterion is used in this project. For node n this criterion is given by: I (n) = 1 − ( p( S | n)) 2 − ( p( I | n)) 2 (2.1) If node 1 has cases with probabilities of p( S | 1) = 0.8 and p ( I | 1) = 0.2 then the value of the Gini diversity index is: I (1) = 1 − (0.8) 2 − (0.2) 2 = 0.32 (2.2) 20 If any node consists of only one class then the probability of that class will be one and the probabilities for the other classes will be zero. This will result in a diversity index of zero which suggests that the particular node has no diversity. Using the Gini criterion, CART ranks the candidate splits and then selects the best one amongst them. The question with the highest score is selected and is called the critical splitting rule (CSR). The parameter used in the critical splitting rule is the critical attribute. Other rules or questions that have a lower score are called competitors. Rules that mimic the action of the critical splitting rule (i.e. have the same score) are called as surrogates. A competitor can generate a good split whereas a surrogate generates exactly the same split as the CSR. As the tree grows the nodes become more and more homogeneous. 2.2.2 Stop-splitting criterion and DT pruning Many of the commercially available decision tree software packages use a goodness-of-split criterion [47] in order to decide when to stop splitting. If all the splits at a particular node do not meet the criterion then that node is considered to be a terminal node. When all branches of the tree reach terminal nodes, the process of growing the tree is considered complete. CART uses a different procedure in order to arrive at the best tree. It uses the data from the learning set in order to build the largest possible tree. This tree is then pruned to generate several sub-trees which are tested using the data from the test set. A commonly used index is the misclassification cost which is calculated using 1 ts R ts = ts ∑ c(i | j ) ⋅ N ij (2.3) N i, j where N ts is the total number of cases in the test set, c(i | j ) is the cost of misclassifying t a class j case as a class i case, N ijs is the number of class j cases that have been misclassified as class i. In order to avoid misclassification of insecure cases, the cost of misclassifying an insecure case as secure is higher than the cost of misclassifying a secure case as insecure (false alarm). Unlike other decision tree software packages CART does not test the splits at each node and therefore avoids missing important information that may be found using further splits. Finally, the decision tree with the least misclassification rate is selected. Other criteria such as size or characteristics of critical attributes may be used to select an optimal tree. 21 2.3 Decision Tree Building (preliminary tasks) Before building a decision tree for transient stability analysis the following tasks need to be carried out [40]: 2.3.1 Generation of operating scenarios The building of the learning set is the most important aspect of the tree building process. In order to capture the essential features of the power system under study, an exhaustive list of operating conditions needs to be generated. Each operating point is described by its characteristic attributes such as generation outputs, line flows and network topology. The generation of different operating scenarios can be done by varying the generation-load profiles within a given range and/or changing network topology by simulating generator and line outages. A detailed discussion regarding the generation of different operating conditions is provided in Section 4 for the Entergy system. 2.3.2 Stability assessment After systematically generating the operating points, various stability tests need to be performed on each operating point. A list of credible contingencies is obtained from past experience or current system status. Each contingency is then analyzed using transient stability assessment software and the results (i.e. ‘Secure’ or ‘Insecure’) are stored. This contingency analysis is carried out for all the generated operating points. 2.3.3 Selection of attributes The selection of candidate attributes should be done such that they adequately characterize the properties of the system. Based on past experience certain parameters such as generator power (real and reactive) outputs, power flows on critical lines, phase angle differences and information about network topology could be used. Since the DT would be used for control purposes, the attributes selected should be controllable. Selection of simple controllable attributes will yield trees that are easy to interpret whereas the use of more complex attributes will result in trees that are more accurate and suitable for the purpose of overall system analysis. 22 2.3.4 Preparation of database The selected attributes and results of the stability analysis for each contingency (for a particular operating condition) are systematically stored in the database. This step is then repeated for all the operating conditions. Therefore each case will be characterized by a certain number of critical attributes and a predetermined (from stability analysis) stability class. Table 2.1 shows a sample database. This set of cases is then separated into a learning set and test set. The learning set (LS) is used to train the DT and the test set (TS) is used to evaluate the performance of the decision tree on unseen cases. Generally a randomly selected fraction of the entire database is used as the test set (20-30 % of all cases). Table 2.1 Sample database Class Secure Insecure Secure Secure Insecure …….. 2.3.5 Fault bus 97745 97747 98654 98634 98767 ……… Other bus 97746 97750 98656 98632 98… ……… P_98537 _98606 (MW) 45 60.5 63 50 ……. ……… PG_98233 (MW) 100 100 110 110 …… …… A_98606 (deg) 15 17 18 17 …… …… Multi-class and multi-contingency DT considerations The number of stability classes to be considered while building the DT needs to decided in advance. It should be kept in mind that the learning set should have an adequate representation of all the considered stability classes. The decision tree performance depends on the quality of the learning set and if sufficient number of cases belonging to every stability class are not present in the database the accuracy and reliability of the DT in classifying unseen cases will be poor. Also, decision trees can be built for a single contingency or for multiple contingencies. In case of multi-contingency decision trees, attributes related to the fault location and type of fault can be used as predictors. The multicontingency DTs will give a list of contingencies that will result in instability for a given operating condition. On the other hand, a single contingency DT will provide information about which operating condition will make the system unstable in the event that the particular contingency (for which the DT is built) occurs. 23 2.4 Use of Phasor Measurement Units for Transient Stability Assessment The system operator is generally provided with a set of operating guidelines that need to be followed in order to maintain secure system operation. These guidelines are obtained by carrying out extensive offline studies. Critical parameters are identified and thresholds for these parameters are evaluated. The operator is then required to maintain system operation in a secure region based on these constraints. Traditionally operators rely on SCADA measurements and the state estimator in order to ensure that the system is operating within the secure region. However since the values obtained from SCADA measurements and state estimators are not synchronized, it is difficult to accurately measure the current operating condition. The main advantage of using PMUs is that they provide GPS-synchronized measurements that can then be used to calculate phase angles, frequency and rate of change of frequency. The time synchronism that can be achieved with the use of PMUs makes it possible to obtain a snapshot of the system condition at any given point in time. Operational nomograms are built offline in order to assess stability in real time. If the current operating condition falls in a particular region, then an accurate assessment of the system stability can be made by using the stored nomograms. The effectiveness of this method depends upon the accuracy and reliability with which operating conditions can be measured. This is where PMUs can be of great help. With accurate synchronized measurements of phase angles and critical power flows throughout the system, PMUs can give system operators the precise information needed in order to characterize the current operating point. 24 3. 3.1 Proposed Scheme DT Based Online Security Assessment Scheme A flow chart describing the DT-based online security assessment scheme is given in Fig. 3.1. There are three main steps involved in building the DT and using it for online security assessment. Figure 3.1 Flowchart describing the proposed scheme 25 The offline DT building consists of creating a large database of cases along with their classification by time domain simulations. This database is then used to build a DT. Depending on the system conditions the DT built offline can be updated on an hourly basis. The updated DT is then used online for transient stability prediction with the help of accurate synchronized measurements from phasor measurement units. 3.1.1 Offline DT building This stage is executed offline to generate a database of cases and build a DT for a 24-hour (day ahead) horizon. Initially, NOC prospective operating conditions (OCs) for the next 24 hours are obtained from short-term load forecast and unit commitment programs. These operating conditions are selected such that they adequately represent probable system behavior with regard to generation and load profiles and network topology. NC contingencies based on type, location and fault duration are either assumed throughout the power system or selected by the operator from a history of critical contingencies. For each OC, detailed time-domain simulations of all the NC contingencies are conducted. Specified security criteria dealing with transient stability, transient voltage dip, transient frequencies, and damping ratios are then checked to determine the security classification for each case. The classifications may be binary (i.e. ‘Secure’ or ‘Insecure’) or multiple, e.g. ‘transient instability’, ‘insufficient damping’, ‘voltage insecurity’, and ‘frequency insecurity’ with different priorities assigned to each criterion. Thus, the class of an insecure case is determined by the criterion with the highest priority. Finally, a database of NC × NOC cases each including a security classification and a vector of attribute values is generated and further divided into a learning set and a test set. A group of predictors is then selected from either fault-dependent (FD) attributes, e.g. the fault type and location, or fault-independent (FI) attributes, e.g. bus voltage angles, MW transfers across lines or interfaces, outputs of generators, etc. Some predictors critical to the system’s security are screened out as critical attributes (CAs) to form the critical splitting rules (CSRs) of the DT. Fault-independent critical splitting rules (FICSRs) are particularly of interest because they define the insecure regions in the attribute space where the current OC can be located. To exactly determine the location of the OC, it is generally necessary to obtain accurately synchronized values of fault-independent critical attributes (FI-CAs). A number of PMUs have been installed in many power systems and they provide accurately synchronized measurements of frequently changing parameters, including the voltages of the buses where they are installed and the currents and power flows of the branches connected to the bus. If impedances of the branches are constant, then voltage angles at the other ends of the branches can also be acquired. Thus, the synchronized measurements from PMUs are good candidates for FI predictors. In addition, predictors can be selected from the parameters monitored by the SCADA system if they do not change frequently, e.g. outputs of generators. In this case the measurements obtained need not be synchronized. 26 Following this, classification tree algorithms are used to build a DT from the learning and test sets. CAs and CSRs are selected from all available attributes and are stored in each inner node. The DT built could be two-class (“insecure” or “secure”) or contain N + 1 classes if N insecure classes are considered. If storage space allows, even N two-class DTs could be built at the same time, each of which focuses on a different type of insecurity. 3.1.2 Periodic OC prediction and DT updating The time horizon is divided into periods of equal length (typically 1 hour). The DT update procedure is executed for each period in order to predict the performance of the DT for the upcoming time period. In case the performance of the DT is not satisfactory, a new DT can be built. Prospective OCs in the next period can be predicted using a short-term load forecast (this could be more accurate than the forecast used previously) and unit commitment rules. If the prospective OCs are close to any of the NOC OCs already considered, the DT will remain unchanged until the next period. If new OCs appear the DT may not perform well. Thus, the NC contingencies will be simulated again at the new OCs to generate new cases, which are used to test the existing DT. If the DT performance is still satisfactory the DT is retained until the next period in the time horizon. Otherwise, the new cases together with the old are used to build a new DT. The speed of building a DT from scratch is discussed below. For the case studied in the next section, 280 ‘N – 1’ contingencies are considered for 56 OCs of the online EMS representation of the Entergy system. Thus, 15680 simulations are needed. Using a PC with a single Pentium-4 3.4GHz CPU, each simulation takes 5~10 seconds and the 280 simulations for each OC are completed within 50 minutes. Since all simulations are independent, they can be executed on a system with multiple parallel CPUs to accelerate the computation time. Finally, CART’s growing and pruning processes are extremely fast such that an iterative DT building procedure may be designed to quickly generate a good DT from the simulation results. As compared to building a new DT, updating an existing DT takes much less time. This is because the simulations need to be carried out for new OCs only and a DT needs to be built from the updated database. 3.1.3 Online security assessment In real time, the control center obtains synchronized measurements of various parameters throughout the system. If the FI-CAs used to build the DT are selected from existing PMU measurements then the operator can perform accurate dynamic security as- 27 sessment for a single contingency or a group of contingencies. To consider a general case, assume that N’C (1 ≤ N’C ≤ NC) contingencies are considered (NC denotes the total number of contingencies considered during simulation). The following procedure can predict security and give guidelines for preventive control: • Prune all branches of the DT that do not correspond to the contingencies of interest. If the operator is only interested in faults at a single bus (say ‘Bus 1’), then only that branch of the DT dealing with ‘Bus 1’ needs to be considered. In such cases the use of fault dependent critical attributes in the initial splits can prove to be helpful. • Use operational nomograms (constructed offline) to ascertain the location of the current operating condition in the attribute space. • If the DT predicts an insecure condition following a particular contingency or set of contingencies, arm appropriate preventive control measures. 3.2 Reliability of the DT When using the DT constructed for online security assessment it is necessary to understand the importance of each predictor. The reliability of a particular CSR depends on the number of cases in the node where it is used. A CSR for a parent node with a large number of cases is obviously more reliable than a CSR for a parent node with a smaller number of cases. Therefore, the first predictor carries the highest weight since it is derived by observing all the cases in the learning set. The predictors used in the later stages of the DT operate on a certain fraction of the data only and hence carry less weight. When operating conditions change, it is likely that the predictors used near the terminal nodes will lose their credibility. In such cases it becomes necessary to use a path based approach in order to ascertain the class of a terminal node. In the normal DT building procedure the class of the terminal node is assigned based on the proportion of secure/insecure cases in the terminal node. However, this does not take into account the path followed by a particular case before entering the terminal node. In order to account for the path traced by any case as it drops down the tree, an insecurity score (S) can be calculated for each path. S= ∑p i =1 K i =1 K i × ωi (3.1) i ∑ω pi ---- Percentage of insecure cases in node i ω i ---- Weight assigned to node i The insecurity score calculates the weighted average of the percentage of insecure 28 cases in the nodes belonging to each path. In Fig. 3.2 the insecurity score for path 1-3 equals the weighted average of the percentage of insecure cases in nodes 1 and 3. In order to classify a terminal node the insecurity score of the path leading to that node is calculated and compared to a threshold. A sample set of insecurity scores are calculated for the tree in Fig. 3.2. The threshold for the insecurity score used here is 30 %. All paths with scores greater than 30 % are considered to be insecure. N = 100 80 20 S I A1>50 A1<=50 N = 70 S I A2>34 65 5 1 S I A2<=34 A3>9 N = 30 15 15 2 A3<=9 N = 51 S I S 50 1 3 S I N = 19 15 4 I 4 S I N = 10 5 5 I 5 S I N = 20 10 10 S 6 Figure 3.2 Sample DT illustrating new path based approach S ---- Secure I ---- Insecure Table 3.1 Insecurity scores of paths for sample tree Path 1-3 1-4 2-5 2-6 Insecurity score 4.96 10.11 50.00 50.00 Predicted class (path based) S S I I Predicted class (terminal node based) S I I S Secure/Insecure 50/1 15/4 5/5 10/10 The prediction accuracy using the traditional terminal node based method (Table 3.2) and the new path based method (Table 3.3) is compared below. 29 30 Table 3.2 Prediction accuracy for traditional method Actual Class S I Total Cases 80 20 Percent Correct 75 45 S N=71 60 11 I N=29 20 9 Table 3.3 Prediction accuracy for path based method Actual Class S I Total Cases 80 20 Percent Correct 81.25 75 S N=70 65 5 I N=30 15 15 Using the path based method; better results are obtained for both secure and insecure cases. The path based classification approach increases the robustness of the DT to variation in operating conditions. This method when used on the Entergy system yields promising results. 3.3 Preventive Control After determining the contingencies that will result in instability (for a particular operating condition) appropriate preventive control action must be initiated if these contingencies occur. Using the decision tree the critical attributes whose thresholds are being violated can be determined. In Fig. 3.2, A1=50 is a threshold value. In order to make more cases enter node 1, the value of A1 can be made greater than 50. However, this approach does not always yield good results. Although the thresholds for the critical parameters are accurate to classify system security in real time, they need not be effective for preventive control. The DT in Fig. 3.2 suggests that A1 ≤ 50 results in insecurity. Assume that attribute A1 denotes the power flow through a line. Increasing the power flow of the line beyond 50 will not necessarily result in a secure operating condition. Given the fact that the power flow through a line can be changed in many different ways, the resulting operating condition might be a new one that has not been considered in the database. Therefore while designing preventive control measures it is necessary to ensure transition to a secure operating condition that has already been analyzed and proven to be secure. 31 3.4 Building a Better DT by Adjusting Penalties for Predictors Both FI-CAs and FD-CAs are used in the DT. Sometimes, FI-CAs are preferable since they are measurable and partially controllable. For the purpose of security analysis, it may be important to separate or change the order in which they are used in the splitting process. FD-CSRs could be used before FI-CSRs to identify insecure areas of the system (by location), or alternatively, FI-CSRs could be used before FD-CSRs to identify insecure operating regions in the attribute space and design preventive control. Either approach can be used by assigning appropriate penalties between 0 and 1 to faultindependent or fault-dependent predictors to reduce their improvement scores. A penalty of 0 and 1 respectively means “no limitation on using” and “not using” the predictor. The higher the penalty a predictor has, the later it is used in the DT. The following procedure can be used to build a DT that uses FD-CSRs after FI-CSRs by adjusting the penalties of fault-dependent predictors: Procedure-1: 1) Assign each fault-dependent predictor a penalty p = p min . 2) Generate a series of DTs by CART’s growing and pruning processes and select the best DT. The best DT is the one with the lowest misclassification cost. The child nodes formed as a result of using FD-CSRs are considered as the terminal nodes and the rest of the tree is pruned. Thus, a DT using FD-CSRs only at the last split is generated and stored. 3) Let p = p + Δp , if p ≥ p max , select the best DT from all DT’s stored and end the procedure; otherwise, go to 2). Heuristic techniques could be applied to find a good value for the penalty. Similar procedures could be designed to generate a DT meeting other requirements. Another approach to build DTs is to only select predictors from critical variables e.g. plant generation outputs of critical plants, power flows of critical interfaces or transmission lines, and phase angle differences between critical buses. It is important to note that given a set of predictors, CART’s DT building algorithm does not guarantee building an optimal DT. At each step of DT growing, i.e. splitting a node, the algorithm simply seeks the CSR that best separates cases of different classes only for that particular node. CART does not take into account the effect that the current split will have on the overall classification accuracy. CART generates a locally optimal split at each stage as opposed to a globally optimal split. Simulations have shown that replacing a CSR by a competitor could totally change the nature of the DT and result in better prediction for unseen cases. Finding an optimal DT given a set of predictors is a complicated problem, and results in a very large number of predictor combinations. However, the following Procedure-2 can be used to quickly find a good DT by recursively executing CART’s DT building algorithm. 32 Procedure-2: 1) Select a group of reasonable predictors and build a DT. Let i = 1. 2) For the CSRs of all nodes at the i-th level starting from the root node, check their non-surrogate competitors. 3) If a CSR at a node has a good non-surrogate competitor whose improvement score is larger than γ times (e.g. γ = 0.9) of its score, then assign its associated predictor the minimum penalty that ensures that the new DT uses the competitor instead of the old CSR at that node. The value of the ‘minimum penalty’ is easily determined by a simple trial and error procedure. Then, go to 5). 4) If no such good non-surrogate competitor exists at the i-th level, go to 6). 5) If the new DT built in 3) is better than the old DT, replace the old DT by the new DT and go to 2). 6) If no parent node exists below the i-th level, then end the procedure; otherwise, let i = i + 1 and go to 2). Using the above procedure an optimal DT for a given set of predictors can be generated. A better understanding of the concepts described in this chapter can be gained after going through the case studies on the actual Entergy system in the following section. 33 4. 4.1 Entergy Case Study System Description The proposed scheme was tested on the Entergy online operational model which consists of around 240 generators, 2000 buses and 2100 lines. In order to carry out the analysis in a systematic manner, the entire network was divided into five areas as shown in Fig. 4.1. These areas were chosen based on the information provided by Entergy regarding their system configuration. Figure 4.1 Entergy transmission system map 34 This case study focuses on the Amite South area. This area with 54-generators, about 400-buses and about 430-lines includes the city of New Orleans. A decision tree (DT) is built for a 24 hour horizon. The different aspects of DT building and analysis are dealt with in detail in the next few sections of this chapter 4.2 4.2.1 Database Generation Operating conditions The power flow data for two days was obtained from Entergy. The data consists of light load and heavy load power flow profiles for July 19, 2006 and July 26, 2006. The network topology and generation-load profile are considered on a day ahead basis for 24 hours of July 19, 2006. A credible stressed situation with a 500 kV line from Wells to Webre being out of service is considered. This is shown by the cross (‘X’) in figures 4.1 and 5.2. The Wells-Webre line connects the WOTAB area and Amite South area. Using this as the starting point, the light load (17819 MW) and heavy load (26600 MW) power flow profiles of July 19, 2006 were used to generate 54 intermediate operating conditions. Percentage scaling was used to vary the Entergy load and the generation was adjusted using certain economic dispatch rules in order to achieve generation load balance. Generators were switched in and out of service depending on a list provided by Entergy. While changing generation, it was ensured that minimum and maximum generation limits were not violated. 4.2.2 Contingencies For each of the 56 operating conditions the following contingencies in the Amite South area were considered: • Three phase faults on both ends of all 500 kV lines with a fault clearing time of 5 cycles • Three phase faults on both ends of all 230 kV lines with a fault clearing time of 6 cycles There are a total of 280 such contingencies for each operating condition. This results in a total of (280 × 56) 15680 cases. Twenty percent (20 %) of these cases (0.20 × 15680 = 3136) are used as the test set (TS) and the rest (12544) of the cases are included in the learning set (LS). Each case is simulated in Powertech Lab’s TSAT software [48] and two security criteria, namely, transient stability and damping are checked. Different DT’s are built for transient stability and damping problems using the CART software developed by Salford Systems [49]. The criteria used are as follows: • Transient stability criterion: After the fault is cleared the transient stability margin of the system should be greater than 5 % and no bus voltage should violate the voltage range of 0.70-1.20 p.u. for more than 20 cycles. Here, the stability margin is es- 35 timated by TSAT’s power swing based algorithm [48]. • Low damping criterion: Entergy requires the damping ratio to be greater than 3 %. This criterion is applied to the inter-area oscillation modes of generator rotor angles (frequency range is 0.25-1.0 Hz). Figure 4.2 Amite South area (showing PMU locations) 4.2.3 Candidate attributes Entergy has installed 9 Phasor Measurement Units (PMUs) to monitor critical interfaces and get a good visibility of the system. The Amite South area has 3 PMUs at buses Waterford 230kV and Fancy Point 230kV (encircled in Fig. 4.2). The PMUs measure voltages of the buses and currents of the following lines: 1) Waterford 230 kV (bus 98537) looking at line to Ninemile 230kV (bus 36 2) 3) 98606) Waterford 230kV (bus 98537) looking at Waterford Auto 500kV (bus 98539) Fancy Point 230kV (bus 98234) looking at Fancy Point Auto 500kV (bus 98233) Therefore the voltage phase angles at bus numbers 98233, 98234, 98537, 98539 and 98606 can be obtained from the PMU measurements. Also the line flows of the three lines mentioned above are available from the PMUs. In order to obtain additional information regarding the potential locations of new PMUs certain additional measurements from other 500 kV buses in the area are also included in the list of critical attributes. There are a total of 8 other 500 kV buses in the area. Two groups of predictors are used to build the decision trees. Group-1 (from candidate PMU measurements): 1) The serial number of the fault bus (FB); 2) MW-flows of the 16 branches (each denoted by P_i_j, where i and j are two bus serial numbers); 3) Angle differences between voltage phase angles of 13 buses (each denoted by A_i_j, where i and j are two bus serial numbers). Group-2 (from existing PMU measurements): 1) The serial number of the fault bus (FB); 2) MW-flows of the 3 branches where PMUs are currently installed (P_98537_98606, P_98537_98539 and P_98233_98234); 3) Angle differences between voltage phase angles of 5 buses currently being monitored by PMUs (98233, 98234, 98537, 98539, 98606). 4.3 Results The classification accuracy obtained from each group of predictors for decision trees built for the transient stability and damping criteria are as follows: 37 4.3.1 Decision trees for transient stability (DT-U) Group-1 (candidate PMU measurements, DT-U1) Table 4.1 Prediction accuracy for Learning Set (DT-U1) Actual Class S U Total Cases 12,263 281 Percent S U Correct N=11722 N=822 95.588 11,722 541 100.000 0 281 Table 4.2 Prediction accuracy for Test Set (DT-U1) Actual Class S U Total Cases 3,062 74 Percent S U Correct N=2943 N=193 96.048 2,941 121 97.297 2 72 Group-2 (actual PMU measurements, DT-U2) Table 4.3 Prediction accuracy for Learning Set (DT-U2) Actual Class S U Total Cases 12,263 281 Percent S U Correct N=11722 N=822 95.588 100.000 11,722 0 541 281 Table 4.4 Prediction accuracy for Test Set (DT-U2) Actual Class S U Total Cases 3,062 74 Percent S U Correct N=2943 N=193 96.048 2,941 121 97.297 2 72 38 FB$ FB$ P _98246_98539 P _98537_98606_P U M P _97301_98430 P _98537_98606_P U M P T_98233_98234_P U M FB$ FB$ FB$ A_198234_198539 FB$ P _98233_98235 P _98235_98390 FB$ P _98537_98606_P U M Figure 4.3 Details of predictors used at each stage of the DT-U1 39 FB$ FB$ P T_98537_98539_P U M P _98537_98606_P U M A_198233_198539 P _98537_98606_P U M P T_98233_98234_P U M FB$ FB$ FB$ A_198234_198539 FB$ P T_98233_98234_P U M P T_98233_98234_P U M FB$ P _98537_98606_P U M Figure 4.4 Details of predictors used at each stage of DT-U2 40 4.3.2 Decision trees for damping (DT-D) Group-1 (candidate PMU measurements, DT-D1) Table 4.5 Prediction accuracy for Learning Set (DT-D1) Actual Class D S Total Cases 2,015 10,529 Percent D S Correct N=2991 N=9553 99.454 2,004 11 90.626 987 9,542 Table 4.6 Prediction accuracy for Test Set (DT-D1) Actual Class D S Total Cases 486 2,650 Percent Correct 97.942 90.981 D S N=715 N=2421 476 10 239 2,411 Group-2 (actual PMU measurements, DT-D2) Table 4.7 Prediction accuracy for Learning Set (DT-D2) Actual Class D S Total Cases 2,015 10,529 Percent D S Correct N=2921 N=9623 98.958 1,994 21 91.196 927 9,602 Table 4.8 Prediction accuracy for Learning Set (DT-D2) Actual Class D S Total Cases 486 2,650 Percent Correct 96.708 91.132 D S N=705 N=2431 470 16 235 2,415 41 P _98246_98390 P _98246_98390 A_98390_198233 A_98246_98487 Surrogate P_98246_98390 P _97301_98430 P _97301_98430 FB$ A_98487_198233 P _98537_98606_P U M FB$ P _98537_98606_P U M FB$ FB$ FB$ FB$ Surrogate P_98246_98539 Figure 4.5 Details of predictors used at each stage of DT-D1 ‘R’ denotes the root node ‘A-N’ denote non terminal nodes ‘1-16’ depict terminal nodes 42 PT_98537_98539_PMU R A 1 P_98537_98606_PMU B PT_98537_98539_PMU PT_98537_98539_PMU C D FB$ P_98537_98606_PMU E G J 2 H A_198233_198537 F 6 L FB$ 3 4 5 I FB$ PT_98537_98539_PM U P_98537_98606_PMU K FB$ P_98537_98606_PMU 7 8 9 10 13 M FB$ N FB$ 11 16 12 14 15 Figure 4.6 Details of predictors used at each stage of DT-D2 43 Fancy Point 230kV (98234) Coly 500kV (98390) Willow Glen 500kV (98246) Waterford 500kV (98539) Waterford 230kV (98537) Nine Mile 230kV (98606) Figure 4.7 Key transmission path 44 The predictors denoting line flows used by DT-D1 and DT-D2 are given by P_98246_98390, P_98246_98539 = P_98539_98537, P_98537_98606. Note that although Fig. 4.5 shows P_97301_98430 as a predictor, the power flows of the key transmission path shown in Fig. 4.7 generate equally good splits. Using the thresholds for the critical attributes identified in DT-D2 (Fig. 4.6) the following operational nomogram is constructed (Fig. 4.8). P-98537-98606 0.15 0.43 256.25 0.55 18.75 0.09 -366.9 -289.85 55.3 268.1 PT-98537-98539 0.30 0.0 0.89 Figure 4.8 Operational nomogram built using DT-D2 The numbers written in the various regions denote the probability that an operating condition lying in the region is insecure. The darker shades depict regions of high insecurity. This nomogram can be used in real time to evaluate the state of the system and arm appropriate control action. For example, if the operator finds that PT_98537_98539 ≥ 55.3 MW and P_98537_98606 ≤ 18.75 for the current operating condition then the chances of any line outage (on 230kV and 500kV lines) resulting in damping problems equals 89 %. 46 4.4 Evaluation of DT reliability In order to assess the performance of the constructed DT’s on unseen cases, further tests were carried out using the data for July 26, 2006 (a week after July 19, 2006). In order to generate different operating scenarios a procedure similar to the one used for the July 19, 2006 cases was followed. Time domain simulations were carried out for contingencies on all lines above 230kV. The results were stored in a database. The DT-D2 was tested on this new set of data. The prediction accuracy obtained is as shown in Table 4.9. Table 4.9 Prediction accuracy for July 26 cases (DT-D2) Actual Class D S Total Cases 1,058 5,614 Percent D S Correct N=1875 N=4797 75.90 803 255 80.90 1,072 4,542 Using the new path based approach proposed in this report the following paths were identified (Table 4.10) and their insecurity scores were calculated. The weight ω i used for each node is assumed to be equal to the number of learning cases present in that node. Since the importance of a particular node depends on the number of cases in that node, this assumption is valid. All paths with an insecurity score higher than 30 % were considered to be insecure. The rows in gray highlight the terminal nodes whose predicted class is different from that given by the terminal node based method. Details regarding the classification results of DT-D2 on the learning set cases (July 19) and the test cases (July 26) are given below (Table 4.11). The total number of cases in each terminal node along with the predicted and actual class is tabulated. The classification accuracy obtained using the path based method is shown in Table 4.12. Using the new approach the accuracy for the insecure cases increases from 75.9 % to 98.72 % and that for the secure cases decreases from 80.90 % to 74.56 %. Since it is more important to ensure that insecure cases are not predicted as secure the results obtained are promising. The robustness of DT-D2 can be enhanced by using a path based approach rather than the traditional terminal node based classification approach. 47 Table 4.10 Insecurity scores of paths for DT-D2 Path R-1 R-A-C-2 R-A-B-D-3 R-A-B-D-4 R-A-B-E-5 R-A-C-F-6 R-A-B-E-G-I-7 R-A-B-E-G-I-8 R-A-B-E-G-J-9 R-A-B-E-G-J-10 R-A-C-F-H-L-11 R-A-C-F-H-L-12 R-A-C-F-H-K-M-13 R-A-C-F-H-K-M-14 R-A-C-F-H-K-N-15 R-A-C-F-H-K-N-16 Insecurity score 0.00 29.32 39.47 41.69 42.07 33.28 45.99 46.58 44.63 44.06 25.50 24.42 26.80 25.04 26.21 25.03 Predicted class (terminal node based) S S D D S D S D D S D S D S D S Predicted class (path based) S S D D D D D D D D S S S S S S July 26, 2006 (Secure/Insecure) 3334/2 296/22 225/13 56/222 302/38 506/544 54/193 8/23 277/1 0/0 0/0 0/0 0/0 0/0 0/0 556/0 48 Table 4.11 Terminal node details for DT-D2 Node 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 Predicted class S S D D S D S D D S D S D S D S Learning Set (July 19, 2006) Total number of S/I (actual cases in node class) 6932 6930/2 253 252/1 199 64/135 226 26/200 552 548/4 696 471/1225 192 192/0 26 7/19 220 126/94 192 192/0 34 7/27 474 474/0 430 171/259 359 359/0 90 55/35 669 655/14 Test Set (July 26, 2006) Total number of S/I (actual cases in node class) 3336 3334/2 318 296/22 238 225/13 278 56/222 340 302/38 1050 506/544 247 54/193 31 8/23 278 277/1 0 0/0 0 0/0 0 0/0 0 0/0 0 0/0 0 0/0 556 556/0 Table 4.12 Prediction accuracy for July 26 cases using path based method (DT-D2) Actual Class D S Total Cases 1,058 5,614 Percent D S Correct N=1875 N=4797 98.72 1,851 24 74.56 1,428 4,186 49 4.5 Conclusions This report proposes an online security assessment scheme based on PMUs and DTs. The scheme was applied to an important operational area of the Entergy system. Off line simulations show that the proposed approach can build DTs with high prediction accuracies to reliably identify stability problems for prospective operating conditions and provide guidelines for preventive control. The scheme also identifies critical security indicators, which could be candidates for new PMU locations. A new paths-based DT classification method is also proposed and compared with the traditional terminal nodesbased method to regulate the DTs prediction reliability for changes in operating conditions. This new method focuses on the development of self-adaptive rules to make better use of available DTs and introduces a shift from “Black-Box” DTs to “White-Box” DTs, or rule-based DT classification systems. The DTs built in this report can reach prediction accuracies of about 97% for insecure cases. This research is a first step towards the goal that in the future, more practical versions of the scheme proposed can achieve a high, accuracy (e.g. 99% or even 99.9%) acceptable by electricity industry. To realize that goal, several issues will have to be carefully analyzed. One aspect that is currently being investigated is the impact of the use of a larger number of operating conditions to increase the number of insecure cases used in DT training. 50 References [1] K.-S. Cho, J.-R. Shin, and S. H. Hyun, "Optimal placement of phasor measurement units with GPS receiver," IEEE Power Engineering Society Winter Meeting, vol. 1, 2001, pp. 258-262, 2001. A. G. Phadke, "Synchronized phasor measurements-a historical overview," IEEE/PES Transmission and Distribution Conference and Exhibition: Asia Pacific, vol. 1, pp. 476-479, 2002. W. Lewandowski, J. Azoubib, and W. J. Klepczynski, "GPS: primary tool for time transfer," Proceedings of the IEEE, vol. 87, pp. 163-172, 1999. X. Dongjie, H. Renmu, W. Peng, and X. Tao, "Comparison of several PMU placement algorithms for state estimation," Eighth IEE International Conference on Developments in Power System Protection, vol. 1, pp. 32-35, 2004. A. G. Phadke, B. Pickett, M. Adamiak, M. Begovic, G. Benmouyal, R. O. Burnett, Jr., T. W. Cease, J. Goossens, D. J. Hansen, M. Kezunovic, L. L. Mankoff, P. G. McLaren, G. Michel, R. J. Murphy, J. Nordstrom, M. S. Sachdev, H. S. Smith, J. S. Thorp, M. Trotignon, T. C. Wang, and M. A. Xavier, "Synchronized sampling and phasor measurements for relaying and control," IEEE Transactions on Power Delivery, vol. 9, pp. 442-452, 1994. R. O. Burnett, Jr., M. M. Butts, and P. S. Sterlina, "Power system applications for phasor measurement units," Computer Applications in Power, IEEE, vol. 7, pp. 813, 1994. S. Zhong and A. Abur, "Combined state estimation and measurement calibration," IEEE Transactions on Power Systems, vol. 20, pp. 458-465, 2005. H. Bai, S. Zhou, and Z. Guo, "Innovation network graph state estimation based PMUs," IEEE/PES Transmission and Distribution Conference and Exhibition: Asia and Pacific, pp. 1-6, 2005. J. Chen and A. Abur, "Improved bad data processing via strategic placement of PMUs," IEEE Power Engineering Society General Meeting, vol. 1 pp. 509-513, 2005. D. Junce and C. Zexiang, "Mixed Measurements State Estimation Based on Wide-Area Measurement System and Analysis," IEEE/PES Transmission and Distribution Conference and Exhibition: Asia and Pacific, pp. 1-5, 2005. C. Rakpenthai, S. Premrudeepreechacharn, S. Uatrongjit, and N. R. Watson, "An Improved PMUs Placement Method for Power System State Estimation," The 7th International Power Engineering Conference, IPEC, pp. 1-4, 2005. [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] 51 [12] L. Zhao and A. Abur, "Multi area state estimation using synchronized phasor measurements," IEEE Transactions on Power Systems, vol. 20, pp. 611-617, 2005. E. Price, "Practical considerations for implementing wide area monitoring, protection and control," 59th Annual Conference for Protective Relay Engineers, pp. 12, 2006. M. Rice and G. T. Heydt, “Power Systems State Estimation Accuracy Enhancement Through the Use of PMU Measurements,” IEEE PES Transmission and Distribution Conference and Exposition, Dallas, TX, May 2006. M. Rice and G. T. Heydt, Phasor Measurement Unit Data in Power System State Estimation, PSerc Intermediate Project Report for "Enhanced State Estimators," April 2005. M. Rice, Phasor Measurement Unit Data in Power System State Estimation, Masters Thesis, Dec. 2004. G. B. Denegri, M. Invernizzi, and F. Milano, "A security oriented approach to PMU positioning for advanced monitoring of a transmission grid," Proceedings of International Conference on Power System Technology, vol. 2, pp. 798-803, 2002. R. Zivanovic and C. Cairns, "Implementation of PMU technology in state estimation: an overview," IEEE AFRICON 4th, vol. 2, pp. 1006-1011, 1996. I. Kamwa and R. Grondin, "PMU configuration for system dynamic performance measurement in large, multiarea power systems," IEEE Transactions on Power Systems, vol. 17, pp. 385-394, 2002. S. B. M. Ingram, S. Matthews, A. P. Meliopoulos, G Cokkinides, "Use of phasor measurements, SCADA and IED data to improve the state estimation," 7th Fault and Disturbance Analysis Conference, April 26-27, 2004. A. P. S. Meliopoulos and G. K. Stefopoulos, "Characterization of state estimation biases," International Conference on Probabilistic Methods Applied to Power Systems, pp. 600-607, 2004. S. Widergren, "Eastern Interconnect Phasor Project." http://phasors.pnl.gov/ "IEEE standard for synchrophasor for power systems. IEEE std. 1344-1995," "IEEE Standard for Synchrophasors for Power Systems," IEEE Std C37.118-2005 (Revision of IEEE Std 1344-1995), pp. 1-57, 2006. [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24] 52 [25] W. Hongxia, N. Hui, and G. T. Heydt, "The impact of time delay on robust control design in power systems," IEEE Power Engineering Winter Meeting, vol.2, pp. 1511-1516, 2002. J. Bertsch, M. Zima, A. Suranyi, C. Carnal, and C. Rehtanz, "Experiences with and perspectives of the system for wide area monitoring of power systems," CIGRE/IEEE PES International Symposium Quality and Security of Electric Power Delivery Systems, pp. 5-9, 2003. D. Karlsson, M. Hemmingsson, and S. Lindahl, "Wide area system monitoring and control - terminology, phenomena, and solution implementation strategies," IEEE Power and Energy Magazine, vol. 2, pp. 68-76, 2004. J. F. Hauer, N. B. Bhatt, K. Shah, and S. Kolluri, "Performance of "WAMS East" in providing dynamic information for the North East blackout of August 14, 2003," IEEE Power Engineering Society General Meeting, vol. 2, pp. 1685-1690, 2004. X. Tong, G. Liao, X. Wang, and S. Zhong, "The analysis of communication architecture and control mode of wide area power systems control," Proceedings Autonomous Decentralized Systems, pp. 59-65, 2005. K. E. Holbert, G. I. Heydt, and H. Ni, "Use of satellite technologies for power system measurements, command, and control," Proceedings of the IEEE, vol. 93, pp. 947-955, 2005. I. Kamwa, R. Grondin, and Y. Hebert, "Wide-area measurement based stabilizing control of large power systems-a decentralized/hierarchical approach," IEEE Transactions on Power Systems, vol. 16, pp. 136-153, K. Tomsovic, D. E. Bakken, V. Venkatasubramanian, and A. Bose, "Designing the next generation of real-time control, communication, and computations for large power systems," Proceedings of the IEEE, vol. 93, pp. 965-979, 2005. P. Kundur, Power System Stability and Control, McGraw-Hill, 1994. M. Pavella, D. Ernst, D. Ruiz-Vega, Transient Stability of Power Systems: A Unified Approach to Assessment and Control, Kluwer Academic Publishers, 2000. L. Wehenkel, Automatic Learning Techniques in Power Systems, Kluwer Academic Publishers, 1998. I. Houben, L. Wehenkel, M. Pavella, “Coupling of K-NN with Decision Trees for Power System Transient Stability Assessment,” Proceedings of the 4th IEEE Conference on Control Applications, pp. 825-832, September 1995. [26] [27] [28] [29] [30] [31] [32] [33] [34] [35] [36] 53 [37] L. Wehenkel, Th. Van Cutsem, M. Pavella, “An Artificial Intelligence Framework for Online Transient Stability Assessment of Power Systems,” IEEE Trans. Power Syst., vol. 4, no. 2, pp. 789-800, May 1989. A. A. Fouad, Tong Jianzhong, “Stability Constrained Optimal Rescheduling of Generation,” IEEE Trans. Power Syst., vol. 8, no. 1, pp. 105-112, February 1993. E. S. Karapidakis, N. D. Hatziargyriou, “Online Preventive Dynamic Security of Isolated Power Systems Using Decision Trees,” IEEE Trans. Power Syst., vol. 17, no. 2, pp. 297-204, May 2002. L. Wehenkel, Th. Van Cutsem, M. Pavella, “Inductive Inference Applied to Online Transient Stability Assessment of Electric Power Systems,” Automatica, vol. 25, pp. 445-451, May 1989. L. Wehenkel, M. Pavella, E. Euxibie, B. Heilbronn, “Decision Tree Based Transient Stability Method A Case Study,” IEEE Trans. Power Syst., vol. 9, no. 1, pp. 459-469, February 1994. L. Wehenkel, M. Pavella, “Decision Trees and Transient Stability of Electric Power Systems,” Automatica, vol. 27, no. 1, pp. 115-134, January 1991. L. Wehenkel, M. Pavella, “Advances in Decision Trees Applied to Power System Security Assessment,” Proc. IEE 2nd International Conference on Advances in Power System Control, Operation and Management, pp. 47-53, December 1993. S. Rovnyak, S. Kretsinger, J. Thorp, D. Brown, “Decision Trees for Real-time Transient Stability Prediction,” IEEE Trans. Power Syst., vol. 9, no. 3, pp. 14171426, August 1994. K. Mei, S. Rovnyak, “Response-Based Decision Trees to Trigger One-Shot Stabilizing Control,” IEEE Trans. Power Syst., vol. 19, no. 1, pp. 531-537, February 2004. N. Senroy, G. T. Heydt, V. Vittal, “Decision Tree Assisted Controlled Islanding,” IEEE Trans. Power Syst., vol. 21, no. 4, pp. 1790-1797, November 2006. L. Breiman, J. Friedman, R. A. Olshen, C. J. Stone, Classification and Regression Trees, Wadsworth, 1984. “TSAT Transient Security Assessment Tool,” [Online] (Date Accessed: April 2007), Available: http://www.dsapowertools.com/downloads/TSAT_Brochure.pdf CART® Tree-Structured Non-Parametric Data Analysis. [Online]. Available: http://www.salfordsystems.com/cart.php [38] [39] [40] [41] [42] [43] [44] [45] [46] [47] [48] [49] 54 Project Publications K. Sun, S. Likhate, V. Vittal, V. Kolluri, and S. Mandal, “An Online Dynamic Security Assessment Scheme using Phasor Measurements and Decision Trees,” Accepted for publication in the IEEE Transactions on Power Systems. 55 A Tool for On-line Stability Determination and Control for Coordinated Operations between Regional Entities Using PMUs Volume 2 ` Project Team Vijay Vittal, Project Leader, Arizona State University Gerald T. Heydt, Arizona State University A.P. Sakis Meliopoulos, Georgia Institute of Technology Information about this project For information contact: Vijay Vittal Ira A. Fulton Chair Professor Department of Electrical Engineering Arizona State University PO Box 875706 Tempe, AZ 85287-5706 Tel: 480-965-1879 Fax: 480-965-0745 Email: [email protected] Power Systems Engineering Research Center This is a project report from the Power Systems Engineering Research Center (PSERC). PSERC is a multi-university Center conducting research on challenges facing a restructuring electric power industry and educating the next generation of power engineers. More information about PSERC can be found at the Center’s website: http://www.pserc.org. For additional information contact: Power Systems Engineering Research Center Arizona State University 577 Engineering Research Center Box 878606 Tempe, AZ 85287-8606 Phone: 480-965-1643 FAX: 480-965-0745 Notice concerning copyright material PSERC members are given permission to copy without fee all or part of this publication for internal use if appropriate attribution is given to this document as the source material. This report is available for downloading from the PSERC website. © 2007 Arizona State University and Georgia Institute of Technology. All rights reserved. Table of Contents 1 Introduction................................................................................................................... 1 1.1 Background........................................................................................................... 1 1.2 Deployment of GPS Synchronized Equipment .................................................... 6 1.3 Characterization of GPS Synchronized Measurements........................................ 7 2 SuperCalibrator Concept............................................................................................. 11 2.1 Introduction ........................................................................................................ 11 2.2 Description of the SuperCalibrator Concept ...................................................... 11 2.3 Dynamic State Estimation .................................................................................. 14 2.3.1 State Vector ................................................................................................. 15 2.3.2 System Model.............................................................................................. 16 2.3.3 Measurement Set ......................................................................................... 18 2.3.4 Dynamic State Estimator (DSE) Algorithm................................................ 22 3 Stability Monitoring.................................................................................................... 24 3.1 Introduction ........................................................................................................ 24 3.2 Introductory Example ......................................................................................... 25 3.3 Generator “Potential” Energy............................................................................. 26 3.4 Generator “Kinetic” Energy ............................................................................... 27 3.5 Generator “Total” Energy................................................................................... 27 4 Case Study .................................................................................................................. 30 4.1 Demonstration Example ..................................................................................... 30 5 References................................................................................................................... 37 5.1 Papers Resulted from this Project....................................................................... 37 5.2 Background Papers............................................................................................. 37 Appendix A. Quadratic Integration Method ..................................................................... 40 A.1 Introduction ........................................................................................................ 40 A.2 Description of Quadratic Integration Method .................................................... 41 A.3 Numerical properties .......................................................................................... 42 Appendix B. Generator Model.......................................................................................... 45 B.1 Introduction ........................................................................................................ 45 B.2 Full Time Domain Model................................................................................... 46 B.3 Quasi-Steady State Model (Mixed Domain) ...................................................... 49 B.4 Frequency Domain Steady State Model ............................................................. 54 i List of Figures Figure 1-1: Evolution of Time Pieces................................................................................. 2 Figure 1-2: Block Diagram of Arun Phadke’s PMS (taken from Phadke) ......................... 3 Figure 1-3: The Macrodyne 1620 PMU (January 1992) .................................................... 5 Figure 1-4: Geographical Location of Some PMUs in the Eastern Interconnection Phasor Project ................................................................................................................................. 7 Figure 1-5: Typical potential and current instrumentation channels. ................................. 8 Figure 2-1. Functional Description of the Distributed State Estimator ............................ 13 Figure 2-2. Conceptual Illustration of the SuperCalibrator concept................................. 13 Figure 2-3: Schematic diagram of dynamic state estimator.............................................. 15 Figure 2-4: Pseudo-measurements of Voltage at Other End of Lines .............................. 20 Figure 2-5: Pseudo-measurements from Kirchoff’s Current Law .................................... 21 Figure 2-6: Pseudo-measurements from Neutral/Shield Wire Current............................. 22 ii 1 Introduction The research effort described in this report explores the capabilities of GPS synchronized measurement for the purpose of developing better tools for monitoring and controlling the stability of an electric power system. It is therefore pertinent to review the GPS synchronized devices and their capabilities as relate to the objectives of this project. We provide an overview of the GPS synchronized measurement technology development, followed by the present deployment activities as well as the characterization of data obtained with GPS synchronized devices. This information is then utilized in subsequent sections to describe the developments under this project. Specifically, the primary objective of this project was to develop applications using synchronized phasor data that are being implemented at several electric utility companies in the US. The application that was specifically addressed in this project dealt with online transient stability assessment using the measurements from phasor measurement units (PMUs). PMUs provide direct measurements of phase angles (on a common time reference – GPS clock). It is well known that the numerical values of phase angles alone do not always provide information regarding the stability of the system. We have investigated methodologies for extracting stability information from PMU data. These methodologies are described in subsequent sections of the report. 1.1 Background It is important to recognize that simply time tagging measurements with a GPS clock is not enough to provide the necessary accuracy for the applications described in this report. While standards for the timing accuracy of GPS synchronized measurements does not exist, in this report we will assume that GPS synchronized measurements mean that the measurements are known with a time accuracy of 1 to 2 microseconds. This accuracy translates into 0.0216 to 0.0432 degrees at 60 Hz. While the technology exists today to achieve this timing accuracy, PMU performance may vary from manufacturer to manufacturer. To better understand the technology, it is expedient to review the need and the historical developments towards phasor measurements. The value of synchronized measurements has been recognized in many applications. As a matter of fact, synchronization of measurements and observations was sought in ancient times. Key to this capability is the development of accurate clocks. In ancient times the “klepsidra” was typically used for timing. Recently it has been discovered that the ancient Greeks developed a clock mechanism (the antikythira mechanism) that is illustrated in Figure 1-1. This mechanism was the basis of the mechanical watches. As technology evolved, the accuracy of clocks has increased. Today the most accurate clock is the Cesium atomic clock with accuracy of 1 microsecond in 100 years. The cost of the Cesium atomic clock is prohibitive for many applications. However, the Federal government has developed and deployed a constellation of satellites that transmit the signal of the Cesium clock throughout the globe and if the signal is received from more than three satellites, with triangulation/estimation methods the signal of the Cesium atomic clock can be reproduced at any point on earth with accuracy of tens of 1 nanoseconds (from the power system applications point of view a timing signal with accuracy better than one microsecond is enough). The Antikythira Mechanism GPS Satellite Figure 1-1: Evolution of Time Pieces In 1983, the first author (the third author) proposed to EPRI/Utah Power company to develop a system of synchronized measurements using timing signals from GOES satellites for the purpose of monitoring the power loop flow through the state of Utah and the swings of this power flow in response of concern by this utility for the frequent loop flow problems that they were experiencing. The GOES signal has an accuracy of better than 100 microseconds and at that time it was the best available clock (other than a Cesium atomic clock). While this proposal was not funded, the authors proposed to NYPA in 1989 to assess and develop technology of synchronized measurements for the purpose of monitoring harmonic distortion on the NYPA transmission system. This work was funded and in 1991 the authors developed a prototype that used signals from a GPS clock and a time vernier method to time tag measurements with accuracy better than 2 microseconds [24]. This device could be inserted in the input of commercially available digital fault recorder and it would time tag the measurements with accuracy better than 2 microseconds. In the period of 1990-1992, Arun Phadke developed the PMS which is illustrated in Figure 1-2 [29]. The PMS used a GPS signal for timing, 720 samples per second sample and hold A/D converter and a front-end anti-aliasing filter with a cutoff frequency of 360 Hz. The combination of the anti-aliasing filter and the multiplexing introduce time delays that are orders of magnitude greater than the accuracy of the GPS clock. Although this device was never tested by independent organizations, the estimated timing errors are more than 50 microseconds. Several PMS were constructed and sold to several utilities 2 (for example three to AEP, one to NYPA, etc.). Despite the use of the GPS clock, the PMS was not capable of performing measurements with comparable accuracy to the GPS clock. Figure 1-2: Block Diagram of Arun Phadke’s PMS (taken from Phadke) (Characteristics: (a) analog antialiasing input filter with a cutoff frequency of 360 Hz, (b) 12 bit sample and hold A/D technology 720 samples per second with analog multiplexing) The first device capable of performing synchronized measurements with accuracy comparable to the GPS clock accuracy was developed by Jay Murphy of Macrodyne and was released in the market in January 1992 [30]. Jay Murphy named the device the Macrodyne 1620 PMU (Phasor Measurement Unit). Jay Murphy introduced the following innovations to achieve the goal of performing synchronized measurements with accuracy comparable to the GPS clock: • • • Individual Channel GPS Sync Common Mode Rejection Filter with Optical Isolation 16 bit A/D Sigma/Delta Modulation converter The Macrodyne 1620 PMU was tested by the authors who verified that the unit is capable of performing measurements with accuracy of one microsecond (or equivalently 3 0.02 degrees for a 60 Hz waveform). The block diagram and the picture of the Macrodyne 1620 PMU is illustrated in Figure 1-3. 4 Input Protection & Isolation Section GPS Antenna A/D Converter (ΣΔ Modulation) Optical Isolation Digitized Data 2880 s/s GPS Receiver Optical Isolation 1PPS IRIGB Analog Inputs V : 300V I : 2V Input Protection & Isolation Section Sampling Clock PLL A/D Converter (ΣΔ Modulation) Optical Isolation Digitized Data 2880 s/s μP Optical Isolation Display & Keyboard RS232 Data Concentrator (PC) Master Workstation (a) (b) Figure 1-3: The Macrodyne 1620 PMU (January 1992) (a) Block Diagram, (b) Photograph 5 Memory For a long period of time (1992 to 2002), the only GPS synchronized equipment was the Macrodyne PMU. Recently additional GPS synchronized equipments were introduced in the market. Yet, standards that determine what the accuracy of the phase measurement should be do not exist. We argue that the accuracy of the phasor measurements should be such that the error in predicting the power flow should not exceed 1%. If we consider a 50 mile long 230 kV line, rated 400 MVA and evaluate the precision in voltage magnitude and phase angle measurements required to achieve a 1% accuracy in power flow then we have the following pairs: • • • • Voltage Magnitude: 0.5%, Phase Angle: 0 degrees Voltage magnitude: 0.4%, Phase Angle: 0.03 degrees Voltage magnitude: 0.3%, Phase Angle: 0.05 degrees Voltage magnitude: 0.2%, Phase Angle: 0.09 degrees Similar analysis could lead to a desired standard. Ignoring other sources of error, most GPS synchronized devices have the capability to measure voltage magnitude with precision 0.1% and the phase angle with precision 0.02 degrees. In this case the expected error in the power flow for the above mentioned line will be: 0.34%. Unfortunately, this precision cannot be achieved because of the other errors that have been mentioned. A year ago we reported on a new concept for correcting for these errors. We introduced the concept of the SuperCalibrator [31, 32]. The SuperCalibrator concept provides a practical approach for correcting the errors arising from instrumentation channels in a practical application. In this paper we propose the use of the SuperCalibrator as a practical and highly precise distributed state estimator. Specifically, the SuperCalibrator is applied at the substation level to provide the substation state which is globally valid if there is a valid GPS-synchronized measurement in the data set. In this report we describe (in later sections) one particular application of the SuperCalibrator for extracting the dynamic state of the power system. 1.2 Deployment of GPS Synchronized Equipment The perceived benefits of GPS synchronized measurements for the electric power system are many including direct measurement of the system state and time alignment of disturbance data for post mortem analysis. In 2003 an initiative was propelled by the August 2003 blackout to utilize PMUs for monitoring the Eastern Interconnection, the EIPP (Eastern Interconnection Phasor Project). Figure 1-4 illustrates some of the PMUs used in the system at the time of the start of this project. The end result of this project is the availability of real time PMU data at a central location. This central location is maintained by TVA and procedures are in place for utilization of the PMU data from all the stakeholders. 6 E ID Mas sena E ID E ID N iagaraobins on R oad R E ID E ID MA R C Y 765 MA R C Y 345 E ID FR A S E R E ID Orange, A E P E ID C allaw ay, A meren E ID R us h Island, A meren E ID K anaw ha R iv er E ID R oc k port, A E P E ID E ID P aradise Foss il P lant S haw nee Fos sil P lant E ID Marshall E ID Jack s ons Ferry, 765 kV E ID Montgomery E ID E ID E ID C umberland Fos sil P lant E ID E ID E ID S ullivan, 500 kV W eak ley W ils on, TN E ID DIEavidson H ook, TN D P in Johns onv ille Foss il P lant Maury E ID E ID V olunteer B ull R un Fos sil R oane E ID E ID E ID J ack son, TN H ayw ood, TN S W S TA E ID S helby C ordov a E ID R oc k S pring Metering B ellefonte N uclear E ID E ID Freeport E ID Madis on E ID E ID P leasant H ill, MS U nion E ID E ID Trico S teel Trinity E ID S heridan, A R E ID E ID Limes tone, A L E ID B ow en E ID E l D orado, A R E ID E IDID E Low ndes, MS W es t P oint Farley E ID Miller E ID Frank lin, MS E ID H artburg, TX E ID W aterford Figure 1-4: Geographical Location of Some PMUs in the Eastern Interconnection Phasor Project 1.3 Characterization of GPS Synchronized Measurements An important issue is the knowledge of the data accuracy. Normally the GPS data represent the power level voltages and currents that are obtained by first transforming the power level voltages and currents to instrumentation level and then the GPS synchronized equipment digitize the reduced level voltages and currents. Assuming an ideal transfer function of the overall instrumentation channel, the power level voltages and currents are obtained. Unfortunately, the instrumentation channel does not have ideal characteristics. We describe typical instrumentation channels and typically expected overall errors from these channels for usual installations of GPS synchronized equipment. 7 The chain of measurement starts from the high voltage or current measurement point and it ends at the digital signal generated by the A/D converters of a GPS synchronized device. The components/devices in between are referred to as the instrumentation channel. Figure 1-5 illustrates the devices forming voltage and current channels typically found in electric power generating stations and substations. Figure 1-5: Typical potential and current instrumentation channels. The first link in the instrumentation channel equipment chain consists of voltage and current transformers, collectively called instrument transformers. These devices transform power system voltages and currents to levels appropriate for driving relays, fault recorders, PMUs and other monitoring equipment. Several instrument transformer technologies are presently in use. The most common traditional technology devices are voltage and current transformers (VTs and CTs), which are based on magnetic core transformer technology. Another type of commonly used voltage transducers are capacitively coupled voltage transformers (CCVT’s). These are based on a combination of capacitive voltage dividers and magnetic core transformers. Recently, voltage and current instrument transformers have been constructed based on the electro-optical and magneto-optical phenomena. These devices are known as EOVT’s (Electro-Optical Voltage Transformers) and MOCT’s (Magneto-Optical Current Transformers). High voltage instrumentation channels introduce errors to phasor measurements. The level of error is dependent upon the type of instrument transformers, control cable type and length and protection circuitry at the input of the A/D converters. As an example Figure 1-6: illustrates the errors for a specific instrumentation channel. Note that the VT introduces a very small error (0.01 degrees), while the 500 ft cable introduces an error of 0.4 degrees. The overall error is more than an order of magnitude higher than the error of a typical PMU. 8 Voltage Measurement IC Substation A, 115 kV Bus Van = 62.53 kV / 27.52 Deg Van = 62.19 V / 27.51 Deg Vbn = 62.96 kV / -92.68 Deg Vbn = 62.61 V / -92.70 Deg Vcn = 62.33 kV / 147.46 Deg Vcn = 61.99 V / 147.45 Deg NORTHBUS3 NBUS3MS Van = 61.63 V / 27.11 Deg Vbn = 63.09 V / -92.85 Deg Vcn = 61.72 V / 148.00 Deg NBUS3MSI RG-8 Cable, 500 ft 69kV:69V Wound Type VT Figure 1-6: Illustration of Instrumentation Channel Errors for a Typical Case Depending on the instrumentation channel, the characterization of these errors may be substantial. It should be recognized that GPS-synchronized equipment may be also connected to existing instrumentation in substations that may be for other purposes, i.e. metering. Many times the instrument transformers are connected in an arrangement that generates a phase shift, for example delta connection. The resulting phase shift must be accounted for. Optical CTs and VTs present a special case. These devices are quite accurate. However, the present designs are equipped with an analog output for compliance with standards. The analog output introduces time latencies. Typical time latencies are in the order of few tens of microseconds. This is orders of magnitude greater than the time errors from typical PMUs. However, there is no need to use the analog output. The data directly from the optical CTs and VTs are digital and should be used directly without the analog output. This will eliminate the time latencies but it will require developing new standards or adopting existing standards for data communications. An evaluation of typical arrangements of instrumentation channels feeding PMUs reveals the following conclusions and comments. GPS-synchronized equipment is very accurate as compared to standard power system instrumentation. However application of this equipment to a practical power system is burdened by the errors that are introduced by the standard instrumentation channels utilized in power systems. The inaccuracies introduced by the instrumentation channels do not present a problem for many applications. However for the application targeted in this project, i.e. stability monitoring, these inaccuracies are important and should be corrected. 9 The most accurate instrumentation channels are current instrumentation channels that use CTs. The length of the control cable is very important in determining the level of errors. The next most accurate instrumentation channels are voltage instrumentation channels that use wound type VTs. The length of the control cable is also very important in determining the level of errors. CCVT based instrumentation channels are relatively accurate when they are well calibrated. They perform well when the frequency is near nominal. A main drawback is that the parameters of the components shift with time and subsequently introduce large errors. In addition, during transients the error is very large because their characteristics deteriorate at frequencies other than the fundamental. On the other hand, due to economic factors, CCVTs are commonly used in high voltage applications. The optical VTs and CTs are high accuracy devices for magnitude measurement but very poor for phase angle measurement. Specifically, they exhibit a time delay in the order of few tens of microseconds that translates to phase angle error in the degree range. The data directly from the optical CTs and VTs are digital and should be used directly. This will eliminate the time latencies but it will require developing new standards or adopting existing standards for data communications. In summary, PMUs are relatively accurate devices. Yet in a typical application the overall accuracy of the PMU data is typically compromised because of the errors introduced by the non-ideal characteristics of the instrumentation channel. For the application at hand, the errors form the instrumentation channels should be corrected. 10 2 SuperCalibrator Concept 2.1 Introduction Recent major outages such as the 2003 U.S.-Canada Black-Out emerged the need for more efficient methods of monitoring the state of power systems and in particular the stability of the system. Traditional approaches for state estimation in modern Energy Management Systems do not provide information about the stability status of the system. Very few utilities have implemented dynamic security assessment that may provide information about critical contingencies that may trigger system instabilities. Dynamic security assessment operates on the present operating condition by performing stability analysis on a number of contingencies. Presently there is no methodology to monitor in real time the stability status of the electric power system. The introduction of GPS synchronized measurements makes it possible to estimate the dynamic state of the system and consequently to monitor the stability of the system. In addition it is possible to decentralize the state estimation process. An approach that will make this task possible has been developed and described in this report. The approach has been named SuperCalibrator because it enables fine filtering and calibration of the PMU data. The SuperCalibrator consists of a filter (state estimation) of the voltage and current data using a three-phase, breaker oriented, instrumentation inclusive model. The measurements obtained from SCADA, relays, PMUs, meters etc are used in a statistical estimation method that fits the data to the detailed model while identification of bad data is also implemented. We describe the basic approach of the SuperCalibrator and then present its application to the dynamic state estimation and stability monitoring. 2.2 Description of the SuperCalibrator Concept The SuperCalibrator is conceptually very simple. The technology is based on a flexible hybrid state estimation formulation [33]. This is a combination of the traditional state estimation formulation and the GPS-synchronized measurement formulation, which uses an augmented set of available data. The basic idea is to provide a model based correction of the errors from all known sources of errors. Specifically, the basic idea is to utilize a detailed model of the substation, (three-phase, breaker oriented model, instrumentation channel inclusive and data acquisition model inclusive). The set of measurements comprises of: a. Traditional, non-synchronized measurements (voltage magnitude, active and reactive line flows and bus injections, and other standard SCADA data), b. GPS-synchronized measurements of voltage and current phasors for each phase, c. GPS synchronized measurements of frequency and rate of frequency change, d. Appropriate pseudo-measurements as they will be explained in the subsequent sections. 11 Typical measurements are illustrated in Table 2-1. Then the measurements obtained with any device, PMU, relay, SCADA etc is utilized in an estimation method that statistically fits the data to the detailed model. Table 2-1: List of Measurements Phasor Measurements Description ~ Voltage Phasor, V ~ Current Phasor, I ~ Current Injection Phasor, I inj Frequency, f Rate of frequency change, df/dt Other Non-synchronized Measurements Description Voltage Magnitude, V Real Power Flow, P f Reactive Power Flow, Q f Real Power Injection, Pinj Reactive Power Injection, Q inj Other This approach leads to a distributed state estimation procedure performed at the substation level. The functional description of the distributed state estimation is illustrated in Figure 2-1 and consists of the following procedures: 1. perform state estimation on each substation using all available data and a three-phase, breaker oriented, instrumentation inclusive model, 2. perform bad data identification and correction (or rejection) as well as topology error identification on each substation, 3. collect the results from all substations at a central location (control center) to construct the system-wide operating state of the system. The advantages of the SuperCalibrator illustrated in Figure 2-2 are: a. Utilization of a high fidelity model (three-phase, breaker-oriented and instrumentation channel inclusive model), b. Utilization of all available data at the substation level, i.e. SCADA, relays, meters, PMUs, DFRs, etc., c. Metrics of state estimation accuracy based on the well known chi-square test and d. Distributed processing of data at each substation (this requires that there is at least one GPS synchronized measurement at each substation) and subsequent transfer of the substation estimation results to the control center for formation of the system wide state estimate. The above SuperCalibrator provides the static state of the system. The approach has been extended to provide the dynamic state of the system or to include a dynamic state estimator. In subsequent paragraphs we describe the dynamic state estimator that provides the dynamic state at each bus of the system. It is assumed that GPS synchronized measurements exist in the substation and all incoming lines. 12 PM Data U SCADA Data Relay Data O ther Substation Level State Estim ation Topology Error Detection and Correction Alarm Processing (Root Cause Analysis) State Estim ation Coordination Positive Sequence State Esim e at Three Phase State Estim ate Expected Error vs Confidence Identified Bad Data Root Cause Event of Alarm s Breaker O riented, Instrum entation Inclusive Three Phase Substation M odel Figure 2-1. Functional Description of the Distributed State Estimator Measurement Layer Phase Conductor v(t) i(t) Current Transformer Potential Transformer i1(t) i2(t) IED Vendor D Relay Vendor C Data Processing Attenuator Burden PMU Vendor A Substation PDC Instrumentation Cables v1(t) v2(t) Attenuator Anti-Aliasing Filters PMU Vendor C Proprietary Protocol Burden Figure 2-2. Conceptual Illustration of the SuperCalibrator 13 2.3 Dynamic State Estimation Dynamic state estimator extends the concept of the static state estimation by using the dynamic states of the power system such as the generator speed (or equivalently the frequency of the voltage) or the generator acceleration (or equivalently the rate of frequency change). For the purpose of generalization, we define the state of the system as the collection of the voltage phasors at each bus of the system and the frequency of the voltage phasor at each bus. Since we apply the dynamic state estimator on the three phase model of the system, then each bus is characterized with the four variables: the three voltage phasors (one for each phase) and the average frequency of these voltages. Note that we make the assumption that the voltages of each phase at a bus will have the same frequency – this assumption is reasonable. The dynamic state estimator is applied to the integrated dynamic model of the power system. The model is based on a detailed three-phase, breaker-oriented and instrumentation inclusive dynamic representation of the substation configuration, the generating units and the interconnected transmission lines. In general, the model is described by a set of algebraic and differential equations, which express the dynamics of the system and is described as follows: dx(t ) = f ( x(t ), y (t ), t ) , dt 0 = g ( x(t ), y(t ), t ) where x and y are the dynamic and algebraic states of the system respectively. The system measurements z are related to the states through the following equation: z = h( x(t ), y(t ), t ) + η , where η represents the measurement error which is often described as random variables with Gaussian probability distribution. The measurement error is directly related to the accuracy of the data acquisition system. Note that variable z can include measurements such as frequency. The overall approach is illustrated in Figure 2-3. 14 zα zm e ˆˆˆ xe, ye , ze ˆ xe Figure 2-3: Schematic diagram of dynamic state estimator. The set of physical measurements za are obtained from PMUs, relays, and other metering devices at the substation level. These values are then compared with the “model” values which are obtained from the dynamic model of the power system, i.e., zm, and with the pseudo-measurements, forming a measurement error e. A standard least square estimation is performed to minimize the sum of the errors squared. The result of the minimization is filtered by a procedure that identifies bad data. If the dynamic state estimator results are not satisfactory, this indicates the presence of bad data. The bad data are identified, removed from the measurement set and the procedure is repeated. Otherwise, the best estimates of the system states and measurements denoted by xe, ye, and ze are calculated. In addition, the estimated values of the dynamic states of the system xe are fed back to the mathematical model to be used as the previous states for the differential equations. These estimated states can also be utilized as inputs to various computational applications such as stability monitoring and prediction. The detailed description of the dynamic state estimator is provided in the following subsections. 2.3.1 State Vector The state of the system for the dynamic state estimator is defined in terms of the voltages and the frequency at the substation as well as the same information at the 15 neighboring substations. The state vector represents the voltage phasors of the three phases and the neutral if applicable: ~ ~ V k , A = V k , A = V k , A, r + jVk , A,i ~ ~ Vk , B = Vk , B = Vk ,B ,r + jVk , B ,i ~ ~ Vk ,C = Vk ,C = Vk ,C ,r + jVk ,C ,i ~ %%% f k (Vk ) : Average frequency of Vk , A ,Vk , A ,Vk ,C 2.3.2 System Model The model of the system is defined in terms of the model of each individual component within the system of interest. The system of interest comprises a substation and all the transmission circuits up to the next substations. Each component of the electric power system is represented with an appropriate model. For example a generating unit is represented with the following model: ζ (t ) = A1v(t ) + A4 z (t ) 0= dx(t ) + B2 x(t ) + B3 y (t ) + B4 z (t ) dt 0 = C1v (t ) + C 2 x (t ) + C 3 y (t ) + C 4 z (t ) 0 = q(v(t ), x(t ), y (t ), z (t )) Assuming that the functions vary quadratically within the integration step [t-h,t], and upon quadratic integration the dynamic model of the generating unit is converted into the following algebraic form (more details on quadratic integration is presented in Appendix A): ζ ( t ) = A1v ( t ) + A4 y ( t ) ζ (t m ) = A1v(t m ) + A4 z(t m ) h 2h h 0 = x (t ) − x (t − h ) + B2 ( x ( t − h ) + x (tm ) + x (t )) + 6 3 6 h 2h h h 2h h B3 ( y (t − h ) + y (tm ) + y (t )) + B4 ( z (t − h ) + z (tm ) + z (t )) 6 3 6 6 3 6 16 0 = x (t ) − x(t − h ) + B2 ( B3 ( −h 24 y (t − h ) + h 3 −h 24 x (t − h) + 5h 24 h 3 x ( tm ) + −h 24 5h 24 x (t )) + h 3 z ( tm ) + 5h 24 z (t )) y ( tm ) + y (t )) + B4 ( z (t − h ) + 0 = C1v(t ) + C2 x(t ) + C3 y (t ) + C4 z (t ) 0 = C1v(t m ) + C 2 x(t m ) + C3 y(t m ) + C 4 z (t m ) 0 = q(v(t ), x(t ), y (t ), z (t )) 0 = q(v(t m ), x(t m ), y(t m ), z(t m )) In above formulation, one of the x variables is frequency and the difference x(t)-x(t-h) represents rate of frequency change. This model is described in Appendix B in greater detail. The dynamical equations represent the mechanical (slow) parts of the system (rotor equations). The electrical transients are neglected. The circuit models such as transmission lines and transformers can be represented by the following expression: ~ ⎡ I S ⎤ ⎡Y11 ⎢~ ⎥ = ⎢ ⎣ I R ⎦ ⎣Y21 ζ (t ) = Ax (t ) ~ Y12 ⎤ ⎡VS ⎤ ⎢~ ⎥ Y22 ⎥ ⎣V R ⎦ ⎦ If ζk,R and ζm,I are the real and imaginary part of phasor ζ respectively, then: ζ = ζ k , R + j ⋅ ζ m, I = a k x(t ) + jam x(t ) ~ ~ where αk, αm are the kth and mth rows of matrix A, respectively. Voltage magnitudes can be computed by the formula: ~= x x 2 k ,R + x 2 m,I while current magnitudes are computed by the formula: 17 ζ = ζ 2 k , R + ζ 2 m, I = (a k x(t )) 2 + (am x(t )) 2 2.3.3 Measurement Set ~ The measurement set is defined based on the available actual measurements in the system and a number of additional pseudo-measurements. 2.3.3.1 Physical Measurements The set of measurements include the voltage phasor measurements at each node (bus), and the current phasors flowing through the lines, transmission lines and so on. Voltage Phasor Measurement: The voltage phasor measurement provides a magnitude and a phase of the corresponding phase voltage. The measurement is expressed as a function of the state of the system: ~~ z = V +η ~ ~ Where V is the corresponding state and η is the measurement error. Frequency Measurement: The frequency measurement provides the frequency of the voltages at a specific bus. The measurement is expressed as a function of the state of the system: ~ z = f +η ~ Where f is the corresponding state (frequency of a specific bus) and η is the measurement error. Electric Current Phasor Measurement: The current phasor measurement provides a magnitude and a phase of the corresponding phase current. An electric current measurement is associated with a circuit. The measurement is expressed as a function of the state of the system utilizing the equation of the circuit. As an example using the circuit model presented earlier an electric current phasor measurement can be expressed as follows: ~ ⎡ I S ⎤ ⎡Y11 z = e ⎢~ ⎥ = ⎢ ⎣ I R ⎦ ⎣Y21 T i ~ Y12 ⎤ ⎡VS ⎤ ~ ⎢ ~ ⎥ +η Y22 ⎥ ⎣VR ⎦ ⎦ where ei is a vector with all entries equal to zero except one that equals 1.0 corresponding ~ to the measured current, and η is the measurement error. 18 Voltage Magnitude Measurement: The voltage magnitude measurement provides a magnitude of the corresponding phase voltage. The measurement is expressed as a function of the state of the system: ~ z = V + η = Vr2 + Vi 2 + η ~ Where V is the corresponding state and η is the measurement error. Electric Current Magnitude Measurement: The current magnitude measurement provides the magnitude of the corresponding phase current. An electric current measurement is associated with a circuit. The measurement is expressed as a function of the state of the system utilizing the equation of the circuit. As an example using the circuit model presented earlier an electric current magnitude measurement can be expressed as follows: ~ ~ ⎡ I ⎤ ⎡Y Y ⎤ ⎡V ⎤ z = eiT ⎢ ~S ⎥ = ⎢ 11 12 ⎥ ⎢ ~S ⎥ + η = I r2 + I i2 + η ⎣ I R ⎦ ⎣Y21 Y22 ⎦ ⎣VR ⎦ where ei is a vector with all entries equal to zero except one that equals 1.0 corresponding to the measured current, Ir and Ii are the real and imaginary parts of the measured current and η is the measurement error. 2.3.3.2 Pseudo-measurements In addition to the physical measurements listed above, a number of pseudo-measurements are introduced as described below. Pseudo-measurements of the voltages at other end of lines (neighboring substations): These pseudo measurements are illustrated in Figure 2-4 and are introduced as follows: ~ ~ Given measurements of VS (all three phases) and I S (all three phases) of a line i at a substation k and a 3-phase model of the line, of the generalized pi-equivalent form: ~ ⎡ I S ⎤ ⎡Y11 ⎢~ ⎥ = ⎢ ⎣ I R ⎦ ⎣Y21 ~ Y12 ⎤ ⎡VS ⎤ ⎢~ ⎥ Y22 ⎥ ⎣V R ⎦ ⎦ Then the voltage pseudo-measurement at the other end of the line (neighboring substation) is defined by: ~ ~ ~ −1 −1 V Rpseudo , m = (I − Z 22Y22 ) Z 21 I S + (I − Z 22Y22 ) Z 22Y21VS ⎡Z where ⎢ 11 ⎣ Z 21 Z 12 ⎤ ⎡Y11 = Z 22 ⎥ ⎢Y21 ⎦⎣ Y12 ⎤ and I is the identity matrix of appropriate dimension. Y22 ⎥ ⎦ −1 19 ~ VR ~ IS ~ IR Li ne i ~ VS Substation k Figure 2-4: Pseudo-measurements of Voltage at Other End of Lines Pseudo-measurements from Kirchoff’s current law: These pseudo-measurements are illustrated in Figure 2-5 and are introduced as follows: In general it holds that ∑k I ~ = 0, where k is the voltage level (kV level), at each internal node of the substation. Therefore, for the substation shown in Figure 11, the following equations define pseudo-measurements: ~~~ I1 + I 2 + I 6 = 0 : Three equations, one for each phase ~~~ I3 + I4 + I5 = 0 : Three equations, one for each phase ~~ ~~ ~ k1 I3 +I4 +k2 I1 +I2 +Im =0:Three equations, one for each phase ( )( ) ~ where I m is the transformer magnetizing current. Since its value is rather small, it can be neglected without significant loss of accuracy. The coefficients k1 and k2 are the two voltage levels that define the transformation ratio of the power transformer. 20 ~ I4 ~ I5 ~ I3 Substation ~ I6 ~ I1 ~ I2 Figure 2-5: Pseudo-measurements from Kirchoff’s Current Law Pseudo-measurements of the individual phase voltages: These pseudo-measurements are introduced if the measurement of a single phase is missing or is not available, and they are defined as follows: Assuming that phase A voltage measurement exists, while phase C voltage measurement does not. In this case the following equation is introduced, defining a pseudo-measurement equation for phase C: 0 ~ ~ Vc pseudo ,m = Va e − j 240 or: 0 ~ ~ Vc pseudo , m − V a e − j 240 = 0 Pseudo-measurements of the neutral/shield wire current: These pseudo measurements are introduced as follows (see Figure 2-6): Given a line model we can compute the ratio of shield/neutral current over the return current as: ~ Is/n α= ~ ~ ~ . − I a + Ib + Ic ( ) Then the following pseudo-measurement is introduced: 21 ~ ~~~ ~ ~~~ I s pseudo,m = −α I a + I b + I c or I s pseudo,m + α I a + I b + I c = 0 . /n /n ( ) ( ) s/n a b c Figure 2-6: Pseudo-measurements from Neutral/Shield Wire Current. Pseudo-measurements of the neutral/ground voltage: These pseudo measurements are introduced as follows: The neutral/ground voltage is computed as the product of the substation ground impedance and the earth current. The earth current is the sum of all earth currents of all transmission lines connected to the substation. The equations are: ~ ~ ~ ~ Vs pseudo,m = − Rg ∑ (1 − α i ) I a ,i + I b,i + I c ,i /n i ( ) where Rg is the substation ground resistance and the coefficient αi is the ratio of shield/neutral current over the return current. Note that under normal operating conditions the neutral/ground voltage will be negligibly small. However during a fault condition or a stuck breaker the neutral/ground voltage may be elevated to a substantial value. 2.3.4 Dynamic State Estimator (DSE) Algorithm The relationship of each measurement to the system dynamic state is expressed in the following general form (see previous description of measurements and states): z m = h( x, y, t ) + η m Due to the existence of dynamic states, the relationship between the states and the measurements is not necessarily linear. However, due to the model quadratization it is in general quadratic of the following form (for linear measurements the matrix Q is zero): ⎡x ⎤ ⎧ z m = H m ⎢ ⎥ + ⎨[x ⎣ y⎦ ⎩ ⎡x ⎤⎫ y ].Qm ⎢ ⎥ ⎬ + η m ⎣ y⎦⎭ 22 The problem is now defined as to calculate the estimate of xe and ye in order to minimize the following objective function which is defined as the weighted sum of the errors squared: J ( x, y ) = η T W η , where: η = z − h ( x, y ) , ⎡1⎤ W = diag ⎢ 2 ⎥ . ⎣σ v ⎦ and W is a diagonal matrix whose non-zero entries are equal to the inverse of the variance of the measurement errors: The state of the system will be derived from the following iterative solution based on the Gauss-Newton method: ˆ ˆ ˆ x j +1 = x j + ( H T WH ) −1 H T W ( z a − h( x j )) , ˆ where x refers to the overall vector of the states: ⎡ xe ⎤ ˆ x = ⎢ ⎥, ⎣ ye ⎦ and H is the Jacobian matrix of the measurement equations, defined in this case as: ⎡x ⎤ H = H m + 2 ⋅ Qm .⎢ ⎥ ⎣ y⎦ Note that above equations are equivalent to a static estimation problem. The reason for this is that the differential equations have been integrated and converted into a set of algebraic equations. Then the estimation algorithm is applied to a set of algebraic equations. This process significantly simplifies the implementation of the dynamic state estimator. 23 3 Stability Monitoring 3.1 Introduction Stability monitoring and control is based on the dynamic state estimator. By necessity the dynamic state estimator is performed at the substation level to avoid time latencies that will be associated with data transfer from substation to control center. It should be understood that the data eventually are transmitted to the control center but with the time delay associated with the transfer of data. The process is more meaningful at generating substations and therefore the discussion below assumes that a generating station is being monitored. The estimated dynamic state of the system provides the necessary information for monitoring the stability of the system and as a consequence it provides an elegant method for (a) monitoring the stability margin in case of a stable system, and (b) predicting out of step conditions in case of an unstable swing of the system. As it will be shown the method predicts instability before it happens and therefore can be used for proper control of the system. The real time dynamic state of the substation includes the full operating condition of the generators and can be utilized to characterize and predict the stability of the system. Specifically this model contains the generator speed, torque angle and acceleration. This information is enough to monitor the dynamics of the generator. We have investigated the use of energy functions to predict the stability of the system. For this purpose we use the output of the real time dynamic model of the system to compute the total energy of the system. Note that the total energy of the system is defined in terms of generator torque angle and speed. The energy function is a quadratic equation of the form (defined in the usual way of the sum of the potential energy plus the kinetic energy): T V ( x (t )) = P1T x N + x N (t ) P2 x N (t ) Note that this form includes the state variables c and s see appendix which are transcendental equations (sines and cosines). The generator model in terms of these variables is presented in Appendix B. The above form is quite complex and the complexity hides the physical meaning of this approach. For this reason, we present a simplified explanation of the approach using classical notation. Te simplified explanation outlines the application of Lyapunov’s direct method to the transient stability problem. We first introduce the concept of generator potential energy, then its kinetic energy and then the total energy. The total energy of a generator is used as the Lyapunov test function. This explanation is not intended to cover this complex topic in great detail but rather to provide a simplified overview of the basic concepts and equations involved. 24 3.2 Introductory Example It is well known that the total energy of a system can be used as a Lyapunov function to study the stability of the system. We introduce this concept with a mechanical system and then apply it to a single generator power system. Example: Consider a metallic ball which can slide frictionless in the inner surface of a cylinder as it is illustrated in Figure 3-1. Determine the stability of this system. r mg h Figure 3-1: Illustration of a Metallic Ball in a Cylinder The total energy of the metallic ball is V (θ (t ), t ) = mgh + 1 mv 2 (t ) 2 The total energy can be expressed as a function of the position θ (t ) (angle from the vertical axis) of the metallic ball as follows. ⎛ θ (t ) ⎞ h = 2 r sin 2 ⎜ ⎟ ⎝2⎠ dθ ( t ) v (t ) = r dt Upon substitution of above expressions: ⎛ θ (t ) ⎞ 1 2 ⎛ dθ (t ) ⎞ V (θ (t ), t ) = 2mgr sin ⎜ ⎟ + mr ⎜ ⎟ ⎝2⎠ 2 ⎝ dt ⎠ 2 2 25 On the other hand, the motion of the system under gravity forces is governed by the equation: m dv (t ) ⎛π ⎞ = −mg cos⎜ − θ (t ) ⎟ dt ⎝2 ⎠ Upon substitution of v in terms of q and some minor manipulations d 2θ ( t ) g = − sin (θ ( t ) ) 2 dt r Now consider the time derivative of V: dV (θ (t ), t ) dθ (t ) d 2θ (t ) ⎛ θ ( t ) ⎞ ⎛ θ ( t ) ⎞ dθ ( t ) cos⎜ + mr 2 ≡0 = 2mgr sin ⎜ ⎟ ⎟ dt dt dt 2 ⎝ 2 ⎠ ⎝ 2 ⎠ dt Note that the time derivative of the Lyapunov function is zero while the function itself is nonnegative and decrescent. By Lyapunov’s theorem, the system is stable. 3.3 Generator “Potential” Energy In general, Lyapunov’s theory does not provide a systematic way for determining the Lyapunov function for a general nonlinear system as that of the electric power system. On the other hand, it is clear that Lyapunov’s functions are energy like functions, i.e. they express total energy of the system. It is natural then to try to define the total energy of a synchronous machine as a Lyapunov function. Obviously the total energy is constant during steady state operation and it is time varying during transients. We introduce the concept of potential energy of a synchronous generator as follows. Assume that the equilibrium position of a generator is at δ= δs. We measure this angle from the center of oscillation (CoO) of the system. The CoO is determined by the entire network. We presented earlier a method for identifying the CoO in real time by an estimation procedure. As we will see next, knowledge of the CoO allows us to express the generator potential energy with one variable only. We assume that the generator position deviates from the equilibrium to an arbitrary position δ and the transition takes place very slow so that the speed of the generator is practically constant equal to the synchronous frequency. At the new position there will be a net torque, T, acting on the generator. In general the torque T is a function of the position of the rotor (at δ= δs the net torque equals zero). The potential energy equals the work done to move the generator from position δs to position δ. This work is given with the integral: E p = ∫ (T − T0 )dδ = ∫ (Pacc )dδ , transitions with ω (t ) = ωs δs δs δ δ 26 The above equation assumes that T and P are expressed in pu. Note that knowledge of the accelerating power as a function of delta, the above integral can be computed. 3.4 Generator “Kinetic” Energy The generator kinetic energy is defined as the total energy stored in the generator rotating part, i.e. the rotor. It is important to note that the classical transient stability problem is formulated in terms of the speed deviation from the synchronous speed. Therefore the variable omega represents the deviation from the synchronous speed. Looking at the problem from the strict mathematical problem, the generator kinetic energy (for the purpose of defining the Lyapunov test function) is: Ek = 1 2H Mω 2 (t ), M = 2 ωs 3.5 Generator “Total” Energy Consider a single generating unit system that experienced a disturbance (for example a fault). At the end of the disturbance, the generator is at a state described with a certain position δ 0 = δ (t = t c ) and certain speed ω 0 = ω (t = t c ) . Further assume that if the system goes to equilibrium at position δ s . This is expressed with the following model: dδ (t ) = ω (t ) dt M dω ( t ) = Pm − f post (δ (t )) dt The initial conditions at time t = t c are: δ (t = t c ) = δ 0 , and ω (t = t c ) = ω 0 . The equilibrium point or the steady state solution is δ s , i.e. Pm − f post (δ s ) = 0 . The Lyapunov test function for this system is defined as the sum of the kinetic and potential energy of this generator, i.e. V (ω (t ), δ (t )) = 1 1 Mω 2 (t ) + ∫ ( f post (δ (t )) − Pm )dδ = Mω 2 (t ) − Pmδ (t ) + ∫ ( f post (δ (t )) )dδ + k 2 2 δs δs δ δ where k is a constant. The total time derivative of the defined function can be computed as follows: 27 dV (ω (t ), δ (t ) ) dω ( t ) dδ ( t ) dδ ( t ) ⎛ dω ( t ) ⎞ = Mω ( t ) − Pm + f post (δ (t ) ) = ω ( t )⎜ M − Pm + f post (δ (t ) )⎟ = 0 dt dt dt dt dt ⎝ ⎠ By virtue of Lyapunov’s theorem, whenever the function V (ω (t ), δ (t )) is positive the system is stable. Therefore the immediate conclusion is that if: V (ω 0 , δ 0 ) > 0 , then the system is stable. In general we can generate equal-value graphs of the function V (ω (t ), δ (t )) in the (δ (t ), ω(t )) plane. The equal-value graph of V (ω(t ), δ (t )) = 0 separates the (δ (t ), ω(t )) plane into two regions (as a matter of fact into multiple regions but we generally are not interested in all regions) one where the margin is positive and another at which the margin is negative. This curve is known as the separatrix. It separates the stable region from the unstable region. Figure 3-2 illustrates a separatrix for a specific system. Figure 3-2: Illustration of the Separatrix for a Specific System Another question that can be answered with the Lyapunov test function is the question of the critical clearing time. For this purpose one needs to consider the system model during the disturbance. During that period, the model should be solved (close form solution or by time domain simulation) to provide the system state (δ (t ), ω(t ) ) as a function of time. By substituting the solution into the Lyapunov test function, one can compute its value as function of time. The time at which the Lyapunov function assumes a value equal to the value of the potential well it will be the critical clearing time. This is illustrated in Figure 3-3. 28 Figure 3-3: Illustration of Monitoring and Identifying the Critical Clearing Time by the Total Energy 29 4 Case Study 4.1 Demonstration Example The proposed real time stability monitoring and instability prediction method has been applied to several systems. In this section we present an example that illustrates the method and the various features of the method. The example system used is kept very small for the purpose of amplifying the features of the method. The test system consists of two generating substations, two step-up transformers and a two section overhead transmission line as is illustrated in Figure 4-1. The parameters of the system are presented in Table 4-1. 1 GEN1 2 GEN1H GEN1H GEN2H 1 2 GEN2 P/Q P/Q P/Q P/Q P/Q P/Q Figure 4-1: Single Line Diagram of Example Test System Table 4-1: Test System Parameters Gen1 100MVA Gen2 300MVA XFMR1 100MVA XFMR2 300MVA Transmission Line Load 1 Load 2 z= 0.001+j0.18 pu z= 0.001+j0.18 pu z=0.001+j0.07 pu z=0.001+j0.08 pu z=0.024+0.2346 pu S=1.33+j0.1 pu S=0.56+j0.1 pu Common Sbase=100 MVA H=2.5 sec H= 3.0 sec 15 kV/115kV 115 kV/18kV 115 kV 115 kV 115 kV It is assume that there is a PMU at the substation to the left of the figure that measures the generator terminal voltages and currents as well as the high side transformer votages and currents and the associated frequency and rate of change of frequency. Through the SuperCalibrator, the data collected at the substation are filtered and in addition an estimate of the voltages and their frequency is estimated at the substation on the right of the figure. Subsequently, the Lyapunov theory is used to monitor and evaluate the stability of the system during normal operating conditions and/or after the occurrence of a 30 fault, using information from the estimated dynamic model of the system of the generator of interest and the next substation. The stability monitoring procedure can be separated into two distinct cases: (a) stability monitoring in the absence of faults and (b) stability monitoring in the presence of faults. These cases are discussed separately. Stability Monitoring in the Absence of Faults: In this case the system operates in near normal condition with small oscillations around the equilibrium point. Of interest here is to determine the stability margin of the system. This margin is defined in terms of the potential energy margin. In particular, the potential energy function is computed as a function of the difference between the angle of the generator of interest and the reference angle of the equivalent system beyond the next substation. The potential energy function is computed from the best estimate of the system dynamic state via a three step procedure: step 1: the electric power model of the system is computed as a function of generator rotor position, step 2: the equilibrium point of the system is computed, step 3: the potential energy function is computed by integrating the electric power model. For the example system we have the following electric power model: 2 H 1 d δ 12 = Pm1 − Pe1 (δ 1 − δ 2 ) ωs dt 2 2 H 2 dδ 22 = Pm 2 − Pe 2 (δ1 − δ 2 ) ωs dt 2 The generator real power equations are evaluated to be: Pe1(δ1 − δ2 ) = P1(δ ) = 0.5222 − 1.75748 ⋅ cos(δ + 100.008) e P 2 (δ1 − δ2 ) = P 2 (δ ) = 0.7248 −1.7592 ⋅ cos(−δ + 100.008) e e Above equations are manipulated to provide: M d 2δ dt 2 =M ωs 2 H1 (P m1 − Pe1 (δ ) ) − M ωs 2H 2 (P m2 − Pe 2 (δ ) ) where M is the two generators equivalent mass: 2 H1 ⋅ H 2 M= ( H1 + H 2 ) ⋅ ωs Upon substitution the following swing equation is obtained. = 0.29563 + 1.375363 ⋅ cos(δ + 100.008) dt 2 −0.0598 − 0.38242 ⋅ cos( −δ + 100.008) M The equilibrium point is computed to be at: d 2δ δs = δs1 − δs 2 = 1.61420 at the synchronous speed. 31 Note that above model is provided in real time from the dynamic state estimator. Next the potential energy is evaluated: E potential = − δ s1 − ∫δ δ [M s2 ωs 2 H1 (P m1 − Pe1 (δ ) ) − M ωs 2H 2 (P m2 − Pe 2 (δ ) )] which for the specific system is given by the following formula: V (0, δ ) = −0.23583 ⋅ δ + 1.73214 −1.375363 ⋅ sin(δ + 100.008) − 0.38242 ⋅ sin( −δ + 100.008) The maximum value of this function is 2.9. Since the potential energy at the equilibrium point is by definition equal to zero, then the margin is 2.9. Similarly the margins can be computed at any other operating condition. The advantage of this approach is that it provides a direct measure of the distance from the instability point and therefore it is a better index for monitoring system stability. It is important to correlate the stability margin to the values of the generator rotor position. This will be the subject of additional investigations. Stability Monitoring in the Presence of Faults: In this case the system operates in near normal condition when a fault occurs in the system. During fault conditions, the system relaying will identify the fault and it will isolate the fault by disconnecting the faulty device. It is important to monitor the stability of the system during this event. It is important to recognize that the sequence of events in this case is very fast and we cannot count on operator action to correct the problem. What is shown though is the fact that by applying the proposed method we can predict the instability (in case of an unstable swing) before the system is at a point beyond return. This prediction can be utilized in automated response systems, for example to provide better relaying action. Of interest here is to determine whether the transient will be a stable transient or an unstable transient. This is achieved by monitoring the total energy of the system versus the potential energy of the system. The potential energy function is computed from the best estimate of the system dynamic state via a four step procedure: step 0: determine the anticipated relay action and the post fault system topology, step 1: the electric power model of the post-fault system is computed as a function of generator rotor position, step 2: the equilibrium point of the post-fault system is computed, step 3: the potential energy function is computed by integrating the electric power model. For the example system we have the following electric power model (for this system we have deliberately selected a fault that the post-fault topology is the same as the prefault): 2 H 1 d δ 12 = Pm1 − Pe1 (δ 1 − δ 2 ) ωs dt 2 32 2 H 2 dδ 22 = Pm 2 − Pe 2 (δ1 − δ 2 ) ωs dt 2 The generator real power equations are evaluated to be: Pe1(δ1 − δ2 ) = P1(δ ) = 0.5222 − 1.75748 ⋅ cos(δ + 100.008) e P 2 (δ1 − δ2 ) = P 2 (δ ) = 0.7248 −1.7592 ⋅ cos(−δ + 100.008) e e Above equations are manipulated to provide: M d 2δ dt 2 =M ωs 2 H1 (P m1 − Pe1 (δ ) ) − M ωs 2H 2 (P m2 − Pe 2 (δ ) ) where M is the two generators equivalent mass: 2 H1 ⋅ H 2 M= ( H1 + H 2 ) ⋅ ωs Upon substitution the following swing equation is obtained. = 0.29563 + 1.375363 ⋅ cos(δ + 100.008) dt 2 −0.0598 − 0.38242 ⋅ cos( −δ + 100.008) M d 2δ The equilibrium point is computed to be at: δs = δs1 − δs 2 = 1.61420 at the synchronous speed. Note that above model is provided in real time from the dynamic state estimator. From this model, the potential energy of the post-fault system is: E potential = − δ s1 − ∫δ δ [M s2 ωs 2 H1 (P m1 − Pe1 (δ ) ) − M ωs 2H 2 (P m2 − Pe 2 (δ ) )] which for the specific system is given by the following formula: V (0, δ ) = −0.23583 ⋅ δ + 1.73214 −1.375363 ⋅ sin(δ + 100.008) − 0.38242 ⋅ sin( −δ + 100.008) The kinetic energy of the system is: Ek = 1 M (ω1 − ω2 ) 2 2 where δ=δ1-δ2 is the difference of the generators’ torque angles and ω1, ω2 are the generators’ speeds. The total energy of the system is: 33 E total = V (0, δ ) + E k Two example cases are presented for this system, one corresponding to a disturbance that results in a stable system and another that leads to an unstable system. Specifically, a three phase fault that will be cleared at 0.5 seconds after the fault initiation results in a stable system. We shall refer to this case as “stable case”. A three phase fault that will be cleared at 0.6 seconds after fault initiation will result in an unstable system. We refer to this case as the “unstable case”. Figure 4-2 illustrates the total energy superposed on the graph of the potential energy function of the “stable case” at the time that the fault was cleared. Note that the total energy is below the highest value of the potential energy, indicating a stable system. Figure 4-3 illustrates the total energy superposed on the graph of the potential energy function of the “unstable case” at the time that the fault was cleared. Note that the total energy is above the highest value of the potential energy, indicating an unstable system. Table 4-2 provides the rotor position and rotor speed at the time of fault clearing for the two cases. Figure 4-2: Total and potential system energy-stable case. 34 Figure 4-3: Total and potential system energy-unstable case Table 4-2: Generators’ Torque Angles & Speeds at Fault Clearing Time. Stable Case δ=δ1-δ2=108.17ο ω=ω1-ω2=7.3953 rad/sec Unstable Case δ=δ1-δ2=155.0280 ω=ω1-ω2=8.8404 rad/sec An important feature of the presented method is the ability to determine the time at which the system becomes unstable since it is a real time monitoring method. Figure 4-4 illustrates the trajectory of the total energy of the system as the disturbance is evolving. The instant that the total energy crosses the maximum value of the potential energy the system has become unstable. The total energy crosses the maximum value of the potential energy at 0.545 seconds after fault initiation. Thus the critical clearing time for this system is 0.545 seconds. 35 Figure 4-4: Test System Total Energy Trajectory. Note that the time involved in the development of the instability is very short and beyond the human response (or control center operator). However this information can be utilized in a relaying scheme to initiate appropriate action. A very important application would be to design a new out of step relay function. As compared to present out of step relaying this approach will identify the out of step condition much earlier than conventional techniques. This application will be further pursued. This example demonstrated the methodology for filtering available data in substations (for example phasor data, relay data and SCADA data) for the purpose of extracting a real time dynamical model of the system. The real time dynamical model is used for monitoring system stability and it is capable to predict possible instability that may arise. The enabling technology is based on the SuperCalibrator and the associated dynamic state estimator that provides the real time dynamic model of the system. The dynamic state estimator uses a three-phase, breaker oriented, instrumentation inclusive and generator dynamics inclusive substation model. The proposed approach provides a precise dynamic state estimator for power systems at the substation level. The estimator is ideally suited for monitoring the stability of the system as well as predicting possible instability before it occurs, utilizing a Lyapunovbased energy function. 36 5 References 5.1 Papers Resulted from this Project Sakis Meliopoulos, Vagelis Farantatos, George Cokkinides, Salman Mohagheghi and George Stefopoulos, ”A New Out of Step Protection Scheme via GPS Synchronized Data”, Paper submitted to the PSCC 2008, Glasgow, Great Britain, July 2008. 5.2 Background Papers [1] A. P. Sakis Meliopoulos and G. J. Cokkinides, ”Visualization and Animation of Instrumentation Channel Effects on DFR Data Accuracy”, Proceedings of the 2002 Georgia Tech Fault and Disturbance Analysis Conference, Atlanta, Georgia, April 29-30, 2002 [2] A. P. Sakis Meliopoulos and George J. Cokkinides, ‘A Virtual Environment for Protective Relaying Evaluation and Testing’, Proceedings of the 34st Annual Hawaii International Conference on System Sciences, p. 44 (pp. 1-6), Wailea, Maui, Hawaii, January 3-6, 2001 [3] A. P. Meliopoulos, F. Zhang, S. Zelingher, G. Stillmam, G. J. Cokkinides, L. Coffeen, R. Burnett, J. McBride, 'Transmission Level Instrument Transformers and Transient Event Recorders Characterization for Harmonic Measurements,' IEEE Transactions on Power Delivery, Vol 8, No. 3, pp 1507-1517, July 1993 [4] B. Fardanesh, S. Zelingher, A. P. Sakis Meliopoulos, G. Cokkinides and Jim Ingleson, ‘Multifunctional Synchronized Measurement Network’, IEEE Computer Applications in Power, Volume 11, Number 1, pp 26-30, January 1998. [5] A. P. Sakis Meliopoulos and G. J. Cokkinides, ”Phasor Data Accuracy Enhancement in a Multi-Vendor Environment”, Proceedings of the 2005 Georgia Tech Fault and Disturbance Analysis Conference, Atlanta, Georgia, April 25-26, 2005 [6] F. Darema, “Dynamic Data Driven Application Systems,” Presentation at Purdue University, May 4, 2004, http://www.cise.nsf.gov/eia/dddas [7] L. Tsoukalas, R. Gao, T. Fieno, X. Wang, “Anticipatory Regulation of Complex Power Systems,” Proc. of European Workshop on Intelligent Forecasting, Diagnosis and Control- IFDICON 2001, Santorini, Greece, June 2001 [8] IEEE Std 1344-1995, IEEE Standard for Synchrophasors for Power Systems, 1995. [9] C37.118, IEEE Standard for Synchrophasors for Power Systems, 2004 Draft. [10] IEEE Std C37.2 – 1991, IEEE Standard Electrical Power System Device Function Numbers and Contact Designations. [11] IEEE Std 1379-1997, IEEE Trial-Use Recommended Practice for Data Communications Between Intelligent Electronic Devices and Remote Terminal Units in a Substation. [12] IEEE C37.111-1999, IEEE Standard Common Format for Transient Data Exchange (COMTRADE) for Power Systems, June 1999. 37 [13] IRIG Standard 200-98, IRIG Serial Time Code Formats, May 1998 [14] IEC 60870, Telecontrol Equipment and Systems, Transmission Protocols. [15] A. P. Sakis Meliopoulos, F. Zhang, and S. Zelingher, 'Power System Harmonic State Estimation,' IEEE Transactions on Power Systems, Vol 9, No. 3, pp 1701-1709, July 1994. [16] A. Arifian, M. Ibrahim, S. Meliopoulos, and S. Zelingher, ‘Optic Technology Monitors HV Bus’, Transmission and Distribution, Vol. 49, No. 5, pp. 62-68, May 1997. [17] T. K. Hamrita, B. S. Heck and A. P. Sakis Meliopoulos, ‘On-Line Correction of Errors Introduced By Instrument Transformers In Transmission-Level Power Waveform Steady-State Measurements’, IEEE Transactions on Power Delivery, Vol. 15, No. 4, pp 1116-1120, October 2000. [18] Evaluating Dynamic Performance of Phasor Measurement Units: Experience in the Western Power System, J. F. Hauer, Ken Martin, and Harry Lee. Interim Report of the WECC Disturbance Monitoring Work Group, partial draft of April 28, 2004. [19] Integrated Monitor Facilities for the Western Power System: The WECC WAMS in 2003, J. F. Hauer, W. A. Mittelstadt, K. E. Martin, and J. W. Burns. Interim report of the WECC Disturbance Monitoring Work Group, June 25, 2003. (Available at http://www.wecc.biz/committees/JGC/DMWG/documents/). [20] "Performance of 'WAMS East' in Providing Dynamic Information for the North East Blackout of August 14, 2003", J. F. Hauer, Navin Bhatt, Kirit Shah, and Sharma Kolluri. IEEE/PES Panel on Major Grid Blackouts of 2003 in North America and Europe, IEEE PES General Meeting, Denver, CO, June 6-12, 2004. [21] "Dynamic Signatures Recorded on the "WAMS East" Monitor Backbone: August 14, 2003 vs. Other Days", J. F. Hauer, Navin Bhatt, Rajesh Pudhota, Sujit Mandal, and Michael Ingram. Presented by J. F. Hauer to the IEEE/PES IEEE Work Group on Power System Dynamics Measurements, IEEE PES General Meeting, Denver, CO, June 7, 2004. [22] Juancarlo Depablos, Virgilio Centeno, Arun G. Phadke, and Michael Ingram, “Comparative Testing of Synchronized Phasor Measurement Units”, Paper presented at the IEEE/PES General Meeting, Denver, CO, June 2004. [23] S. Zelingher, G.I. Stillmann, A. P. Sakis Meliopoulos, “Transmission System Harmonic Measurement System: A Feasibility Study," Proceedings of the Fourth International Conference on Harmonics in Power Systems (ICHPS IV), pp. 436-444, Budapest, Hungary. October 1990. [24] A. P. Sakis Meliopoulos, F. Zhang, and S. Zelingher, "Hardware and Software Requirements for a Transmission System Harmonic Measurement System," Proceedings of the Fifth International Conference on Harmonics in Power Systems (ICHPS V), pp. 330-338, Atlanta, GA. September 1992. [25] Real-time applications task team, EIPP goals and objectives [26] EIPP evaluating dynamic performance of phasor measurement units, March 2004 [27] EIPP Standards and Performance Task Team preliminary report, August 2004 [28] NASPI PRTT Testing and Interoperability Guide, 2007. 38 [29] Arun Phadke, EIPP Presentation, Washington, DC, October 2006. [30] R. Jay Murphy, “Fault Recorder Integration”, Proceedings of the Georgia Tech Fault and Disturbance Analysis Conference, Atlanta, GA, May 3-4, 1999. [31] A. P. Meliopoulos, G. J. Cokkinides, Floyd Galvan, Bruce Fardanesh and Paul Myrda, ”Advances in the SuperCalibrator Concept – Practical Implementations”, Proceedings of the of the 40st Annual Hawaii International Conference on System Sciences, Kona, Hawaii, January 3-6, 2007. [32] A. P. Meliopoulos, G. J. Cokkinides, Floyd Galvan and Bruce Fardanesh, “Distributed State Estimnator – Advances and Demonstration”, Proceedings of the of the 41st Annual Hawaii International Conference on System Sciences, Kona, Hawaii, January 7-10, 2008. [33] Mike Ingram, Sandra Bell, Sherica Matthews, A. P. Sakis Meliopoulos and G. J. Cokkinides, ”Use of Phasor Measurements, SCADA and IED Data to Improve State Estimation Procedures”, Proceedings of the 2004 Georgia Tech Fault and Disturbance Analysis Conference, Atlanta, Georgia, April 26-27, 2004. 39 Appendix A. Quadratic Integration Method A.1 Introduction This section of the proposal describes a new proposed approach for power system time domain simulation. The new methodology, referred to as quadratic integration method, is based on the following two innovations: (a) the nonlinear system-model equations (nonlinear differential or differential-algebraic equations) are reformulated to a fully equivalent system of linear differential and quadratic algebraic equations, by introducing additional state variables and algebraic equations, and (b) the system model equations are integrated using a numerical integration scheme that assumes that the system states vary quadratically within an integration time step, as functions of time (quadratic integration). As a comparison, in the trapezoidal rule it is assumed that the system functions/states vary linearly throughout a time step. This basic concept in the derivation of the quadratic integration method is illustrated in Figure A. 1. Note that within an integration time step of length h , defined by the interval [t − h, t ] , the two end points, x(t − h) , x(t ) , and the midpoint x(t m ) = x(t − h / 2) fully define the quadratic function in the interval [t − h, t ] . This quadratic function is integrated in the time interval [t − h, t ] resulting in a set of algebraic equations for this integration step. The solution of the equations is obtained via Newton’s method, in the general, nonlinear system case, or via a direct solution in the linear case. Note that by virtue of the first step of quadratization the resulting algebraic equations are either linear or quadratic. The proposed method demonstrates improved convergence characteristics of the iterative solution algorithm. It is an implicit numerical integration method (it can be easily observed that it makes use of information at the unknown point x(t ) ) and therefore demonstrates the desired advanced numerical stability properties compared to explicit methods. It is also fourth order accurate and therefore much more precise compared to all the traditionally used methods in power system applications. Furthermore, the proposed method does not suffer from the numerical oscillation problem, contrary to the trapezoidal rule. 40 x x(t) xm Actual Curve Trapezoidal x(t-h) Quadratic t t-h 0 tm h/2 t h τ Figure A. 1 Illustration of the Quadratic Integration Method - Basic Assumption. A.2 Description of Quadratic Integration Method This section presents the key features of the quadratic integration method. The method is based on two innovations: First, the nonlinear system-model equations (nonlinear differential or differential-algebraic equations) are reformulated to a fully equivalent system of linear differential and quadratic algebraic equations, by introducing additional state variables and additional algebraic equations. This step aims in reducing the nonlinearity of the system to at most quadratic in an attempt to improve the efficiency of the solution algorithm, as well as facilitate the method implementation, especially for large complex models. It is independent of the integration method and thus can be applied in combination with any numerical integration rule. Second, the system model equations are integrated using the implicit numerical scheme that was conceptually described in the introductory section of this chapter. The quadratic integration method belongs to the category of implicit, one-step, RungeKutta methods. More specifically it is an implicit Runge-Kutta method based on collocation and it can be alternative derived based on the collocation theory. The basic idea is to choose a function from a simple space, like the polynomial space, and a set of collocation points, and require that the function satisfy the given problem equations at the collocation points [21]-[23]. The method has three collocation points, at x(t − h) , x m = x(t m ) , and x(t ) . It uses the Lobatto quadrature rules and is a member of the Lobattomethods family. Any Lobatto method with s collocation points has an order of accuracy 41 of 2s − 2 , and therefore the method is order four accurate [21]-[23]. Assuming the general nonlinear, non-autonomous dynamical system: & x = f (t , x) (B.1) the algebraic equations at each integration step of length h , resulting from the quadratic integration method, are: h h 5h x m − f (t m , x m ) + f (t , x(t )) = x(t − h) + f (t − h, x (t − h)) 3 24 24 h h 2h x(t ) − f (t m , x m ) − f (t , x(t )) = x(t − h) + f (t − h, x(t − h)) 3 6 6 Solution of the system, via Newton’s method in practice, yields the value of the state vector x(t) . Note that the value at the midpoint, x m , is simply an intermediate result and it is discarded at the end of the calculations at each step. For the special case of a linear system, & x = Ax + Bu , the algebraic equations at each time step become: ⎡⎛ 5h ⎞ ⎤ h⎤ ⎡h ⎡h I+ A⎟ B ⎢ 24 A I − 3 A⎥ ⎡ x(t ) ⎤ ⎢⎜ ⎢− 24 ⎠⎥ ⎝ ⎥ ⋅ x(t − h) + ⎢ 24 ⎢ ⎥⋅⎢ ⎥=⎢⎛ 2h h h h ⎢I − A − ⎢B A ⎥ ⎣ xm ⎦ ⎢ ⎜ I + A⎞ ⎥ ⎟⎥ ⎢⎝ 6 3⎦ 6 ⎣ ⎣ 6 ⎠⎦ ⎣ 5h B 24 h B 6 h ⎤ ⎡ u (t ) ⎤ B 3 ⎥ ⋅ ⎢u (t − h)⎥ , ⎥ 2h ⎥ ⎢ B⎥ ⎢ ⎥ 3 ⎦ ⎣ um ⎦ where I is the identity matrix of proper dimension and h the length of the integration step. The proposed integration approach has the following advantages: (a) improved accuracy and numerical stability, and (b) free of fictitious numerical oscillations. Details about the numerical properties of the method are discussed next. A.3 Numerical properties The basic numerical stability properties of a numerical integration method are studied using the first order test equation: & x = ax . Applying the quadratic integration method yields at each time step: ⎡⎛ 12 + 6ah + a 2 h 2 ⎞⎤ ⎟⎥ ⎢⎜ ⎡ x(t )⎤ ⎢⎜ 12 − 6ah + a 2 h 2 ⎟⎥ ⎝ ⎠ x(t − h) ⎢ x ⎥ = ⎢⎛ 22 ⎣ m ⎦ ⎜ 12 − 0.5a h ⎞ ⎥ ⎟⎥ ⎢⎜ 12 − 6ah + a 2 h 2 ⎟ ⎦ ⎠ ⎣⎝ 42 and therefore: 12 + 6ah + a 2 h 2 ⋅ x (t − h ) , 12 − 6ah + a 2 h 2 where h is the integration step. Setting z = ah yields the characteristic polynomial for the method: z 2 + 6 z + 12 R( z ) = 2 z − 6 z + 12 Note that the eigenvalue a of the system can be complex, so z is in general a complex number. x (t ) = The region of absolute stability is given by the set of values z such that R( z ) ≤ 1 . A method is called A-stable if the region of absolute stability in the complex z-plane contains the entire left half plane. This means that independently of the step size h > 0 , a stable eigenvalue a of the original continuous time system, with Re(a) < 0 , will be still represented as a stable mode in the discrete time system, and thus the discrete system mimics accurately the behavior of the original system, in terms of stability. Note that for Re( z ) < 0 it follows that R ( z ) ≤ 1 . Therefore, the proposed method is A-stable. Furthermore, the absolute stability region is exactly the left-hand half complex plane. This property is called strict A-stability. If the dynamical system under study includes an unstable mode, then, irrespectively of the integration step-size, this mode will remain unstable in the descretized system. This is not the case for other methods, for example, the backward Euler, or the BDF linear multi-step methods, where the numerical stability domain extends in the right-hand plane, where Re( z) > 0 . In this case, if the real dynamical system includes an unstable mode, this mode could appear as stable for some step size, in the discrete system. Comparing the quadratic and the trapezoidal integration methods the following hold: 1. Both the trapezoidal method and the quadratic integration method are strictly A2+ z stable. The characteristic polynomial for the trapezoidal method is R ( z ) = , 2− z and it holds that R( z ) ≤ 1 in the whole left-hand complex plane, i.e., Re( z ) < 0 . 2. The trapezoidal method is second order accurate. The quadratic integration is a fourth order method. Therefore, in terms of accuracy, quadratic integration is much preferable. This however, comes with the drawback of additional computational expense. 3. It has been observed in applications that the trapezoidal method can provide an oscillatory solution even for systems that have exponential solutions as the simple 2+ z test equation above. This is apparent if one considers the term R ( z ) = for a 2− z physically stable system. Note that it is possible to select the integration time step 43 ( z = ah ), so that this term is negative (for example any real value for z , with z < −2 ). This can occur when larger integration steps are selected. In this case the solution will be oscillatory, oscillating around the true solution of the problem. In the z 2 + 6 z + 12 case of the quadratic integration, the corresponding term R ( z ) = 2 can z − 6 z + 12 never be negative as long as Re(z) is negative, i.e. as long as the physical system is stable. Therefore, it appears that this method is free of such fictitious oscillations. This can be a very nice characteristic in many applications. The behavior of the quadratic integration method in terms of this issue is still under investigation and will be part of the future work associated with this research. As described, the proposed quadratic integration method is a one-step implicit RungeKutta method based on collocation. The trapezoidal integration method can be also viewed as a member of this category; however, trapezoidal rule uses two collocation points, while the proposed method uses three. This provides a great advantage in terms of accuracy. As every numerical integration method, the quadratic integration directly converts the system of differential equations to a set of algebraic equations, at each integration step. The formulation of these equations is straightforward and the procedure can also be automated. This can facilitate the process in more complicated models. However, the number of algebraic equations of the quadratic integration scheme is double compared to that of the trapezoidal rule, due to the additional collocation point. The end-result is increased computational effort compared to the trapezoidal method per iteration (approximately double when sparsity techniques are used). Nonetheless, the improved method accuracy (order four, compared to order two of the trapezoidal method) allows the use of larger time-steps, so that the total computational effort becomes less than that of trapezoidal integration, while the accuracy remains significantly higher. The trade-off between accuracy and computational speed applies also to higher order implicit Runge-Kutta methods. As the number of collocation points, and thus the order, increase, the computational effort also increases. It appears that the quadratic integration method achieves a good balance between accuracy and computational speed. The proposed method also appears to possess better numerical properties and be more accurate when compared to linear, multi-step methods commonly used in power system transient analysis. The use of such methods is usually restricted to order two accurate methods. Detailed comparison of the quadratic integration and linear, multi-step methods used will be reported in future work of this research. 44 Appendix B. Generator Model B.1 Introduction This appendix presents the generator model, as an example of the dynamic quadratized equations. A two axes generator model is used as it is illustrated in figure B1. The quadratized full time domain model is presented first. Consequently, the time domain model is simplified and converted to a quadratic quasi-steady state model, suitable for stability studies, by neglecting the fast electrical transients and assuming that the electric network operates under sinusoidal steady state conditions of fundamental frequency. Therefore, phasor representation is assumed, at the fundamental frequency. These phasors are time-varying phasors, since their magnitude and phase change through time. Finally a frequency domain quadratic steady steady model is also presented. It should be noted that the presented generator model is a physically based model that assumes all quantities are phase quantities. No reference frame transformation is used at any stage. Phase A Magnetic Axis (Reference) ia(t) ra q-ax is θ d-ax is D of irec Ro tio ta n tio n ω Rotor iD(t ) Stator vf(t) if(t) iQ(t ) in(t) rc rb ib(t) ic(t) Figure B-1. Two Axes Generator Model 45 B.2 Full Time Domain Model The quadratized equations of the model are (separated into extrernal, internal differential, internal linear and internal quadratic): External equations: i a = i a (t ) i b = i b (t ) i c = i c (t ) in = −ia (t ) − ib (t ) − ic (t ) i f = i f (t ) i fn = −i f (t ) Tm = Tm (t ) State: v(t ) = va (t ) vb (t ) vc (t ) vn (t ) v f (t ) v fn (t ) Tm (t ) Internal equations: Differential Equations: dλ a (t ) = ea (t ) dt dλb (t ) = eb (t ) dt dλc (t ) = ec (t ) dt dλ f (t ) = e f (t ) dt d λ D (t ) = e D (t ) dt dλQ (t ) = e Q (t ) dt d θ m (t ) = ω m (t ) dt d ω m (t ) 1 = Tacc (t ) dt J dc (t ) = y1 (t ) dt ds (t ) = y 2 (t ) dt T State: x(t ) = λa (t ) λb (t ) λc (t ) λ f (t ) λD (t ) λQ (t ) θ m (t ) ωm (t ) c(t ) s (t ) [ T [ 46 Linear equations: 0 = Tacc (t ) − Te (t ) − Tm (t ) − T fw (t ) p θ m (t ) 2 p 0 = ω (t ) − ω m (t ) 2 0 = θ (t ) − 0 = δ (t ) − θ (t ) + ω s t + π 2 0 = ea (t ) + ra ia (t ) − v a (t ) + v n (t ) 0 = eb (t ) + rb ib (t ) − vb (t ) + v n (t ) 0 = ec (t ) + rc ic (t ) − vc (t ) + v n (t ) 0 = e f (t ) + r f i f (t ) − v f (t ) + v fn (t ) 0 = eD (t ) + rD i D (t ) 0 = eQ (t ) + rQ iQ (t ) State: y (t ) = Tacc (t ) θ (t ) ω (t ) δ (t ) ea (t ) eb (t ) ec (t ) e f (t ) eD (t ) eQ (t ) [ T Nonlinear Equations: 1 (ia (t )λb (t ) − ia (t )λc (t ) + ib (t )λc (t ) − ib (t )λ a (t ) + ic (t )λ a (t ) − ic (t )λb (t ) ) 0 = Te (t ) + 3 2 0 = Twf (t ) − a − bω m (t ) − cω m (t ) 0 = Pe (t ) − ea (t )ia (t ) − eb (t )ib (t ) − ec (t )ic (t ) − e f (t )i f (t ) − eD (t )i D (t ) − eQ (t )iQ (t ) 0 = y1 (t ) + s(t ) ⋅ ω(t ) 0 = y 2 (t ) − c(t ) ⋅ ω (t ) 0 = λa (t ) − Laa (t )ia (t ) − Lab (t )ib (t ) − Lac (t )ic (t ) − Laf (t )i f (t ) − LaD (t )i D (t ) − LaQ (t )iQ (t ) 0 = λb (t ) − Lba (t )ia (t ) − Lbb (t )ib (t ) − Lbc (t )ic (t ) − Lbf (t )i f (t ) − LbD (t )i D (t ) − LbQ (t )iQ (t ) 0 = λc (t ) − Lca (t )ia (t ) − Lcb (t )ib (t ) − Lcc (t )ic (t ) − Lcf (t )i f (t ) − LcD (t )i D (t ) − LcQ (t )iQ (t ) 0 = λ f (t ) − L fa (t )ia (t ) − L fb (t )ib (t ) − L fc (t )ic (t ) − L f i f (t ) − M RiD (t ) − L fQiQ (t ) 0 = λD (t ) − LDa (t )ia (t ) − LDb (t )ib (t ) − LDc (t )ic (t ) − M Ri f (t ) − LD iD (t ) − LDQiQ (t ) 0 = λQ (t ) − LQa (t )ia (t ) − LQb (t )ib (t ) − LQc (t )ic (t ) − LQf i f (t ) − LQDiD (t ) − LQ iQ (t ) 0 = Laa (t ) − Ls − Lm c(t ) 2 + Lm s(t ) 2 2π 2π 2π 0 = Lbb (t ) − Ls − Lm cos( )c(t ) 2 + Lm cos( ) s (t ) 2 + 2 Lm sin( )c(t ) s (t ) 3 3 3 2π 2π 2π 0 = Lcc (t ) − Ls − Lm cos( )c(t ) 2 + Lm cos( ) s (t ) 2 − 2 Lm sin( )c(t ) s (t ) 3 3 3 π π π 0 = Lab (t ) + M s + Lm cos( )c(t ) 2 − Lm cos( ) s (t ) 2 − 2 Lm sin( )c(t ) s (t ) 3 3 3 47 7π 7π 7π )c(t ) 2 − Lm cos( ) s (t ) 2 + 2 Lm sin( )c(t ) s (t ) 3 3 3 2 2 0 = Lbc (t ) + M s + Lm cos(π )c(t ) − Lm cos(π ) s(t ) + 2 Lm sin(π )c(t ) s(t ) 0 = Lac (t ) + M s + Lm cos( 0 = Laf (t ) − M F c(t ) 2π 2π )c(t ) − M F sin( ) s (t ) 3 3 4π 4π 0 = Lcf (t ) − M F cos( )c(t ) − M F sin( ) s (t ) 3 3 0 = LaD (t ) − M D c(t ) 2π 2π 0 = LbD (t ) − M D cos( )c(t ) − M D sin( ) s (t ) 3 3 4π 4π 0 = LcD (t ) − M D cos( )c(t ) − M D sin( ) s (t ) 3 3 0 = LaQ (t ) − M Q s(t ) 0 = Lbf (t ) − M F cos( 2π 2π ) s (t ) + M Q sin( )c(t ) 3 3 4π 4π 0 = LcQ (t ) − M Q cos( ) s (t ) + M Q sin( )c(t ) 3 3 State: 0 = LbQ (t ) − M Q cos( ⎡Te (t ) T fw (t ) Pe (t ) y1 (t ) y 2 (t ) i a (t ) ib (t ) ic (t ) i f (t ) i D (t ) iQ (t ) Laa (t ) Lbb (t ) Lcc (t )⎤ z (t ) = ⎢ ⎥ ⎢ Lab (t ) Lac (t ) Lbc (t ) Laf (t ) Lbf (t ) Lcf (t ) LaD (t ) LbD (t ) LcD (t ) LaQ (t ) LbQ (t ) LcQ (t ) ⎥ ⎣ ⎦ T Implicit Equations (linear – will be eliminated – do not introduce additional states): 0 = Lba (t ) − Lab (t ) 0 = Lcb (t ) − Lbc (t ) 0 = Lca (t ) − Lac (t ) 0 = L fa (t ) − Laf (t ) 0 = L fb (t ) − Lbf (t ) 0 = L fc (t ) − Lcf (t ) 0 = LDa (t ) − LaD (t ) 0 = LDb (t ) − LbD (t ) 0 = LDc (t ) − LcD (t ) 0 = LQa (t ) − LaQ (t ) 0 = LQb (t ) − LbQ (t ) 0 = LQc (t ) − LcQ (t ) 48 Note that L fQ = LQf = 0 and LDQ = LQD = 0 since the Q and D windings are perpendicular. In compact matrix notation the model is: i (t ) = A1v(t ) + A4 z (t ) dx(t ) 0= + B2 x(t ) + B3 y (t ) + B4 z (t ) dt 0 = C1v(t ) + C2 x(t ) + C3 y (t ) + C4 z (t ) 0 = q(v(t ), x(t ), y(t ), z (t )) B.3 Quasi-Steady State Model (Mixed Domain) The quadratized equations of the model are (separated into external, internal differential, internal linear and internal quadratic) as follows. Note that time dependence is explicitly shown for time domain quantities, but it’s not shown, though implied for phasor quantities. External equations: ~~ Ia = Ia ~~ Ib = Ib ~~ Ic = Ic ~ ~~~ I n = −I a − I b − I c i f = i f (t ) i fn = −i f (t ) Tm = Tm (t ) ~~~~ State: v (t ) = [Va Vb Vc Vn v f (t ) v fn (t ) Tm (t ) ] T Internal equations: Differential Equations: dλ f (t ) = e f (t ) dt d λ D (t ) = e D (t ) dt dλQ (t ) = e Q (t ) dt 49 d δ (t ) = ω m (t ) − ω s dt dω m (t ) 1 = Tacc (t ) dt J dc (t ) = y1 ( t ) dt ds(t ) = y 2 (t ) dt State: x(t ) = λ f (t ) λ D (t ) λQ (t ) δ (t ) ω m (t ) c(t ) s (t ) Linear equations: [ T ~ ~ jωΛ a = Ea ~ ~ jωΛ b = Eb ~ ~ jωΛ c = Ec 0 = Tacc (t ) − Te (t ) − Tm (t ) − T fw (t ) 0 = ω (t ) − p ω m (t ) 2 ~ ~~~ 0 = Ea + ra I a − Va + Vn ~ ~~~ 0 = Eb + rb I b − Vb + Vn ~ ~~~ 0 = Ec + rc I c − Vc + Vn 0 = e f (t ) + r f i f (t ) − v f (t ) + v fn (t ) 0 = eD (t ) + rD iD (t ) 0 = eQ (t ) + rQ iQ (t ) ~ State: y (t ) = Λ a [ ~ Λb ~ Λc Tacc (t ) ω (t ) ~ Ea ~ Eb ~ Ec e f (t ) e D (t ) eQ (t ) T Nonlinear Equations: 0 = ω m (t )Te (t ) − Pe (t ) ~ ~* ~ ~* ~ ~* 0 = Pe (t ) − Re E a I a − Eb I b − Ec I c − e f (t )i f (t ) − e D (t )i D (t ) − eQ (t )iQ (t ) { } 0 = y1 (t ) + s (t ) ⋅ (ω m (t ) − ω s ) 0 = y 2 (t ) − c(t ) ⋅ (ω m (t ) − ω s ) ~~ ~~ ~~~ ~ ~ ~ ~ ~ ~ ~ ~ ~ 0 = Λa − Ls Ia − F Laa * Ia + Ms Ib + F Lab * Ib + Ms Ic + F Lac * Ic − Laf i f (t) − LaDiD(t) − LaQiQ(t) 1 1 1 ~~ ~~ ~~~ ~ ~ ~ ~ ~ ~ ~ ~ ~ 0 = Λb + Ms Ia + F Lba * Ia − Ls Ib − F Lbb * Ib + Ms Ic + F Lbc * Ic − Lbf if (t) − LbDiD(t) − LbQiQ(t) 1 1 1 ′ 0 = T wf ( t ) + D wf ω m ( t ) + D wf ω m ( t ) 2 {} {} {} {} {} {} 50 ~~ ~~ ~~~ ~ ~ ~ ~ ~ ~ ~ ~ ~ 0 = Λc + Ms Ia + F Lca * Ia + Ms Ib + F Lcb * Ib − Ls Ic − F Lcc * Ic − Lcf if (t) − LcDiD(t) − LcQiQ(t) 1 1 1 ~~ ~~ ~~ 0 = Λ f − F0 Lfa * Ia − F0 Lfb * Ib − F0 Lfc * Ic − Lf i f (t) − MRiD (t) − LfQiQ (t) ~~ ~~ ~~ 0 = ΛD − F0 LDa * Ia − F0 LDb * Ib − F0 LDc * Ic − MRi f (t) − LDiD (t) − LDQiQ (t) ~~ ~~ ~~ 0 = ΛQ −F LQa*Ia −F0 LQb*Ib −F LQc*Ic −LQfif (t) −LQDiD(t) −LQiQ(t) 0 0 { { { }{ }{ } }{ }{ } }{ }{ } y1 (t ) {} {} {} State: z (t ) = Te (t ) [ Pe (t ) Twf (t ) y2 (t ) ~ Ia ~ Ib ~ Ic i f (t ) iD (t ) iQ (t ) T In many cases of quasi-steady state analysis it makes sense for all the electrical transients to be neglected, since they are too fast compared to the slower mechanical transients. Therefore, in this case, all the electrical states become algebraic states, and all the corresponding differential equations become algebraic equations. In this case the field and damper winding quantities are constant DC values. Note that, for convenience, time dependence is not explicitly shown, though implied, except for the states associated with the dynamical equations. External equations: ~~ Ia = Ia ~~ Ib = Ib ~~ Ic = Ic ~ ~~~ I n = −I a − I b − I c If = If I fn = − I f Tm = Tm (t ) ~~~~ State: v (t ) = [Va Vb Vc Vn V f V fn Tm (t ) T Internal equations: Differential Equations: 0 = Ef 0 = ED 0 = EQ d δ (t ) = ω m (t ) − ω s dt dωm (t ) 1 = Tacc (t ) dt J dc (t ) = y1 (t ) dt 51 ds(t ) = y 2 (t ) dt State: x(t ) = f Linear equations: ~ ~ jωΛ a = Ea ~ ~ jωΛ b = Eb ~ ~ jωΛ c = Ec [ ΛD Λ Q δ (t ) ωm (t ) c(t ) s (t ) T 0 = Tacc (t ) − Te (t ) − Tm (t ) − T fw (t ) 0 = ω (t ) − p ω m (t ) 2 ~ ~~~ 0 = Ea + ra I a − Va + Vn ~ ~~~ 0 = Eb + rb I b − Vb + Vn ~ ~~~ 0 = Ec + rc I c − Vc + Vn 0 = E f + rf I f − V f + V fn 0 = ED + rD I D 0 = EQ + rQ I Q ~ State: y (t ) = Λ a [ ~ Λb ~ Λc Tacc (t ) ω (t ) ~ Ea ~ Eb ~ Ec Ef ED EQ T Nonlinear Equations: 0 = ω m (t )Te (t ) − Pe (t ) ~ ~ * ~ ~* ~ ~* 0 = Pe (t ) − Re Ea I a − Eb I b − Ec I c − E f I f − ED I D − EQ I Q { } 0 = y1 (t ) + s (t ) ⋅ (ω m (t ) − ω s ) 0 = y 2 (t ) − c(t ) ⋅ (ω m (t ) − ω s ) ~~ ~~ ~~~ ~ ~ ~ ~ ~ ~ ~ ~ ~ 0 = Λa − Ls Ia − F Laa * Ia + Ms Ib + F Lab * Ib + Ms Ic + F Lac * Ic − Laf I f − LaDID − LaQIQ 1 1 1 ~~ ~~ ~~~ ~ ~ ~ ~ ~ ~ ~ ~ ~ 0 = Λb + Ms Ia + F Lba * Ia − Ls Ib − F Lbb * Ib + Ms Ic + F Lbc * Ic − Lbf I f (t) − LbDID − LbQIQ 1 1 1 ~~ ~~ ~~~ ~ ~ ~ ~ ~ ~ ~ ~ ~ 0 = Λc + Ms Ia + F Lca * Ia + Ms Ib + F Lcb * Ib − Ls Ic − F Lcc * Ic − Lcf I f − LcDID − LcQIQ 1 1 1 ~~ ~~ ~~ 0 = Λ f − F0 Lfa * Ia − F0 Lfb * Ib − F0 Lfc * Ic − Lf I f − MR I D − LfQIQ ~~ ~~ ~~ 0 = ΛD − F0 LDa * Ia − F0 LDb * Ib − F0 LDc * Ic − MR I f (t) − LDID − LDQIQ ~~ ~~ ~~ 0 = ΛQ −F0 LQa*Ia −F0 LQb*Ib −F0 LQc*Ic −LQfif (t) −LQDID −LQIQ ′ 0 = T wf ( t ) + D wf ω m ( t ) + D wf ω m ( t ) 2 { { { {} {} {} }{ }{ } }{ }{ } }{ }{ } {} {} {} {} {} {} State: 52 z (t ) = Te (t ) [ Pe (t ) Twf (t ) y1 (t ) y2 (t ) ~ Ia ~ Ib ~ Ic If ID IQ T Note that L fQ = LQf = 0 perpendicular. ζ ( t ) = A1v ( t ) + A4 y ( t ) 0= dx(t ) + B2 x(t ) + B3 y (t ) + B4 z (t ) dt 0 = C1v(t ) + C2 x(t ) + C3 y(t ) + C4 z(t ) and LDQ = LQD = 0 since the Q and D windings are 0 = q (v (t ), x (t ), y (t ), z (t )) Note also that the single “~” notation indicates phasors at the fundamental frequency, while the double “~” notation indicates phasors at twice the fundamental frequency. No “~” indicates DC quantities. The “*” symbol indicates convolution in the frequency domain. Note that the results of these convolutions include components at several frequencies. More specifically, the stator fluxes include components at fundamental and third harmonic, while the field and damper winding fluxes include DC components and second harmonic components. These higher than the fundamental components are neglected. Thus, F1 indicates the fundamental frequency component of the involved expression, while F0 indicates the DC component. Note also that the time varying self and mutual inductances and their phasor representations are defined as follows: Laa (t ) = Ls + Lm cos(2θ (t ) ) Lbb (t ) = Ls + Lm cos 2θ (t ) − 2π Lcc Lab s m Lbc (t ) = Lcb (t ) = − M S − Lm cos(2θ (t ) − π ) ( ) 3 (t ) = L + L cos(2θ (t ) + 2π ) 3 (t ) = L (t ) = −M − L cos(2θ (t ) + π ) 3 ba S m Lca (t ) = Lac (t ) = −M S − Lm cos 2θ (t ) − π Laf ( t ) = L fa ( t ) = M F cos θ( t ) ( 3 ) cos θ (t ) + 2π Lbf ( t ) = L fb ( t ) = M F cos θ( t ) − 2π Lcf (t ) = L fc (t ) = M F LaD ( t ) = LDa ( t ) = M D cos θ( t ) ( ) 3 cos(θ (t ) − 4π ) = M 3 ( 3 F ( 3 ) LbD (t ) = LDb (t ) = M D cos θ(t ) − 2π ) 53 LcD (t ) = LDc (t ) = M D cos θ (t ) + 2π LaQ (t ) = LQa (t ) = M Q LbQ (t ) = LQb (t ) = M Q LcQ (t ) = LQc (t ) = M Q cos(θ (t ) − π + 2π ) 2 3 and therefore ~ ~ Laa = Lm e j 0 2π ~ −j ~ Lbb = Lme 3 2π ~ j ~ Lcc = Lme 3 π ~ ~ j ~ ~ Lab = Lba = Lme 3 ~ ~ ~ ~ Lbc = Lcb = Lm e − jπ π ~ ~ −j ~ ~ Lca = Lac = Lme 3 ~ ~ Laf = L fa = M F e j 0 −j ~ ~ Lbf = L fb = M F e 2π 3 ( ) 3 cos(θ (t ) − π ) 2 cos(θ (t ) − π − 2π ) 2 3 j ~ ~ Lcf = L fc = M F e 3 ~ ~ LaD = LDa = M D e j 0 2π −j ~ ~ LbD = LDb = M D e j ~ ~ LcD = LDc = M F e 2π 3 2π 3 −j ~ ~ LaQ = LQa = M Q e 2 − j⎜ + ~ ~ LbQ = LQb = M Q e ⎝ 2 − j⎜ − ~ ~ LcQ = LQc = M Q e ⎝ 2 ⎛π ⎛π 2π ⎞ ⎟ 3⎠ 2π ⎞ ⎟ 3⎠ π B.4 Frequency Domain Steady State Model The quadratized equations of the model are (separated into external, internal differential, internal linear and internal quadratic) as follows: 54 External equations: ~~ Ia = Ia ~~ Ib = Ib ~~ Ic = Ic ~ ~~~ I n = −I a − I b − I c If = If I fn = − I f Tm = Tacc − Te − T fw ~~~~ State: v = [Va Vb Vc Vn V f V fn Tm ] T Internal equations: Linear equations: ~ ~ jωΛ a = Ea ~ ~ jωΛ b = Eb ~ ~ jωΛ c = Ec 0 = Ef 0 = ED 0 = EQ ωm = ω s 0 = Tacc p ωm 2 ~ ~~~ 0 = Ea + ra I a − Va + Vn ~ ~~~ 0 = Eb + rb I b − Vb + Vn ~ ~~~ 0 = Ec + rc I c − Vc + Vn 0=ω− 0 = E f + r f I f − V f + V fn 0 = E D + rD I D 0 = EQ + rQ I Q ~ State: y = Λ a [ ~ Λb ~ Λc Λf ΛD ΛQ ωm Tacc ω ~ Ea ~ Eb ~ Ec Ef ED EQ T Nonlinear Equations: 0 = ω mTe − Pe ′2 0 = T wf + D wf ω m + D wf ω m 55 ~ ~* ~ ~* ~ ~* 0 = Pe − Re E a I a − Eb I b − E c I c − E f I f − E D I D − EQ I Q ~~ ~~ ~~~ ~ ~ ~ ~ ~ ~ ~ ~ ~ 0 = Λa − Ls Ia − F Laa * Ia + Ms Ib + F Lab * Ib + Ms Ic + F Lac * Ic − Laf I f − LaDID − LaQIQ 1 1 1 ~~ ~~ ~~~ ~ ~ ~ ~ ~ ~ ~ ~ ~ 0 = Λb + Ms Ia + F Lba * Ia − Ls Ib − F Lbb * Ib + Ms Ic + F Lbc * Ic − Lbf I f − LbDID − LbQIQ 1 1 1 ~~ ~~ ~~~ ~ ~ ~ ~ ~ ~ ~ ~ ~ 0 = Λc + Ms Ia + F Lca * Ia + Ms Ib + F Lcb * Ib − Ls Ic − F Lcc * Ic − Lcf I f − LcDID − LcQIQ 1 1 1 ~~ ~~ ~~ 0 = Λ f − F0 Lfa * Ia − F0 Lfb * Ib − F0 Lfc * Ic − Lf I f − MR I D − LfQIQ ~~ ~~ ~~ 0 = ΛD − F0 LDa * Ia − F0 LDb * Ib − F0 LDc * Ic − MR I f − LDID − LDQIQ ~~ ~~ ~~ 0 = ΛQ −F0 LQa*Ia −F0 LQb*Ib −F0 LQc*Ic −LQf I f −LQDID −LQIQ { { { { {} {} {} }{ }{ } }{ }{ } }{ }{ } ~ Ia ~ Ib ~ Ic If ID {} {} {} } {} {} {} State: z = Te [ Twf Pe IQ T Note that L fQ = LQf = 0 perpendicular. and LDQ = LQD = 0 since the Q and D windings are Note also that the single “~” notation indicates phasors at the fundamental frequency, while the double “~” notation indicates phasors at twice the fundamental frequency. The quantities without the sign “~” indicate DC quantities. The “*” symbol indicates convolution in the frequency domain. Note that the results of these convolutions include components at several frequencies. More specifically, the stator fluxes include components at fundamental and third harmonic, while the field and damper winding fluxes include DC components and second harmonic components. These higher than the fundamental components are neglected. Thus, F1 indicates the fundamental frequency component of the involved expression, while F0 indicates the DC component. Note also that the time varying self and mutual inductances and their phasor representations are defined as follows: Laa (t ) = Ls + Lm cos(2θ (t ) ) Lbb (t ) = Ls + Lm cos 2θ (t ) − 2π Lcc Lab s m Lbc (t ) = Lcb (t ) = − M S − Lm cos(2θ (t ) − π ) ( ) 3 (t ) = L + L cos(2θ (t ) + 2π ) 3 (t ) = L (t ) = −M − L cos(2θ (t ) + π ) 3 ba S m Lca (t ) = Lac (t ) = −M S − Lm cos 2θ (t ) − π Laf ( t ) = L fa ( t ) = M F cos θ( t ) ( 3 ) Lbf ( t ) = L fb ( t ) = M F cos θ( t ) − 2π ( 3 ) 56 Lcf (t ) = L fc (t ) = M F cos θ (t ) − 4π LaD ( t ) = LDa ( t ) = M D cos θ( t ) ( 3 )= M ) F cos θ (t ) + 2π ( 3 ) LbD (t ) = LDb (t ) = M D cos θ(t ) − 2π LcD (t ) = LDc (t ) = M D cos θ (t ) + 2π LaQ (t ) = LQa (t ) = M Q LbQ (t ) = LQb (t ) = M Q LcQ (t ) = LQc (t ) = M Q cos(θ (t ) − π + 2π ) 2 3 and therefore ~ ~ Laa = Lm e j 0 2π ~ −j ~ Lbb = Lme 3 2π ~ j ~ Lcc = Lme 3 π ~ ~ j ~ ~ Lab = Lba = Lme 3 ~ ~ ~ ~ Lbc = Lcb = Lm e − jπ π ~ ~ −j ~ ~ Lca = Lac = Lme 3 ~ ~ Laf = L fa = M F e j 0 −j ~ ~ Lbf = L fb = M F e 2π 3 ( ) 3 cos(θ (t ) − π ) 2 cos(θ (t ) − π − 2π ) 2 3 ( 3 j ~ ~ Lcf = L fc = M F e 3 ~ ~ LaD = LDa = M D e j 0 2π −j ~ ~ LbD = LDb = M D e j ~ ~ LcD = LDc = M F e 2π 3 2π 3 −j ~ ~ LaQ = LQa = M Q e 2 − j⎜ + ~ ~ LbQ = LQb = M Q e ⎝ 2 − j⎜ − ~ ~ LcQ = LQc = M Q e ⎝ 2 ⎛π ⎛π 2π ⎞ ⎟ 3⎠ 2π ⎞ ⎟ 3⎠ π 57 The presented steady state model assumes that, at steady state, an excitation signal (voltage or current) is applied to the field winding from an exciter model (steady state model) specifying the field winding current and voltage and thus the voltage output of the unit and a steady state torque signal is applied by a prime mover (steady state model), specifying the produced active power by the unit. Therefore, the unit is operating, in fact, at a PV mode, with P and V defined implicitly and externally by another model. A more conventional and probably intuitive approach would be to specify the three possible operating modes of the generator (a) PV mode, (b) PQ mode, and (c) slack mode with additional equations. These equations will not introduce additional states in this model, but will introduce new states, (in fact make inputs become states) in the prime mover and exciter models respectively. Therefore, the generator, exciter and prime mover models will not be square models, but the whole generating unit system will be. (a) PV mode: ~~ ~~ ~~ 0 = Pspec + Re Va I a* + Vc I c* + Vc I c* ~ ~ 0 = Vspec − V1 , where V1 is the positive sequence component of the terminal voltage ~~~ (computed directly as a weighted sum of Va , Vb , Vc ). { } Alternatively, if only a specific phase if monitored by the AVR, e.g. phase A, the second equation will become: ~ 0 = Vspec − Va Similarly, if the voltage between two phases is monitored, e.g., then the second equation will become: ~~ 0 = Vspec − Va − Vb (b) PQ mode: ~~ ~~ ~~ 0 = Pspec + Re Va I a* + Vc I c* + Vc I c* ~~ ~~ ~~ 0 = Qspec + Im Va I a* + Vc I c* + Vc I c* { { } } (c) Slack mode: ~ 0 = Vspec − V1 ~ 0 = Im V1 ~ where V1 is the positive sequence component of the terminal voltage (computed directly ~~~ as a weighted sum of Va , Vb , Vc ). {} Alternatively, if only a specific phase if monitored by the AVR, e.g. phase A, the equations will become: 58 ~ 0 = Vspec − Va ~ 0 = Im Va Similarly, if the voltage between two phases is monitored, e.g., then the equations will become: ~~ 0 = Vspec − Va − Vb ~~ 0 = Im Va − Vb {} { } 59 ...
View Full Document

This note was uploaded on 01/29/2011 for the course ENGR 52 taught by Professor Mcmillan during the Spring '10 term at Baylor Med.

Ask a homework question - tutors are online