Observer-Based Anomaly Detection of Synchronous Generators for Power Systems Monitoring

This paper proposes a rigorous anomaly detection scheme, developed to spot power system operational changes which are inconsistent with the models used by operators. This novel technique relies on a state observer, with guaranteed estimation error convergence, suitable to be implemented in real time, and it has been developed to fully address this important issue in power systems. The proposed method is fitted to the highly nonlinear characteristics of the network, with the states of the nonlinear generator model being estimated by means of a linear time-varying estimation scheme. Given the reliance of the existing dynamic security assessment tools in industry on nominal power system models, the suggested methodology addresses cases when there is deviation from assumed system dynamics, enhancing operators’ awareness of system operation. It is based on a decision scheme relying on analytical computation of thresholds, not involving empirical criteria which are likely to introduce inaccurate outcomes. Since false-alarms are guaranteed to be absent, the proposed technique turns out to be very useful for system monitoring and control. The effectiveness of the anomaly detection algorithm is shown through detailed realistic case studies in two power system models.

Vector of maximum measured input noise values. υ Vector of maximum measurement noise values. n Transformer off-nominal ratio. w Vector of maximum process disturbance values. δ Rotor angle in rad. Vector of state estimation errors. ν Column vector of measured input noise variables. ω, ω B Rotor speed in p.u. and its base value in rad/s. φ I y Measured stator current phase. Φ Discrete-time form of system transition matrix. ψ 1d Subtransient emf due to d-axis damper coil in p.u. ψ 2q Subtransient emf due to q-axis damper coil in p.u. I n Identity n × n matrix. θ HV bus voltage phase in rad. Ratio (X d − X ls ) / (X d − X ls ).

K d2
Ratio (X d − X d ) / (X d − X ls ). K q 1 Ratio X q − X ls / X q − X ls . K q 2 Ratio X q − X q / X q − X ls . M Inertia constant in p.u. P y , υ P Measured active power value and its associated noise. Q y , υ Q Measured reactive power value and its associated noise. M ODERN power networks are undergoing considerable operational and structural modifications, resulting from significant technology advances with respect to forms of electric power generation and communications infrastructure, let alone the energy market liberalisation [1], [2]. Power system complexity increases, given the current trends. Specifically, advanced control systems associated with newly integrated system components may give rise to unpredictable system behaviour of nonlinear nature. Moreover, aging equipment, given the longstanding power network operation, is prone to failures, whereas, lack of investments in transmission systems may lead to stressed system operation [2]. In this context, it is widely recognised that wide area monitoring systems (WAMS) contribute to operators' knowledge of system's operational status, and dynamic state estimation (DSE) may prove to be a very useful tool [3].
Real-time DSE is model based, thus, good knowledge of the power system dynamic model is essential to obtain accurate results. In this context, given power networks' high complexity, methods detecting deviation from routine conditions or 'nominal' models play a significant role in enhancing operators' awareness of the system operation. Such anomaly detection techniques have been reported in power systems literature, engaging Kalman filter variants, and relying on thresholds, the values of which are typically based on empirical criteria [4]- [9]. In this respect, Kalman filter utilization requires knowledge and assumptions regarding noise probability distributions (e.g. Gaussianity), while divergence issues are likely to arise when there is a high level of mismatch between the assumed estimation system model and the real one [10].
On the other hand, observer-based anomaly detection techniques constitute a significant class of methods identifying deviations from nominal models, being popular in various research areas [11]- [18]. Observers are developed based on the principle of guaranteed estimation error convergence (whereas Kalman filter based estimators primarily rely on the trace minimization of the state error covariance) [19], [20]. There have been several observer based research studies, involving power system models requiring many assumptions and simplifications [21], [22]; in cases where more advanced power system models are utilized, they are linearised with respect to one specific operating point, with the observer's filtering matrix term being optimized with reference to this particular point [23]- [26], restricting the applicability of the observer in a continuously changing environment like operation of a power network.
To address these issues, the work conducted in [27] is extended to consider nonlinear output equations, leading to the development of a discrete-time observer, based on time-varying linearisation with respect to the estimate at every time step, in a similar approach as in Extended Kalman filtering (EKF) [10]. The filtering term, guaranteeing the observer's convergence, changes at every time step, and its computation is simple and suitable for real-time implementations. This research effort leads to the following main contributions: r to propose a novel observer in the context of realistic power system models; r to establish an observer-based anomaly detection technique for models with nonlinear output equations; r to formulate a systematic approach for threshold calculation, guaranteeing the absence of false alarms, unlike the case where empirical criteria are used. The paper is organised as follows: In the next section, the synchronous generator models dealt with in this research are presented and thoroughly analysed. Section III includes the analysis of the observer-based anomaly detection algorithm. In Section IV, the proposed methodology is applied to two study systems; the 9-bus, 3-machine system used in [28] and the IEEE benchmark 68-bus, 16-machine system, corresponding to the New England (NETS) and New York (NYPS) power systems, along with three neighbouring areas [29]. Section V gives some concluding remarks.

