A Data Driven Approach to Robust Event Detection in Smart Grids Based on Random Matrix Theory and Kalman Filtering ^†

Abstract

Increasing levels of complexity, due to growing volumes of renewable generation with an associated influx of power electronics, are placing increased demands on the reliable operation of modern power systems. Consequently, phasor measurement units (PMUs) are being rapidly deployed in order to further enhance situational awareness for power system operators. This paper presents a novel data-driven event detection approach based on random matrix theory (RMT) and Kalman filtering. A dynamic Kalman filtering technique is proposed to condition PMU data. Both simulated and real PMU data from the transmission system of Great Britain (GB) are utilized in order to validate the proposed event detection approach and the results show that the proposed approach is much more robust with regard to event detection when applied in practical situations.

1. Introduction

The development and increasing complexity of modern power systems is leading to growing concerns over power system operational security, with network management and analysis becoming increasingly challenging. Traditional monitoring devices, employed in supervisory control and data acquisition (SCADA) systems, have become inefficient for this challenge [1]. As a result, phasor measurement units (PMUs) are being deployed to measure phasors of bus voltages and currents, enabling real-time state monitoring in power systems [2,3]. Due to PMUs’ higher sampling rates and accurate time synchronization, a huge amount of synchronous data is now being collected from power systems globally [4,5], facilitating the detection of transient events.

At present, event detection is carried out following either a model-driven approach (supervised) or a data-driven approach (unsupervised). Model-driven approaches include neural networks [6,7], decision trees [8,9], and principal component analysis (PCA) [10,11,12,13,14]. However, for supervised learning, it is argued that if the training dataset is selected inappropriately or insufficiently, it will have a significantly negative effect on the performance of the approach when it is applied online. Therefore, much attention has been paid to the employment of unsupervised learning in this field [15,16,17,18], including wavelet-based approaches, energy function-based approaches, and frequency difference-based approaches.

More recently, as a data-driven approach, random matrix theory (RMT) was introduced and applied to early event detection in power systems [19]. However, PMU data are normally collected with significant noise levels and plenty of missing measurements, due to communication issues between PMUs and phasor data concentrators (PDCs) [20,21] and further issues in the measurement processes. When the above approach is faced with significantly noisy and missing data, it will be unable to accurately detect abnormal events.

A large amount of research has been presented in the field of event detection for power systems. In [7], Vasilic and Kezunovic presented an approach for transmission line fault classification building on an adaptive resonance theory (ART) neural network and fuzzy logic. This work applies a fuzzy decision rule to the output of the ART neural network to improve selectivity for a variety of real events. In [9], Adhikari et al. presented a method to effectively classify traditional power system contingencies and cyber-attacks in real time. It combines Hoeffding adaptive trees with a drift detection method and adaptive windowing to update the model continuously.

In recent years, PCA has been widely applied for early event detection in power systems. For example, Xie et al. [10] analyzed the dimensionality of PMU data under both normal and abnormal conditions and employed PCA to detect abnormal events based on the change of core subspaces of PMU data at the occurrence of an anomaly. Meanwhile, this model has an adaptive training mechanism that updates the model when an event happens. A similar work was carried out in [11] that uses a sliding window PCA to detect and classify multiple events, including islanding, loss of load, and loss of generation. However, different from the work presented in [10], this method uses a sliding window to update the PCA model in real time to adapt to the time varying behaviors of power systems. Another similar event detection method was developed by Guo et al. [12], which focuses on islanding detection for distributed generation systems. In this work, a recursive PCA algorithm is employed to analyze the time-varying behaviors of power systems for the purpose of reducing false alarms. Although a wide range of machine learning techniques have been applied to event detection and classification, these methods are mainly based on model-driven supervised learning. As mentioned above, when the sample space is chosen improperly or insufficiently, an ill-trained model may be obtained.

In order to address the limitations of supervised learning techniques, a few works have considered unsupervised learning techniques. For instance, Kim et al. designed a wavelet-based detection method using PMU data [15]. The key idea is to monitor the energy of coefficients within a sliding window in real time and to determine whether the energy value exceeds the threshold. Apart from this work, another data-driven analytics method, which relies on the rules created by PMU data, was proposed for power system fault detection in [18]. However, it should be noted that these unsupervised approaches need to pre-define an appropriate threshold in detection of events, which is a challenging issue. The RMT-based event detection approach presented in [19] does not need to pre-define a threshold, but this approach does not take into account the quality aspects of PMU data in dealing with noisy and missing PMU measurements.