A. Model Development -State Equations
Synchronous generator modelling is one of the basic requirements of power system dynamic modelling. Depending on the targets of each research and the amount of information available about the modelling detail, various models have been reported in literature [28], [30]. In the context of multi-machine systems, decentralisation has gained popularity, since it reduces the computational burden and enables the estimation procedure to be conducted on a local basis [31]- [35]. The decentralisation procedure relies on system partitioning (in terms of the estimation calculations) and the use of measured inputs, which are measurements on the 'boundary' of the assumed subsystem [31], [32]. This partitioning based approach is usually called 'playback process', which is also useful in generator model validation procedures [36]. Here, the 'boundary' is the high voltage (HV) side of the transformer connected to the terminal bus of the respective generator (henceforth termed as 'HV bus'), to enable local anomaly detection based on measurements at the bus which corresponds to this side of the transformer, as in [36]. The assumed estimation model boundary is depicted in Fig. 1. In several power system studies, the noise component of the measured inputs has been taken into account in the context of the analysis [32], [33], whereas, in some other studies, these inputs are considered as noise-free [34], [35], [37]. Furthermore, the consideration of different measured inputs leads to different decentralised synchronous generator models, with diverse models having been utilized for this purpose [32]- [35], [37]. Two synchronous generator models have been used in the context of the case studies here, the transient and the subtransient models [28], [38], [39]. Both models are valid for transient conditions, with the subtransient being more detailed and capturing phenomena taking place immediately after an event occurrence [30]. The subtransient synchronous generator model used here, including the respective transformer model elements, is described by the following equations [28], [38]: and To obtain the discrete form of these equations, it can be as- This model is built based on the principle that d-axis leads q-axis. Similarly to the analysis in [33], the internal rotor angle (α) is used instead of the rotor angle (δ), since the utilization of the latter relies on the multi-machine context, requiring the knowledge of the global reference frame angle, which is impossible when estimation is performed based on local information and measurements only. Here, the phase angles of all quantities are defined with respect to the HV bus voltage phasor. With regard to this, the internal rotor angle (α) replaces the rotor angle (δ) in this decentralised context, as carried out in [33], defined as α = δ − θ, with (1) characterizing its dynamics [33]. The rate of change of the HV bus voltage phase is approximated by the equation below (divided by ω B , to obtain the p.u. value) [28]: Analogously to the procedure described in [33], the measured inputs used are the HV bus voltage magnitude (V y ) and the rate of change of its phase (f y ). Thus, in the context of the synchronous generator subtransient model, the state vector has the following form: Also, the noisy input vector is: whereas: where u, ν correspond to the noise-free and noise components of the noisy input vector, respectively (i.e.ũ = u + ν). It can be noted that the first two elements of ν are equal to 0, since T m and V ref are considered to be perfectly known.
Further to the subtransient model, the transient model can be obtained by considering X q = X q and X d = X d , meaning that K q 1 = K d1 = 1 and K q 2 = K d2 = 0; therefore, in this case, the state vector does not include ψ 2q and ψ 1d ((5), (6) are omitted).

B. Output Equations
This work considers a decentralized generator model, with a Phasor Measurement Unit (PMU) at the HV side of the transformer which comes after the generator terminal bus. The current magnitude (I y ), its phase with reference to the HV bus voltage phasor (φ I y ), the active power output (P y ) and the reactive power output (Q y ) at the HV bus are the measurement outputs. These are governed by the following equations: where I q , I d are given by (9), whereas υ I , υ φ I are the measurement noise terms, associated with I y and φ I y , respectively, P y and Q y are the active and reactive power, respectively, measured at the HV substation of the transformer of the corresponding generator, and υ P , υ Q are their associated noise terms, respectively.
Frequency measurement (f sy s y ) has also been considered, since it is closely related to speed [40], and its p.u. value is: where υ f is the associated measurement noise. Therefore, the output measurement vector is: y = f sy s y I y φ I y P y Q y .

A. The Anomaly Detection Logic
The proposed anomaly detection method is inspired by the work conducted in [27], which has been significantly extended to consider the nonlinear output equations. (1)- (16) and (22)- (26) constitute the nonlinear synchronous generator state space model, which can be discretized and described by the following general discrete-time model: where f represents the discretized system state (1)-(16) at time step k, h corresponds to the discretized output measurement (22)- (26),w k is the process disturbance vector accounting for modelling uncertainty, whereasυ k is the output measurement noise vector.
In order for anomalies in the system to be detected, a model based decision scheme is proposed. A predictionŷ k of the output measurements y k is computed with reference to the nominal model (28). The prediction error r k = y k −ŷ k is compared with a suitably defined threshold to decide whether the current prediction error is just caused by the process and the output disturbances or also by the additional influence of an anomaly causing a significant discrepancy between the measured output and the one predicted via the nominal model.
The following time-varying observer is defined to compute the output prediction: with A k , B k , C k , D k ,b k ,c k defined as: wherex is the state estimation, and K a matrix collecting the time-varying observer parameters. Matrix K k is computed at each time step as in [41], guaranteeing the convergence of the state estimation error k = x k −x k (see the theoretical analysis illustrated in the Appendix B). The prediction error is compared component-by-component with a suitably defined threshold (see Section III-C). If, for at least one component, the residual crosses the corresponding threshold, this designates an anomaly (i.e., the measurements cannot be "explained" by the nominal dynamics). However, before the threshold analysis is presented, the system model (28) has to be formulated in an appropriate way with respect to the aforementioned observer equations, as illustrated in the following section.

B. Reformulation of the Nominal Model
From (1)- (16) and (22)-(26), it is clear that the synchronous generator's state and output equations are nonlinear. In order to allow dynamical analysis for the rigorous definition of a suitable anomaly detection threshold, the model (28) is reformulated as a linear time-varying model. For this purpose, the discrete-time state and output equations are linearised around the estimate at every time step (a similar approach is described in [10]), to obtain the prediction for the next time step. Thus, the previously presented synchronous generator models are rewritten according to the following discrete-time state space form: where w k and υ k also include the linearisation errors. In the context of the analysis conducted here, the process disturbance vector w lies within a range of values corresponding to nominal operation, and, as far as the measured input and output measurement noises (i.e. ν and υ, respectively) are concerned, they are characterized by maximum errors, as specified by standards, with which the measurement devices used have to comply, and in this case the IEEE Standard C37.118.1-2011 for phasor measurement units (PMUs) [42]. This means that: If any of the above inequalities is violated, this denotes deviation from nominal system operation and, thus, it is considered as an anomaly.

C. The Anomaly Detection Threshold
The anomaly detection threshold is designed in order to act as an upper bound for the output estimation error in the presence of disturbances. The threshold is analytically designed so to guarantee the following: i) to always be an upper bound for the residual; ii) absence of false-alarms; iii) stability of the threshold. To derive a suitable threshold, the output estimation error is analysed: Thanks to the result in Proposition B.1 (see Appendix B) guaranteeing the convergence of the estimator, and, from (31), the output estimation error (32) can be bounded at each time step component-wise using the triangular inequality: The state estimation error k = x k −x k is then analysed (see (38) in Appendix B). It can be written as: where Φ k,0 is the transition matrix from time 0 to time k: . Since the unforced system describing the dynamics of the estimation error (Appendix B - (35)) is exponentially asymptotically stable thanks to the result in Appendix B, in [43] it is proved that there exist a > 0 and λ ∈ [0, 1] so that: Therefore, the state estimation error can be bounded by: where¯ 0 is the initial estimation error bound, properly defined thanks to the knowledge of the system, andd l = |B l |ν +w + |K l | (|D l |ν +ῡ). As a consequence, we can define the anomaly detection threshold for the residual signal r: Thanks to the way it is designed, the following inequality is a necessary condition for the 'healthy' status of the monitored system: It is sufficient that at least one component of the residual signal r (i) crosses the corresponding thresholdr (i) for at least one time-step to state that the monitored system cannot be explained by the nominal dynamics (30), i.e., something has changed in the dynamics of the system from its modelled dynamics.