This paper presents an event detection approach that improves the performance of the work presented in [19], making it much more robust when dealing with corrupted and significantly noisy samples of PMU data. To be specific, a robust event detection approach is developed by combining RMT and Kalman filtering. Furthermore, a dynamic Kalman filter is developed to condition PMU data. Finally, the proposed event detection approach is tested and validated on the standard IEEE 118-bus system and the transmission system of Great Britain (GB). The results show that the proposed approach is much more robust in practical situations that include significant levels of noisy or missing PMU measurements.

The remainder of this paper is organized as follows. Section 2 briefly introduces random matrix theory and formulates a dynamic quadratic prediction model based Kalman filtering technique. Section 3 combines RMT with Kalman filtering for event detection. Section 4 evaluates the performance of the proposed event detection approach on both simulated and real PMU data. Finally, Section 5 concludes the paper and points out some future work.

2. Random Matrix Theory and Kalman Filtering

Random matrix theory is first introduced briefly and then a dynamic quadratic prediction model based Kalman filtering technique is formulated.

2.1. Random Matrix Theory

RMT has proven to be effective for power system analysis, and a wide range of laws and theorems, such as Ring Law and Marchenko–Pastur Law, are included in RMT, which have been applied in power systems for different purposes [19,22,23]. In this paper, for the purpose of event detection, Ring Law, as the theoretical basis of data processing of event detection, and mean spectral radius (MSR), as an event indicator, will be briefly introduced.

2.1.1. Ring Law

Let ${\mathbf{X}}_{\mathrm{s}}\in {\mathbf{C}}^{N\times T}$ be a standard non-Hermitian random matrix, whose entries are independent and identically distributed (i.i.d.) variables with

$$\mu \left({\mathbf{x}}_{\mathrm{s},i}\right)=0,{\sigma }^2\left({\mathbf{x}}_{\mathrm{s},i}\right)=1\hspace{1em}\left(i=1,2,\dots ,N\right)$$

(1)

where ${\mathbf{x}}_{\mathrm{s},i}$ is the ith row vector of ${\mathbf{X}}_{\mathrm{s}}$. For L standard non-Hermitian random matrices ${\mathbf{X}}_{\mathrm{s},i}\left(i=1,2,\dots ,L\right)$, a matrix product is defined as

$$\mathbf{Z}={\displaystyle \prod }_{i=1}^L{\mathbf{X}}_{\mathrm{u},i}$$

(2)

where ${\mathbf{X}}_{\mathrm{u},i}$ is the singular value equivalent of ${\mathbf{X}}_{\mathrm{s},i}$. The matrix product $\mathbf{Z}$ can be transformed to the standard matrix product ${\mathbf{Z}}_{\mathrm{s}}$, whose ${\sigma }^2\left({\mathbf{z}}_{\mathrm{s},i}\right)=1/N$ in each row. Thus, the empirical spectral density (ESD) of ${\mathbf{Z}}_{\mathrm{s}}$ converges almost surely to the limit given by

$$f\left({\lambda }_{{\mathbf{Z}}_{\mathrm{s}}}\right)=\{\begin{array}{ll}\frac{1}{\pi cL}{\left|{\lambda }_{{\mathbf{Z}}_{\mathrm{s}}}\right|}^{\frac{2}{L}-2}& {(1-c)}^{\frac{L}{2}}\le \left|{\lambda }_{{\mathbf{Z}}_{\mathrm{s}}}\right|\le 1\\ 0& \mathrm{otherwise}\end{array}$$

(3)

As $N,T\to \infty$ with a constant ratio $c=N/T\in \left(0,1\right]$, N is the number of variables, and T is the time length.

2.1.2. MSR

Linear eigenvalue statistics indicate the statistical features of random matrices. A linear eigenvalue statistic of a random matrix X is defined as

$${\mathrm{S}}_n\left(\varphi \right)={\displaystyle \sum }_{i=1}^n\varphi \left({\lambda }_i\right)$$

(4)

where ${\lambda }_i\left(i=1,2,\dots ,n\right)$ are eigenvalues of X, and $\varphi (\cdot )$ is a test function.

MSR is used to indicate the eigenvalue distribution of a random matrix. For the standard matrix product ${\mathbf{Z}}_{\mathrm{s}}$ (as mentioned above in Ring Law), MSR is formulated as

$$\mathrm{MSR}=\frac{1}{N}{\displaystyle \sum }_{i=1}^N\left|{\lambda }_{{\mathbf{Z}}_{\mathrm{s}},i}\right|$$

(5)

where ${\lambda }_{{\mathbf{Z}}_{\mathrm{s}},i}\left(i=1,2,\dots ,N\right)$ are eigenvalues of ${\mathbf{Z}}_{\mathrm{s}}$, and $\left|{\lambda }_{{\mathbf{Z}}_{\mathrm{s}},i}\right|$ is the radius of ${\lambda }_{{\mathbf{Z}}_{\mathrm{s}},i}$ on the complex plane.

2.1.3. Data Processing of Ring Law

Within a raw data source Ω, a raw random matrix $\mathbf{X}\in {\mathbf{C}}^{N\times T}$ can be formed through a sliding window. Then, X can be transformed to a standard non-Hermitian matrix ${\mathbf{X}}_{\mathrm{s}}\in {\mathbf{C}}^{N\times T}$ row by row as follows:

$${\mathbf{x}}_{\mathrm{s},i}=\left({\mathbf{x}}_i-\mu \left({\mathbf{x}}_i\right)\right)\frac{\sigma \left({\mathbf{x}}_{\mathrm{s},i}\right)}{\sigma \left({\mathbf{x}}_i\right)}+\mu \left({\mathbf{x}}_{\mathrm{s},i}\right)\hspace{1em}\left(i=1,2,\dots ,N\right)$$

(6)

where ${\mathbf{x}}_i$ and ${\mathbf{x}}_{\mathrm{s},i}$ are the ith row vectors of X and ${\mathbf{X}}_{\mathrm{s}}$, respectively, and $\mu \left({\mathbf{x}}_{\mathrm{s},i}\right)=0$, ${\sigma }^2\left({\mathbf{x}}_{\mathrm{s},i}\right)=1$.

Afterward, the matrix ${\mathbf{X}}_{\mathrm{u}}\in {\mathbf{C}}^{N\times N}$ is introduced as the singular value equivalent of ${\mathbf{X}}_{\mathrm{s}}$ by

$${\mathbf{X}}_{\mathrm{u}}=\sqrt{{\mathbf{X}}_{\mathrm{s}}{\mathbf{X}}_{\mathrm{s}}{}^{\mathrm{T}}}\mathbf{U}$$

(7)

where ${\mathbf{X}}_{\mathrm{s}}{}^{\mathrm{T}}$ is the transpose of ${\mathbf{X}}_{\mathrm{s}}$, and $\mathbf{U}\in {\mathbf{C}}^{N\times N}$ is a Haar unitary matrix, ${\mathbf{X}}_{\mathrm{u}}{\mathbf{X}}_{\mathrm{u}}{}^{\mathrm{T}}\equiv {\mathbf{X}}_{\mathrm{s}}{\mathbf{X}}_{\mathrm{s}}{}^{\mathrm{T}}$.

However, for a series of arbitrary non-Hermitian random matrices ${\mathbf{X}}_i\left(i=1,2,\dots ,L\right)$ from the raw data source Ω, the matrix product $\mathbf{Z}={\displaystyle \prod }_{i=1}^L{\mathbf{X}}_{\mathrm{u},i}\in {\mathbf{C}}^{N\times N}$ is obtained. Then, it is converted to the standard matrix product ${\mathbf{Z}}_{\mathrm{s}}\in {\mathbf{C}}^{N\times N}$ through the following formula:

$${\mathbf{z}}_{\mathrm{s},i}=\frac{{\mathbf{z}}_i}{\sqrt{N}\sigma \left({\mathbf{z}}_i\right)}\hspace{1em}\left(i=1,2,\dots ,N\right)$$

(8)

where ${\mathbf{z}}_i$ and ${\mathbf{z}}_{\mathrm{s},i}$ are the ith row vectors of Z and ${\mathbf{Z}}_{\mathrm{s}}$, respectively.

2.2. Dynamic Kalman Filtering Technique

As a state estimation technique, Kalman filtering is a popular method for conditioning PMU data in power systems [24,25,26]. Jones et al. [27] proposed a quadratic prediction model-based Kalman filtering method to cleanse and recover PMU data easily. However, within this method, the measurement noise covariance matrix R is constant, thus leading to inaccurate PMU data conditioning in some situations. As a result, this paper makes an improvement by dynamically adapting the measurement noise covariance matrix R to the conditioning process.

Following [27], the classic model of Kalman filtering is defined as follows:

$$\{\begin{array}{l}\mathbf{x}\left(t+1\right)=\mathsf{\Phi }\left(t+1,t\right)\mathbf{x}\left(t\right)+\mathsf{\Gamma }\left(t+1,t\right)\mathbf{w}\left(t\right)\\ \mathbf{z}\left(t+1\right)=\mathbf{H}\left(t+1\right)\mathbf{x}\left(t+1\right)+\mathbf{v}\left(t+1\right)\end{array}$$

(9)

For the purpose of optimization, the classic model of Kalman filtering can be expressed in the recursive way:

$$\{\begin{array}{l}\widehat{\mathbf{x}}(t+1|t+1)=\mathsf{\Phi }\left(t+1,t\right)\widehat{\mathbf{x}}(t|t)\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}+\mathbf{K}\left(t+1\right)(\mathbf{z}\left(t+1\right)-\widehat{\mathbf{z}}(t+1|t))\\ \widehat{\mathbf{z}}(t+1|t)=\mathbf{H}\left(t+1\right)\widehat{\mathbf{x}}(t+1|t)\end{array}$$

(10)

According to [28], the quadratic prediction model is formulated by

$$\widehat{x}(t+1|t)=3\widehat{x}(t|t)-3\widehat{x}(t-1|t-1)+\widehat{x}(t-2|t-2)$$

(11)

In addition, each state vector is formed by a sliding window, which contains three successive snapshots of the system states and moves forward only one snapshot at a time. Thus, in order to estimate the next state, the state vectors $\widehat{\mathbf{x}}(t|t)$ and $\widehat{\mathbf{x}}(t+1|t)$ are respectively expressed as follows [27]:

$$\widehat{\mathbf{x}}(t|t)={\left[\begin{array}{cc}\begin{array}{cc}\widehat{x}(t|t)& \widehat{x}(t-1|t-1)\end{array}& \widehat{x}(t-2|t-2)\end{array}\right]}^{\mathrm{T}}$$

(12)

$$\widehat{\mathbf{x}}(t+1|t)={\left[\begin{array}{cc}\begin{array}{cc}\widehat{x}(t+1|t)& \widehat{x}(t|t)\end{array}& \widehat{x}(t-1|t-1)\end{array}\right]}^{\mathrm{T}}$$

(13)

where $\widehat{x}(t+1|t)$, $\widehat{x}(t|t)$, $\widehat{x}(t-1|t-1)$, and $\widehat{x}(t-2|t-2)$ are all complex values, including magnitude and phase angle.

As for Φ(t + 1, t) and H(t + 1), they can be expressed as constant matrices shown in (14) and (15):

$$\mathsf{\Phi }\left(t+1,t\right)=\left[\begin{array}{ccccccccc}3& -3& 1;& 1& 0& 0;& 0& 1& 0\end{array}\right]$$

(14)

$$\mathbf{H}\left(t+1\right)=\left[\begin{array}{ccc}1& 0& 0\end{array}\right]$$

(15)

K(t + 1) is specified by the following formulas:

$$\{\begin{array}{l}\mathbf{K}\left(t+1\right)=\mathbf{P}(t+1|t){\mathbf{H}}^{\mathrm{T}}\left(t+1\right)\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\times {[\mathbf{H}\left(t+1\right)\mathbf{P}(t+1|t){\mathbf{H}}^{\mathrm{T}}\left(t+1\right)+\mathbf{R}\left(t+1\right)]}^{-1}\\ \mathbf{P}(t+1|t)=\mathsf{\Phi }\left(t+1,t\right)\mathbf{P}(t|t){\mathsf{\Phi }}^{\mathrm{T}}\left(t+1,t\right)\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}+\mathsf{\Gamma }\left(t+1,t\right)\mathbf{Q}\left(t\right){\mathsf{\Gamma }}^{\mathrm{T}}\left(t+1,t\right)\\ \mathbf{P}(t+1|t+1)=\left[\mathbf{I}-\mathbf{K}\left(t+1\right)\mathbf{H}\left(t+1\right)\right]\mathbf{P}(t+1|t)\end{array}$$

(16)

where R and Q are constant scalar values in this method. Furthermore, it is assumed that the Kalman filtering process starts at t = 3, $\widehat{\mathbf{x}}(3|3)$, and P(3|3) are the matrixes with 1 s and 0 s, respectively. This is to make the Kalman filter track the optimal estimate fast within a short period.

As the measurement noise covariance matrix R has a significant influence on the performance of Kalman filtering, this parameter will be adjusted dynamically. First of all, a residue between the a posteriori estimate and the actual measurement is calculated at each time t as follows [29]:

$$r\left(t\right)=\widehat{x}(t|t)-z\left(t\right)$$

(17)

where r(t) is the residue at time t, and z(t) = z(t). For a period, a residue vector is formed as:

$$\mathbf{r}\left(t+1\right)=\left[\begin{array}{ccc}\begin{array}{cc}r\left(t-W+1\right)& r\left(t-W+2\right)\end{array}& \begin{array}{cc}\cdots & r\left(t-1\right)\end{array}& r\left(t\right)\end{array}\right]$$

(18)

where r(t + 1) is the residue vector at time t + 1, and W is the time length. Then, the variance of the residue vector is calculated as the measurement noise covariance matrix shown below:

$$\mathbf{R}\left(t+1\right)=\mathrm{Var}\left(\mathbf{r}\left(t+1\right)\right)$$

(19)

where R(t + 1) is the measurement noise covariance matrix at time t + 1. Here, it is noted that Var(r(t + 1)) is the variance of a series of complex values. According to (18) and (19), it can be easily seen that R(t + 1) is obtained on the basis of W previous residues—namely, from time t − W + 1 to time t.

However, when the noise is weak, R(t + 1) approaches 0. As a result, when $\widehat{x}(t+1|t+1)$ is calculated, z(t + 1) is almost trusted fully, while $\widehat{x}(t+1|t)$ is hardly trusted, thus failing to reduce noise effectively. To solve this problem, R_min is introduced as the lower bound for R(t + 1):

$$R_{\min =\alpha }$$

(20)

where α is a small value. When R(t + 1) is less than R_min, R(t + 1) is set to R_min by default. Additionally, it is supposed that R(t + 1) = R_min within the initial period W.

Similarly, when missing data exist, which is processed as 0 in this paper, R_max is defined as the upper bound for R(t + 1) to enhance the trust in $\widehat{x}(t+1|t)$:

$$R_{\max =\beta }$$

(21)

where β is a large value. Specifically, in the presence of missing data at time t + 1, R(t + 1) is set to R_max to recover the missing data. This practice is based on the fact that z(t + 1) is far from the real value at time t + 1 and can be hardly trusted.

3. Real-Time Event Detection

To perform real-time analysis, an N × T real-time sliding window is utilized to obtain the raw matrix X from the raw PMU data source Ω. Let N denote the number of voltage magnitudes and phase angles and T denote the time period. At time t, the raw matrix X(t) is formed as follows:

$$\mathbf{X}\left(t\right)=\left[\begin{array}{ccc}\mathbf{x}\left(t-T+1\right)& \mathbf{x}\left(t-T+2\right)& \begin{array}{cc}\begin{array}{cc}\cdots & \mathbf{x}\left(t-1\right)\end{array}& \mathbf{x}\left(t\right)\end{array}\end{array}\right]$$

(22)

where $\mathbf{x}\left(t\right)={\left[\begin{array}{ccc}{x}_1\left(t\right)& x_2\left(t\right)& \begin{array}{cc}\cdots & x_N\left(t\right)\end{array}\end{array}\right]}^{\mathrm{T}}$ is the measurement data at time t.

In order to add white Gaussian noise into the original data source O to obtain the raw data source Ω, the signal-to-noise ratio (SNR) is defined as

$$\mathrm{SNR}=10\lg (\frac{\mathrm{Tr}\left(O{O}^{\mathrm{T}}\right)}{\mathrm{Tr}\left(G{G}^{\mathrm{T}}\right)\times m^2})$$

(23)

where Tr(·) denotes the matrix trace, G is the white Gaussian noise matrix, whose entries follow the standard normal distribution, m is the magnitude of white Gaussian noise, and G^T and O^T are the transposes of G and O, respectively.

During a period T_max, the raw data source Ω is obtained through the following formulas:

$$\mathsf{\Omega }=\mathbf{O}+m\mathbf{G}$$

(24)

$$m=\sqrt{\frac{\mathrm{Tr}\left(O{O}^{\mathrm{T}}\right)}{\mathrm{Tr}\left(G{G}^{\mathrm{T}}\right)\times {10}^{\mathrm{SNR}/10}}}$$

(25)

where Ω, O, and G are all N × T_max matrixes.

As the event indicator, MSR is calculated at each time interval and finally visualized to detect events. For simplicity, we set L = 1 when forming the matrix product Z.

To mitigate against false positives, the dynamic Kalman filter also serves as the data conditioner. Thus, after data conditioning, the conditioned raw matrix $\widehat{\mathbf{X}}\left(t\right)$ is formed as follows:

$$\widehat{\mathbf{X}}\left(t\right)=\left[\begin{array}{ccc}\widehat{\mathbf{x}}\left(t-T+1\right)& \widehat{\mathbf{x}}\left(t-T+2\right)& \begin{array}{cc}\begin{array}{cc}\cdots & \widehat{\mathbf{x}}\left(t-1\right)\end{array}& \widehat{\mathbf{x}}\left(t\right)\end{array}\end{array}\right]$$

(26)

where $\widehat{\mathbf{x}}\left(t\right)={\left[\begin{array}{ccc}{\widehat{x}}_1(t|t)& {\widehat{x}}_2(t|t)& \begin{array}{cc}\cdots & {\widehat{x}}_N(t|t)\end{array}\end{array}\right]}^{\mathrm{T}}$ is the conditioned data at time t. As a result, the procedure of real-time event detection based on random matrix theory and Kalman filtering is shown in

Table 1. Procedure of real-time event detection based on random matrix theory and Kalman filtering [30]. (Reproduced with permission from Fujia Han, IEEE Proceedings; published by IEEE 2018.).

Kalman Filter	1: Initialization
	2: Calculate ${\widehat{x}}_i(t|t)$ using (10), $1\le i\le N$
	3: Update ${\mathbf{R}}_i\left(t\right)$ using (19)–(21)
	4: Update ${\mathbf{K}}_i\left(t\right)$ using (16)
Event Detector	5: Form $\widehat{\mathbf{X}}\left(t\right)$ using (26)
	6: Transform $\widehat{\mathbf{X}}\left(t\right)$ to ${\mathbf{X}}_{\mathrm{s}}\left(t\right)$ using (6)
	7: Transform ${\mathbf{X}}_{\mathrm{s}}\left(t\right)$ to ${\mathbf{X}}_{\mathrm{u}}\left(t\right)$ using (7)
	8: Form $\mathbf{z}\left(t\right)$ using (2)
	9: Transform $\mathbf{z}\left(t\right)$ to ${\mathbf{z}}_{\mathrm{s}}\left(t\right)$ using (8)
	10: Calculate MSR using (5)
	11: If $t<T_{\max }$, repeat 2–10 at the next time t + 1; otherwise, go to 12
	12: Visualize MSR

.

4. Experimental Results and Discussion

In this section, the improved event detection approach is tested and compared with the work [19] using both simulated and real PMU data. The simulated data were generated based on the standard IEEE 118-bus system using MATPOWER [31], while the real PMU data were collected from the GB transmission system. Besides, voltage magnitude (in p.u.) and phase angle (in rad) were used, respectively, to validate the feasibility and effectiveness of the improved event detection approach. The SNRs of voltage magnitude and phase angle were kept the same in the tests. All the experiments were conducted on a common desktop using the matlab programming language.

For all the cases, we set Γ = 0.3I (I is a 3 × 3 matrix of 1 s), Q = 1 × 10⁻², α = 1 × 10⁻², β = 1 × 10⁶, W = 50, N = 118, and T = 240.

4.1. Data Conditioning

4.1.1. Noise Reduction

A comparison was performed between the dynamic and original Kalman filters. During the whole period T_max = 1500, there was a discrete event (an increase of active power demand at bus 60) happening at t = 601. In this case, the constant R of the original Kalman filter was 1 × 10⁻². The results are shown in

Table 2. Performance comparison of noise reduction between the dynamic and original Kalman filters at bus 60 (T_max = 1500).

SNR	RMSE of Voltage Magnitude			SNR	RMSE of Voltage Phase Angle
	Raw Values	Original KF	Dynamic KF		Raw Values	Original KF	Dynamic KF
5	17.519	8.6245	3.7744	5	6.6478	4.7623	3.8785
15	5.6129	2.7867	2.3817	10	3.5212	2.0497	1.5254
25	1.7739	0.8944	0.8906	15	2.0960	1.0613	0.9556
35	0.5534	0.2807	0.2765	20	1.1541	0.5932	0.5699
45	0.1738	0.0886	0.0877	25	0.6867	0.3389	0.3274
55	0.0550	0.0278	0.0276	30	0.3650	0.1867	0.1762
65	0.0174	0.0090	0.0088	35	0.2182	0.1118	0.1020
75	0.0055	0.0029	0.0027	40	0.1178	0.0663	0.0655
85	0.0018	0.0010	0.0009	45	0.0670	0.0424	0.0409
95	0.0006	0.0006	0.0005	50	0.0357	0.0326	0.0295

.

In order to eliminate the effect of the error, root mean squared error (RMSE) is calculated during t = 501–1500. In Table 2, it is clearly shown that the RMSE value of voltage magnitude of the dynamic Kalman filter is almost the same as that of the original Kalman filter when SNR is greater than or equal to 25, while the RMSE value of voltage magnitude of the dynamic Kalman filter is less than that of the original Kalman filter when SNR is less than or equal to 15. The RMSE value of voltage phase angle follows similar to that of voltage magnitude.

4.1.2. Missing Data Recovery

A further comparison was performed between the dynamic and original Kalman filters in terms of missing data recovery, as shown in

Table 3. Performance comparison of missing data recovery between the dynamic and original Kalman filters at bus 60 (T_max = 1500).

Percentage of Missing Data	RMSE of Voltage Magnitude (SNR = 95)			RMSE of Voltage Phase Angle (SNR = 35)
	Raw Values	Original KF	Dynamic KF (×10−4)	Raw Values	Original KF (×10−2)	Dynamic KF (×10−2)
5%	6.96	3.08	6	2.63	18.22	16.82
10%	10.02	4.62	7	3.71	18.27	16.84
15%	12.08	6.17	8	4.62	18.46	17.40
20%	14.12	7.50	9	5.33	18.63	17.87
25%	15.63	8.84	9	5.95	18.78	18.29
30%	17.17	10.37	10	6.52	18.90	18.79
35%	18.46	11.74	11	7.05	19.52	19.11
40%	19.85	12.96	12	7.53	19.88	19.68
45%	21.18	14.41	13	8.01	20.30	20.01
50%	22.11	16.03	13	8.42	20.82	20.71

. In this case, the same discrete event happened to bus 60, and RMSE was still calculated during t = 501–1500, and the constant R of the original Kalman filter remained 1 × 10⁻².

In Table 3, for voltage magnitude, it can be seen that the RMSE of the dynamic Kalman filter is significantly lower than that of the original Kalman filter. It is worth noting that as the percentage of missing data goes up, the RMSE of the dynamic Kalman filter only increases slightly from 6 × 10⁻⁴ to 13 × 10⁻⁴, while the RMSE of the original Kalman filter increases dramatically from 3.08 to 16.03. However, for voltage phase angle, it is shown clearly in Table 3 that the RMSE of the dynamic Kalman filter is only slightly lower than that of the original Kalman filter. Moreover, when the percentage of missing data rises, both cases go up rather slowly.

4.2. Event Detection Using Voltage Magnitudes

Event detection was conducted in situations with heavy noise and missing data using voltage magnitudes. Different signals were used to test the performance of the improved and original event detection approaches. The defined signals are shown in

Table 4. Defined signals for event detection.

Signal	Bus Number	Sampling Time	Active Power Demand (MW)
swell and sag signals	60	t = 1–600	fluctuation around 80
		t = 601–700	fluctuation around 120
		t = 701–1050	fluctuation around 80
		t = 1051–1150	fluctuation around 40
		t = 1151–1500	fluctuation around 80
	others	t = 1–1500	no change
multiple signals	60	t = 1–600	fluctuation around 120
		t = 601–700	fluctuation around 150
		t = 701–800	fluctuation around 180
		t = 801–900	fluctuation around 210
		t = 901–1000	fluctuation around 240
		t = 1001–1100	fluctuation around 270
		t = 1101–1200	fluctuation around 300
		t = 1201–1300	fluctuation around 330
		t = 1301–2000	fluctuation around 360
	99	t = 1–600	fluctuation around 60
		t= 601–700	fluctuation around 100
		t = 701–1200	fluctuation around 60
		t = 1201–1300	fluctuation around 20
		t = 1301–2000	fluctuation around 60
	others	t = 1–2000	no change

.

4.2.1. Heavy Noise

Figure 1 shows the MSR-t curves of the improved and original event detectors. It is noted that event detection starts after a short period because the sliding window needs to be filled up with samples and the dynamic Kalman filter takes a short time to track the real values.

In Figure 1a, it can be observed from the MSR-t curve of the improved event detector that signals can be detected easily. Specifically, during t = 260–600, MSR remains stable around 0.77, which means that the system state is normal without any signals. Then, at t = 600, MSR starts to decrease gradually until t = 770 (from 0.7733 to 0.7462). After that, MSR begins to increase to the original stable level (0.7674) until t = 940, then MSR remains almost constant. During this period, a U-shaped curve is observed, which indicates that there are signals changing the system state during t = 601–940. In addition, the effect of a signal on MSR extends T extra instants, due to historical data in the sliding window. As a result, the actual duration of the signals can be calculated as 940 – 601 + 1 – T = 100 instants, and the occurrence time of the signals is the starting time of the U-shaped curve, which is t = 601. That corresponds exactly to the fact that there is an increase of active power demand at bus 60 during t = 601–700. As a result, event detection is conducted in this way. Likewise, there is a second U-shaped curve from t = 1051 to t = 1390, which indicates that the signals occur at t = 1051, and the duration of the signals is 100 instants. However, for the original event detector, signals cannot be detected under the same situation. In other words, there is no clear U-shaped curve during the whole period because MSR changes slightly around 0.86 during t = 600–1390.

Apart from single signals, multiple signals can be also discovered by the improved event detector in a situation with heavy noise, as shown in Figure 1b. According to the MSR-t curve of the improved event detector, two types of signals can be detected. Firstly, two U-shaped curves are found during t = 601–940 and t = 1201–1540, and it indicates continuous signals during t = 601–700 and t = 1201–1300, respectively. Secondly, a third U-shaped curve is found during t = 601–1540, which indicates continuous signals during t = 601–1300. However, the same conclusion cannot be drawn from the MSR-t curve of the original event detector.

4.2.2. Missing Data

Figure 2 depicts the MSR-t curves of the improved and original event detectors when the percentage of missing data is 30%. As the dynamic Kalman filter also needs to take a short time to track the real values in the situation with missing data, event detection starts after a short period.

In Figure 2a, according to the MSR-t curve of the improved event detector, signals can be detected in the same way mentioned above. To be specific, two U-shaped curves are observed during t = 601–940 and during t = 1051–1390. However, the original event detector is unable to detect signals under the same situation because U-shaped curves cannot be found during the signal occurrence period. In addition, as shown in Figure 2b, multiple signals can still be discovered by the improved event detector in the situation with missing data. By contrast, the same signals cannot be discovered by the original event detector.

4.3. Event Detection Using Voltage Phase Angles

Event detection was also conducted using voltage phase angles with respect to the slack bus in situations with heavy noise and missing data. However, as voltage magnitude-based event detection is quite similar to voltage phase angle-based event detection, only swell and sag signals (as shown in Table 4) were employed to test both approaches.

Figure 3a describes the MSR-t curves of the improved and original event detectors when SNR is 20, while Figure 3b shows the MSR-t curves of the improved and original event detectors when SNR is 35 and the percentage of missing data is 30%. In Figure 3a,b, regardless of heavy noise or missing data, a U-shaped curve can be observed from t = 601 to t = 940 and from t = 1051 to t = 1390, respectively, for the improved event detector. By contrast, for the original event detector, there is no U-shaped curve found during the whole period in both situations.

4.4. Event Detection Using Real PMU Data

Real PMU data were also utilized to further validate the robustness of the proposed event detection approach. The PMU data were collected from 13 PMUs of the GB transmission system [32,33], and the sampling rate of the PMUs was 50 Hz. Besides, the total time length of the PMU data was 1 min, and thus there were 50 × 60 = 3000 samples in total for each PMU. What is more, the missing samples in the PMU data accounted for 12.5%. In this case, N was 13 and T was 500. The results are shown in Figure 4 and Figure 5.

In Figure 4a or Figure 5a, it can be observed that there is a U-shaped curve during t = 1777–2387 in the MSR-t curves of the improved event detector. This indicates that there are signals occurring at t = 1778 and lasting for 2387 − 1778 + 1 − T = 110, which is 110/50 = 2.2 s. According to Figure 4b or Figure 5b, the voltage magnitude and phase angle of a certain PMU change dramatically during t = 1778–1887, which confirms the above inference. In fact, there is indeed a tripping event from t = 1778 to t = 1887 across the GB transmission system. However, the original event detector fails to detect this event since no U-shaped curve is seen throughout the period.

5. Conclusions

This paper presents a robust event detection approach based on random matrix theory and Kalman filtering. In addition, as the data conditioner, a dynamic Kalman filtering technique was developed through the adjustment of the measurement noise covariance matrix in order to reduce noise and recover missing samples in PMU data. The experimental results show that the dynamic Kalman filter outperforms the original one, in terms of both noise reduction and missing data recovery. Furthermore, the improved event detection approach is much more robust than the original one, especially in practical situations where PMU data exhibit significant noise levels and large numbers of missing samples. The proposed event detection approach is also capable of performing computationally efficient event detection in real time in practice, without pre-training a model.

In this paper, the proposed approach only aimed to identify the occurrence time and the duration period of an event, without locating and classifying the event. Hence, in future research, the authors will focus on the location and classification of events to further develop and demonstrate the proposed robust approach in terms of event detection.