A. 9-Bus, 3-Machine System
The anomaly detection methodology has been applied to a 9-bus, 3-machine power system model, the details of which can be found in [28], and it is shown in Fig. 2. In this model, all synchronous machines are designed using their transient models, as previously presented, whereas, fast static exciters have been used for all the machines, described by (7), as opposed to the DC excitation systems which are utilized in [28], [39], to enable fast system response to contingencies [4]. With reference to (7), the parameter details of the excitation system used are K A = 200, and T A = 0.01 s, for all generators. Regarding the estimation procedure, IEEE Standards-compliant PMUs with reporting rate of 120 frames per second are used, to obtain measurements at the HV side of the transformer which comes after the terminal bus of the synchronous generator of interest. The maximum error of the measured voltage magnitude is 9 · 10 −3 p.u., and, the one corresponding to the voltage phase is 2 · 10 −3 rad, whereas the maximum errors for the measured active and reactive power values is 6 · 10 −3 p.u. This errors are set in accordance with the PMU related requirement of maximum 1% total vector error (TVE), dictated by the IEEE Standard C37.118.1-2011, blending together three possible sources of error for each phasor: phasor magnitude, angle and time synchronisation [42]. In the same context, the maximum frequency error is 0.005 Hz [42]. It has to be noted that the process disturbance is comprised by the noise coming from the measured inputs, and the estimation model is considered to be a highly accurate representation of the real one. Power system modelling is MATLAB/Simulink based, and all simulations last for 10 s. Three case studies have been considered: r Case Study 1C: The overexcitation limiter (OEL) of Gen.
1 gets activated at the time instant t = 3.34 s. The OEL used is of takeover type, designed according to the error signal substitution scheme, the design characteristics and the block diagram of which can be found in [4]. This OEL type is likely to lower system stability margins, therefore, it is important for its activation to be detected by the proposed anomaly identification scheme [4]. The measurements are obtained at bus 4. The results showing the residual and the threshold values for all measurement outputs of the case studies 1A, 1B, 1C, and Gens. 3, 2, 1, respectively, are illustrated in Figs. 3, 4 and 5, respectively. It can be clearly noticed that the proposed method is able to capture all simulated events, with at least one residual value exceeding the value of the corresponding threshold. Specifically, the events are detected at time instants 2.01 s, 2.01 s and 3.47 s, for case studies 1A, 1B, and 1C, respectively, denoting in this way success of the anomaly detection strategy.

B. IEEE Benchmark 68-Bus, 16-Machine System
The successful performance of the suggested anomaly detection scheme on the previous study system, incentivized its application to a larger, realistic power system. For this purpose, the IEEE benchmark 68-bus, 16-machine NETS-NYPS system has been used [29], illustrated in Fig. 6. The details of this system can be found in [44]. Here, the synchronous machines, which are not characterized by manual excitation (i.e. the first 12), are equipped with fast static exciters, described by (7), to enable fast system response to contingencies, in a similar approach as in the previous system, and to facilitate the application of OELs [4]. The exciter design details are the same as in the previous study system. Moreover, same measurement noise levels as in the previous test system related case studies have been considered. Three case studies have been examined:      studies 2A, 2B, and 2C, respectively. It has to be noted that, the proposed anomaly detection strategy is able to capture the OEL operation, without the need for excitation voltage measurements, as in [4], where the practice of OEL activation detection is primarily based on excitation voltage measurements.

V. CONCLUSION
An observer-based anomaly detection scheme has been presented, which is able to trace deviations from nominal power system operation. This methodology relies on a linear timevarying observer, addressing the need for an estimation methodology tailored to a nonlinear system like a power network, and it can be implemented in real time. The method makes use of measurements at the high voltage side of the transformer connected to the terminal bus of the respective generator only, facilitating the utilization of a local approach. The proposed algorithm has been tested on two power system models, a small 9-bus 3-machine system, and the large, realistic IEEE benchmark 68-bus 16-machine system, where the results validate the success of its implementation. It is important to highlight that the anomaly detection scheme is based on a rigorous threshold calculation, without requiring the consideration of noise probability distributions or empirical criteria, making this method highly practical. This technique, being effective in detecting the activation of limiting devices such as overexcitation limiters, may serve as a valuable tool for model updates in the context of power system dynamic security assessment algorithms. The suggested methodology can prove to be extremely beneficial to operators in the context of power system monitoring, given the model uncertainty characterizing modern power network operation.

APPENDIX A ANOMALY DETECTION ALGORITHM
The anomaly detection scheme is summarized in Alg. 1.
Algorithm 1: Anomaly detection. Setx 0 ,ŷ 0 ,¯ 0 k = 1 repeat Collect measurements and inputs y k , u k Compute state and output estimatesx k +1 andŷ k (29) Compute output error r k = y k −ŷ k Compute estimation error bound¯ k (35) Compute detection thresholdr k (36) until r The proposed observer's (29) convergence analysis requires the consideration of the state estimation error, whose dynamics can be described by the following linear time-varying model: In order to guarantee the convergence of the state estimation error (38), the result in [41] is used to define the matrix K k at each discrete time-step k as: where t is a positive constant integer, Φ k,k−t−1 is the transition matrix from time k − t − 1 to time k: with Φ k,k = I, and G k,k−t−1 is the observability Gramian, defined as: We have the following result: Proposition B.1: The state estimation error (38), with the filtering matrix K k defined as in (39), represents the dynamics of a Bounded Input Bounded Output stable system.
Proof: In [41] it is demonstrated that the proposed filtering matrix K k allows to guarantee the exponentially asymptotically stability of the unforced system k +1 = (A k − K k C k ) k under the assumption that the pair [A k , C k ] is uniform with respect to complete observability. Since the uncertainties are all bounded, then the time-varying system is Bounded Input Bounded Output stable.