Online Measurement-Based Estimation of Dynamic System State Matrix in Ambient Conditions

TL;DR

This paper introduces a measurement-based, real-time method to estimate the dynamic system state matrix in power systems, independent of network and generator models, effective even under noise and missing data.

Contribution

It presents a novel recursive algorithm leveraging the Ornstein-Uhlenbeck process for accurate, model-independent estimation of the system state matrix in ambient conditions.

Findings

Accurately estimates the system state matrix during topology changes.

Robust against measurement noise and missing data.

Effective with higher-order generator models and control devices.

Abstract

In this paper, a purely measurement-based method is proposed to estimate the dynamic system state matrix by applying the regression theorem of the multivariate Ornstein-Uhlenbeck process. The proposed method employs a recursive algorithm to minimize the required computational effort, making it applicable to the real-time environment. One main advantage of the proposed method is model independence, i.e., it is independent of the network model and the dynamic model of generators. Among various applications of the estimated matrix, detecting and locating unexpected network topology change is illustrated in details. Simulation studies have shown that the proposed measurement-based method can provide an accurate and efficient estimation of the dynamic system state matrix under the occurrence of unexpected topology change. Besides, various implementation conditions are tested to show that the proposed method can provide accurate approximation despite measurement noise, missing PMUs, and the implementation of higher-order generator models with control devices.

Tables

Table 1. TABLE I : Notations for the power system dynamic model

${\delta }_i$	generator rotor angle
${\omega }_i$	generator angular frequency
$M_i$	inertial constant
$P_{mi}$	mechanical power input
$P_{ei}$	real power injection
$D_i$	damping coefficient
$n$	total number of buses
$E_i$	the voltage magnitude behind the transient reactance
$G_{ij}$	equivalent conductance between bus $i$ and $j$
$B_{ij}$	equivalent susceptance between bus $i$ and $j$

Equations

$$\dot{\mathbf{x}}$$

$$0$$

$$\dot{\mathbf{\delta }}$$

$$M\dot{\mathbf{\omega }}$$

$$P_{e_i}=\sum _{j=1}^n{E}_i{E}_j\left[G_{ij}\cos ({\delta }_i-{\delta }_j)+B_{ij}\sin ({\delta }_i-{\delta }_j)\right]$$

$$\dot{\mathbf{\delta }}$$

$$M\dot{\mathbf{\omega }}$$

$$\displaystyle\left[\begin{array}[]{c}{\Delta}\dot{\bm{\delta}}\\ {\Delta}\dot{\bm{\omega}}\end{array}\right]$$

$$\displaystyle+\left[\begin{array}[]{c}0\\ -M^{-1}E^{2}G\Sigma\end{array}\right]\bm{\xi}$$

$$\begin{gathered}\bm{x}=\left[\begin{array}[]{c}{{\Delta}\bm{\delta}}\\ {{\Delta}\bm{\omega}}\end{array}\right],\enskip A=\left[\begin{array}[]{cc}{{0}}&{I_{n}}\\ -M^{-1}\frac{\partial{\bm{P_{e}}}}{\partial{\bm{\delta}}}&-M^{-1}D\end{array}\right],\\ B=\left[\begin{array}[]{c}{0}\\ {-M^{-1}E^{2}G\Sigma}\end{array}\right]\end{gathered}$$

$$\dot{\mathbf{x}}=A\mathbf{x}+B\mathbf{\xi }$$

$$C=\left\langle\left[\bm{x}(t)-\mu_{\bm{x}}\right]\left[\bm{x}(t)-\mu_{\bm{x}}\right]^{T}\right\rangle=\left[\begin{array}[]{cc}C_{\bm{\delta}{\bm{\delta}}}&C_{\bm{\delta}{\bm{\omega}}}\\ C_{\bm{\omega}{\bm{\delta}}}&C_{\bm{\omega}{\bm{\omega}}}\end{array}\right]$$

$$G(\tau )=⟨\left[\mathbf{x}(t+\tau )-{\mu }_{\mathbf{x}}\right]{\left[\mathbf{x}(t)-{\mu }_{\mathbf{x}}\right]}^T⟩$$

$$\frac{d}{d\tau }\left[G(\tau )\right]=-AG(\tau )$$

$$A=\frac{1}{\tau }\log \left[G(\tau )C^{-1}\right]$$

$$\bar{\mathbf{x}}=\frac{1}{N}\sum _{i=1}^N{\mathbf{x}}_i$$

$$\hat{C}=\frac{1}{N}\left(F-\bar{\mathbf{x}}{1}_N^T\right){\left(F-\bar{\mathbf{x}}{1}_N^T\right)}^T$$

$$\hat{G}(△t)=\frac{1}{N}(F_{2:N}-\bar{\mathbf{x}}{1}_{N-1}^T)(F_{1:N-1}-\bar{\mathbf{x}}{1}_{N-1}^T{)}^T$$

$$A=\frac{1}{△t}\log \left[\hat{G}(△t){\hat{C}}^{-1}\right]$$

$${\bar{\mathbf{x}}}_N$$

$$=\left(1-\frac{1}{N}\right){\bar{\mathbf{x}}}_{N-1}+\frac{1}{N}{\mathbf{x}}_N$$

$${\hat{C}}_N^{-1}$$

$$=\frac{1}{(1-\alpha )}\left[{\hat{C}}_{N-1}^{-1}-\alpha \frac{{\hat{C}}_{N-1}^{-1}{\mathbf{z}}_N{\mathbf{z}}_N^T{\hat{C}}_{N-1}^{-1}}{1+\alpha {\mathbf{z}}_N^T{\hat{C}}_{N-1}^{-1}{\mathbf{z}}_N}\right]$$

$$A_0=\frac{1}{△t}\log \left[{\hat{G}}_0(△t){\hat{C}}_0^{-1}\right]$$

$$\begin{array}{c}{\bar{\mathbf{x}}}_j=(1-\alpha ){\bar{\mathbf{x}}}_{j-1}+\alpha {\mathbf{x}}_j\\ {\hat{G}}_j(△t)=(1-\alpha )\left[{\hat{G}}_{j-1}+\alpha ({\mathbf{x}}_j-{\bar{\mathbf{x}}}_j)({\mathbf{x}}_{j-1}-{\bar{\mathbf{x}}}_{j-1}{)}^T\right]\\ {\hat{C}}_j^{-1}=\frac{1}{(1-\alpha )}\left[{\hat{C}}_{j-1}^{-1}-\alpha \frac{{\hat{C}}_{j-1}^{-1}{\mathbf{z}}_j{\mathbf{z}}_j^T{\hat{C}}_{j-1}^{-1}}{1+\alpha ({\mathbf{z}}_j^T{\hat{C}}_{j-1}^{-1}{\mathbf{z}}_j-1)}\right]\end{array}$$

$$A_j=\frac{1}{△t}\log \left[{\hat{G}}_j(△t){\hat{C}}_j^{-1}\right]$$

$${\dot{\tilde{\delta }}}_i$$

$$M_i{\dot{\tilde{\omega }}}_i$$

$$+\frac{M_i}{M_T}\sum _{j=1}^3{E}_j^2{G}_{jj}{\sigma }_j{\xi }_j$$

$$\begin{array}{c}{M}_T=\sum _{j=1}^3{M}_j\\ {\delta }_0=\frac{1}{M_T}\sum _{j=1}^3{M}_j{\delta }_j,{\omega }_0=\frac{1}{M_T}\sum _{j=1}^3{M}_j{\omega }_j\\ {\tilde{\delta }}_i={\delta }_i-{\delta }_0,{\tilde{\omega }}_i={\omega }_i-{\omega }_0,i=1,2,3\end{array}$$

$$\begin{array}{c}{P}_{e_i}=\sum _{j=1}^3{E}_i{E}_j(G_{ij}\cos ({\tilde{\delta }}_i-{\tilde{\delta }}_j)+B_{ij}\sin ({\tilde{\delta }}_i-{\tilde{\delta }}_j))\\ P_{coi}=\sum _{i=1}^3(P_{m_i}-P_{e_i}),i=1,2,3\end{array}$$

$$A=\left[\begin{array}[]{cc|cc}0&0&1&0\\ 0&0&0&1\\ \hline\cr\lx@intercol\hfil\hbox{\multirowsetup$-M^{-1}\left(\frac{\partial\bm{P_{e}}}{\partial\bm{\tilde{\delta}}}\right)_{coi}$}\hfil\lx@intercol\vline&-\frac{D_{1}}{M_{1}}&0\\ &&0&-\frac{D_{2}}{M_{2}}\end{array}\right]$$

$${\left(\frac{\partial {\mathbf{P}}_{\mathbf{e}}}{\partial \tilde{\mathbf{\delta }}}\right)}_{coi}=\frac{\partial {\mathbf{P}}_{\mathbf{e}}}{\partial \tilde{\mathbf{\delta }}}+M\frac{1}{M_T}\frac{\partial P_{coi}}{\partial \tilde{\mathbf{\delta }}}$$

$${\left[{\left(\frac{\partial {\mathbf{P}}_{\mathbf{e}}}{\partial \tilde{\mathbf{\delta }}}\right)}_{coi}\right]}_{ii}$$

Full text

I Introduction

The primary goal of the power system operators is to maintain the system in the secure, reliable and efficient operating state despite varying operating conditions. Achieving this goal necessitates continuous monitoring of the system condition, conducting the periodic security analysis, and determining necessary control actions. The power flow Jacobian matrix and the dynamic system state matrix (also termed as the system dynamic state Jacobian in [1]) play essential roles in power system analysis since they contain invaluable information regarding the existence of the steady-state equilibrium point and the stability of the power system dynamic model, respectively. The critical difference between these two matrices is that the dynamic system state matrix entails the dynamic information such as the stability [1], the oscillation modes [2], the coherency groups [3] of the dynamical power system, which cannot be achieved from the power flow Jacobian matrix for general cases [1]. Note that these matrices can be easily obtained if the power system model and topology are known and accurate. However, the parameters of components may be inaccurate due to the errors introduced in the modeling stage (e.g., inaccurate manufacturing data and line length) or maintenance stage. Moreover, the topology may be corrupted due to erroneous telemetry from a remotely monitored circuit breaker. The above factors may lead to inaccurate or even incorrect results of static and dynamic security analysis.

Broad adoption of phasor measurement units (PMUs) and wide-area measurement system (WAMS) technology [4] provides a great opportunity to bypass this issue and gives birth to a new class of methods termed as measurement-based methods for power system monitoring, analysis, and control. The unique feature of these methods is that they rely mainly on high-speed, time-synchronized PMU data and are partially or entirely model-free. Significant attention have been giving to the PMU-based applications such as robust state estimation (e.g., [5]), generator dynamic state estimation (e.g., [6]), fault detection and location (e.g., [7], [8]), inter-area oscillations (e.g., [9], [10]), dynamic equivalent (e.g., [11], [12]), model calibration (e.g.,[13]), etc.

The broad adoption of PMUs also offers an alternative way to estimate the matrices mentioned above solely from measurement data. In previous literature, authors of [14] utilized high-speed synchronized voltage and current phasor data to determine the power flow Jacobian matrix through linear total least-square estimation; authors of [15] exploited the sparsity of the power flow Jacobian matrix and transformed its estimation into a sparse-recovery problem via compressed sensing; authors of [16] estimated the largest singular value of the inverse Jacobian matrix utilizing PMU data, which serves as a voltage stability indicator. However, the matrix considered in these works is the power flow Jacobian matrix where no dynamic information of the system is embedded. In our previous work[17], [18], we have proposed a hybrid measurement and model-based method to estimate the dynamic system state matrix which carries significant information about the dynamics of the system. The proposed approach is not purely measurement-based as we assume the inertia constants $M$ are known, but no information of network topology and parameters is needed.

In this paper, we make one further step by proposing a purely measurement-based method for estimating the dynamic system state matrix. Compared to [17],[18], the proposed approach is entirely model-free, i.e., it does not require any information about the network model and the dynamic generator model. This assumption relief is significant as the inertia and damping constants are also subject to inaccuracy and unavailability. To remove this assumption, the regression theorem of multivariate Ornstein-Uhlenbeck process [19] is utilized to estimate the dynamic state Jacobian matrix from the stationary time correlation and covariance matrix of measurements. To enable online application, a recursive estimation algorithm is developed to determine the matrix in near real-time. In sum, the contributions of the paper are presented as below:

• •

A model-free, purely measurement-based method is proposed to estimate the dynamic system state matrix by applying the regression theorem of the multivariate Ornstein-Uhlenbeck process.

• •

A fast recursive estimation algorithm is developed to minimize the required computational effort, enabling the application in the near real-time environment.

• •

An accurate estimate of the dynamic system state matrix can be obtained despite measurement noise, missing PMUs and the implementation of higher-order generator models with control devices. The topology change can be identified and located to the closest generator bus in terms of electrical distance.

The rest of the paper is organized as follows. Section II introduces the stochastic power system dynamic model and elaborates the proposed recursive measurement-based method for dynamic system state matrix. Section III presents the simulation studies to illustrate the feasibility and accuracy of the proposed method under various implementation conditions and demonstrate the application of the estimated matrix in detecting and locating network topology change.

II The system model and the proposed methodology

The general power system dynamic model can be described as below:

[TABLE]

where (1) describes dynamics of generators, excitation systems, load dynamics, etc., and (2) describes the transmission system and the static behaviors of passive devices. $f$ and $g$ are continuous functions; vectors $\bm{x}\in\mathbb{R}^{n_{\bm{x}}}$ and $\bm{y}\in\mathbb{R}^{n_{\bm{y}}}$ are the corresponding state variables (e.g., generator rotor angles and speeds) and algebraic variables (e.g., terminal bus voltage magnitudes and angles) [20].

In this paper, we are interested in small signal stability when the power system operates around its steady state, i.e., in ambient conditions. Under this condition, the dynamics of generators or aggregated generators can be represented by the swing equations. Thus (1)-(2) take the following form [21]:

[TABLE]

with

[TABLE]

where the notations are as in Table I.

In this paper, we make the common assumption [22] that the active powers of loads are perturbed by independent Gaussian variation from their base loadings. As shown in [23],[24], the load variations can be described by the stochastic perturbations at the diagonal elements of the reduced admittance matrix by $Y(i,i)=G_{ii}+jB_{ii}=Y_{ii}(1+\sigma_{i}{\bm{\xi}})\angle{\phi_{ii}}$, where $\bm{\xi}=[\dot{W_{1}},...,\dot{W_{n}}]^{T}$ is a vector of independent standard Gaussian random variable, and $\sigma_{i}^{2}$ denotes the variance of stochastic variation. Similar to the approach in [23], we assume that the power factor (and subsequently $\phi_{ii}$) remains unchanged, i.e., the load variations only reveal in the magnitude of admittance and not in the angle such that the nonlinearity does not kick in. Therefore, the stochastic power system model considering stochastic load variations can be described by:

[TABLE]

where $P_{e_{i}}$ is given in (5), $E=\mbox{diag}([E_{1},...,E_{n}])$, $G=\mbox{diag}([G_{11},...,G_{nn}])$, $\Sigma=\mbox{diag}([\sigma_{1},...,\sigma_{n}])$. Note that in the above formulation, system loads are modeled as constant impedances, since we assume that the system dynamics are dominated by generator angles. More realistic load models such as ZIP loads may be incorporated by adding additional terms in (5)[25], which on the other hand, does not affect the overall formulation.

Linearizing (6)-(7) around the steady state, we obtain:

[TABLE]

where $I_{n}$ is an identity matrix of size $n$. If we set

[TABLE]

then (8) can be represented by the following set of stochastic differential equations which is indeed a vector Ornstein-Uhlenbeck process [19]:

[TABLE]

The stationary covariance matrix is defined as:

[TABLE]

where $\mu_{\bm{x}}$ denotes the mean; the $\tau$-lag time correlation matrix is defined as:

[TABLE]

The regression theorem of multivariate Ornstein-Uhlenbeck process states that the following relation holds [19]:

[TABLE]

if the dynamic system state matrix $A$ is stable. In other words, if the system is in ambient conditions such that $A$ is stable, the dynamic system state matrix can be solved by:

[TABLE]

Equation (14) indeed is the theoretical basis of the proposed measurement-based method, which links the statistical properties of the measurements tactfully up with the physical model knowledge. Note that estimating $A$ via (14) is purely measurement-based since no knowledge of dynamic model parameters (e.g., damping and inertia constants) and network model (e.g., topology and line parameters) is assumed. Particularly, the proposed method requires no information of inertia constants compared to the hybrid method in [17]. It is also worth mentioning that the proposed method is different from system identification methods in the sense that the estimated matrix $A$ corresponds to the true dynamic system state matrix of the physical model whereas the matrix $A$ acquired from typical system identification methods does not necessarily correspond to the true dynamic system state matrix.

In Section II-B, we will develop a fast recursive estimation method based on (14) to estimate the dynamic system state matrix in near real-time. Before that, we will discuss the estimation of the time correlation matrix $G(\tau)$ and the covariance matrix $C$ from limited data set.

II-A Estimation of $G(\tau)$ and $C$

Ideally, $C$ and $G$ are calculated from an infinite number of data which is impossible in practice. Hence, the sample mean $\bar{\bm{x}}$, sample covariance matrix $\hat{C}$ and sample correlation matrix $\hat{G}$ calculated from finite data set is typically used as follows:

[TABLE]

[TABLE]

[TABLE]

where $N$ is the sample size, $F=[\bm{x}_{1},\bm{x}_{2},\cdots,\bm{x}_{N}]$ is a $M_{s}\times N$ matrix assuming there are $M_{s}$ state variables, $F_{i:j}$ denotes the submatrix of $F$ from $i$ to $j$ columns, $\boldsymbol{1}_{N}$ is an $N$ by 1 vector of ones, and $\triangle t$ is the sampling time step. Finally, the system state matrix can be estimated from:

[TABLE]

II-B The Recursive Estimation Algorithm

In order to enable online estimation, the recursive method is preferred which does not require calculating the basic equations over and over again. Once the initial $\hat{C}$ and $\hat{G}$ are estimated from a pre-defined length of data, they can be updated recursively using new data. For instance, the sample mean satisfies:

[TABLE]

Similar recursion for $\hat{C}$ and $\hat{G}$ can also be obtained.

Typically, we may generalize the recursion formulas by replacing $\frac{1}{N}$ by $\alpha$. How to choose the smoothing factor $\alpha$ will be discussed in the subsequent section. Moreover, the inverse covariance matrix can be calculated without inverting directly using Sherman-Morrison formula [26]:

[TABLE]

Before elucidating the algorithm, we make the following assumptions. First, PMUs are installed at all power plants where generators are connected, which is a reasonable assumption according to the basic guidelines for placing PMUs[27]. Second, we can use the PMUs to estimate the values of rotor angle $\bm{\delta}$ and rotor speed $\bm{\omega}$ in ambient conditions as explained in [28],[29],[30]. Note that such task may be difficult in transient condition and thus require dynamic system state estimation techniques, yet is achievable in ambient conditions [28]. Based on these assumptions, the following recursive algorithm is able to estimate the dynamic system state matrix in near real-time.

Step 1. Given a pre-defined window length $N$ data, initialize $\bar{\bm{x}}_{0}$, $\hat{C}_{0}$ and $\hat{G}_{0}(\triangle t)$ for $j=0$ by using (15)-(17), compute the inverse matrix $\hat{C}_{0}^{-1}$ and calculate the initial estimation for the dynamic system state matrix $A_{0}$:

[TABLE]

Step 2. For $j=1,2,\cdots$, perform the following:

• •

a) Update the sample mean, the sample correlation matrix and the inverse of the sample covariance matrix:

[TABLE]

where $\bm{z}_{j}=(\bm{x}_{j}-\bar{\bm{x}}_{j-1})$.

• •

b) Compute the dynamic system state matrix:

[TABLE]

II-C Implementation Issues

Like other measurement-based methods, the performance of the proposed recursive algorithm may be affected by parameter values such as the initialization window length and the smoothing factor $\alpha$.

Suppose the system is around the steady-state operating point, the larger initialization window will give more accurate estimations. Given the study in [17], the estimation error of the algorithm in [17] does not decrease substantially as the window length increases beyond 200s for the IEEE 39-bus system. Hence, we use 200s initialization window for the same 39-bus system in this paper which shows reasonably good accuracy (See Fig. 5 for more details). In practice, the selection of initial data length may be decided by the off-line study of the system of interests and the PMU data quality.

Another critical factor is the smoothing factor $0\leq\alpha\leq 1$. With $\alpha=\frac{1}{N}$, we assume the system operates in steady state, and hence all observations have the same weight independent on the time of occurrence; otherwise, a fluctuation is accepted such that the new observations may have a higher weight than the old ones which is indeed exponential smoothing, and the observation frame is approximately $\frac{1}{\alpha}$. Values of $\alpha$ close to 1 have a less smoothing effect and give greater weight to recent observations, while values of $\alpha$ closer to 0 have a greater smoothing effect and are less responsive to recent changes. In the limiting case with $\alpha=1$, the output is just the current observation. Generally, there is no reliable way of choosing $\alpha$, and it is mostly decided using heuristic methods [31]. In practice, $\alpha$ may be chosen adaptively based on the operating condition of the system. For instance, if the system is in steady state, let $\alpha=\frac{1}{N}$; if a contingency happens such that the system jumps to another operation state, choose a larger $\alpha$ to forget old samples more quickly and gradually decrease it until it reaches $\frac{1}{N}$.

In this paper, the adaptive exponential smoothing shown in Algorithm 1 is applied and integrated into the proposed algorithm in Section II-B. Particularly, $\beta$ and the change rate $w$ are two tunable parameters that affect how fast the algorithm could provide accurate estimation under a sudden change of the system, and describe the trade-off between the tracking capability under a sudden change and the accuracy of the estimation in steady state.

III Case Studies

III-A Numerical Example I: Estimating the Dynamic System State Matrix

First, we use a small system—the standard Western System Coordinating Council (WSCC) 3-generator, 9-bus system model to validate the proposed method, i.e., the accuracy of (18). The stochastic dynamic model integrating the load variations in the center-of-inertia (COI) formulation can be represented as [17], [23]:

[TABLE]

where

[TABLE]

[TABLE]

The model-based system state matrix is as follows:

[TABLE]

where

[TABLE]

The diagonal element and off-diagonal element of $\left(\frac{\partial\bm{P_{e}}}{\partial\bm{\tilde{\delta}}}\right)_{coi}$ can be calculated by

[TABLE]

Note that the conventional model-based method of estimating $(\frac{\partial\bm{P_{e}}}{\partial\bm{\tilde{\delta}}})_{coi}$ and thus $A$ heavily depends on the network topology and parameter values embedded in the admittance matrix $Y$. In contrast, the proposed technique, i.e., (18), is purely measurement-based, completely independent of the network model and the generator parameters.

Assuming $\sigma_{1}=\sigma_{2}=0.01$ p.u. describing the standard deviations of the load variations, Fig. 1 presents the emulated $\tilde{\delta}_{1}$ and $\tilde{\omega}_{1}$ estimated from the PMU data with a sampling rate of 50 Hz, from which we see that the system’s trajectories are oscillating around their steady-state solutions (denoted by the red lines) because of the stochastic load variations.

If there is no stochastic load variation, i.e., $\sigma_{1}=\sigma_{2}=0$ p.u., the dynamic system state matrix $A$ would be a constant matrix calculated from (27)-(29) using the steady-state operating point:

[TABLE]

To validate the proposed algorithm, we will show that the dynamic system state matrix estimated by the proposed algorithm is close to the deterministic system state matrix (30). To this end, we first use the $200$s data to estimate the sample covariance matrix and the sample one-lag correlation matrix by (15)-(17):

[TABLE]

[TABLE]

and then estimate the dynamic system state matrix $A$ by (18):

[TABLE]

where ⋆ denotes the estimation by the proposed method. It is observed that the estimated dynamic system state matrix (31) is close to the deterministic matrix shown in (30). Particularly, we use the normalized Frobenius norm to quantify the difference between the actual and estimated matrices:

[TABLE]

where $\|\|_{F}$ denotes the Frobenius norm. Note that we only use 200s data for demonstration. However, the accuracy will improve as the sampling rate increases and/or the length of measurement data grows. There is always a trade-off between the computational effort and the accuracy.

In addition, we compare the accuracy of the proposed method with that of the method developed in [17]. The method in [17] are developed based on another important property of the stationary covariance matrix $C$ that $C$ satisfies the Lyapunov equation in ambient condition:

[TABLE]

After algebraic reduction, we have:

[TABLE]

Assuming that the contributions from the matrix $C_{\bm{\delta\omega}}$ are negligible, it yields $\frac{\partial{\bm{P_{e}}}}{\partial{\bm{\delta}}}\approx MC_{\bm{\omega\omega}}C^{-1}_{\bm{\delta}\bm{\delta}}$. It is obvious that this method is not purely measurement-based since the dynamic system state matrix $\frac{\partial{\bm{P_{e}}}}{\partial{\bm{\delta}}}$ depends on the inertia constants $M$. If $A$ needs to be estimated, $D$ has to be known. In this numerical example, the estimation error defined in (32) is $3.73\%$ using (34), which is less than that of the proposed method. Nevertheless, since the method in this paper and the one in [17] are developed based on different properties of the stochastic dynamical system (10), there is no direct and clear conclusion that which method should possess a higher accuracy, even though the method proposed in this paper requires less model information and is completely model-free.

III-B Numerical Example II: Detecting Topology Change and Identifying the Location

In this example, we apply the recursive estimation method to a larger system—the IEEE 39-bus 10-generator test system to demonstrate that the proposed methodology is able to detect the topology change and more importantly, identify the corresponding location. This information is of great significance for power system operators to adopt corrective control strategies to fix the error and maintain the security of power grids.

For this purpose, we assume that an undetectable line tripping between Bus 22 and Bus 23 occurs at 400s due to, for instance, communication errors. As shown in Fig. 2, the topology change location is close to Generator 6 and 7. The trajectories of the rotor angles of Generator 6 and 7 are shown in Fig. 3, from which we see that there is an abrupt change in $\tilde{\delta}_{6}$ at 400s. Nevertheless, the post-fault system still maintains the stability.

Suppose the change of $\bm{\delta}$ and $\bm{\omega}$ can be detected by the PMU measurements, yet the network topology change is not detected. Denote the dynamic system state matrix estimated by the model-based method without accounting for the topology change as $A^{\diamond}$, it is found that the normalized distance between $A^{\diamond}$ and the true dynamic system state matrix $A$ after the change is $17.90\%$. The big discrepancy is due to the incorrectly assumed network model. In contrast, the normalized distance between the dynamic system state matrix $A^{\star}$ estimated by the proposed method and the true one is equal to $2.72\%$. Particularly, the data on [410s 1200s] is used to estimate the covariance matrix and the one-lag correlation matrix. It is obvious that the proposed method provides more accurate estimation than the model-based method under the occurrence of undetectable topology change. In addition, the big difference between the model-based and the measurement-based matrices indicates that the assumed system model is inaccurate and identifying the location of the error is stringent. It will be shown next that the proposed recursive matrix estimation algorithm can also help identify the location of the topological error.

To implement the recursive algorithm, we set the initial window length to be 200s, i.e., $N=10000$ for a sampling rate of $50$ Hz; $\beta=200$ and $w=2$ in the adaptive algorithm for $\alpha$ (Algorithm 1). A small $\beta$ means that a large smooth factor $\alpha$ is applied to quickly forget the previous samples and $w$ describes a gradual decrease in $\alpha$ until it reaches its steady-state value $\frac{1}{N}$ or when the system returns to the steady state, after which $\alpha$ is reset to be $\frac{1}{N}$. The normalized Frobenius norm between the estimated dynamic system state matrix and the true one is shown in Fig. 4(a), from which we see that the proposed method is able to provide an accurate estimation for the dynamic system state matrix both before and after the line tripping. It should be noted that the normalized distance shown in Fig. 4 (a) by the recursive algorithm is higher than $2.72\%$ when the data on [410s 1200s] is used. Such difference lies in the fact that different window lengths are applied, i.e., 790s vs. 200s. The impact of window length will be further discussed in Fig. 5. Besides, the errors in the transient period 400s-600s are due to the fact that it takes some time for the recursively estimated covariance matrix, the one-lag correlation matrix, and thus the system state matrix to settle down to their new values corresponding to the new steady state after the topology change. In addition, the normalized Frobenius norm between the estimated dynamic system state matrix by the model-based method (i.e., without detecting the topology change) and the true one is shown in Fig. 4(b), showing that the model-based method matrix becomes inaccurate after the undetected topology change. Given the fact that the true dynamic system state matrix is typically unknown in practice, we also compare the normalized distance between the dynamic system state matrix estimated by the proposed recursive method and the one obtained from the conventional model-based method. As shown in Fig. 4(c), the two estimated matrices are close to each other with a difference around 3% (both the model-based method and the proposed measurement-based method work well), yet they deviate from each other after the line tripping due to the incorrectly assumed network model used by the model-based method. By monitoring the difference between the two estimated matrices, power system operators can quickly detect that an error occurs in the assumed network model at 400s and locating the error is imperative.

To justify that a window length of 200s is reasonable, we compare the accuracy of the estimation when the window length changes from 10s to 900s as shown in Fig. 5. It can be seen that the estimation accuracy does not improve dramatically when the window length goes beyond 200s. Considering the compromise between accuracy and efficiency, we choose $200$s to be the length of the initial window and the moving window in the recursive algorithm.

Regarding the efficiency, the proposed recursive algorithm is not computationally expensive since it involves only basic arithmetic operations as seen from (21)-(23). In our simulation study, the average time consumed for each time instant $j$ in Step 2 where the main computational effort lies in is $1.8592\times 10^{-4}$s for the 3-generator system (Numerical Example I in Section III-A), and $8.7911\times 10^{-4}$s for this 10-generator system. Note that there is always one generator being the reference. It can be observed that the time consumption grows approximately linearly with the number of generators. The efficiency can be further improved to conduct the recursive estimation less often in the steady state, whereas more often after the system’s evolution due to, for instance, topology change.

To better visualize the difference between the true matrix, the estimated matrix by the proposed method, and the estimated matrix by the model-based method, we use the 3D bar graph to show the difference between each component of the matrices. Denote $A^{\diamond}(t)$ as the matrix estimated by the model-based method using the measurement at $t$, and $A^{\star}(t)$ as the matrix estimated by the recursive measurement-based method at time $t$, then Fig. 6 presents the difference between $|a_{ij}-a_{ij}^{\diamond}(1200)|$ and $|a_{ij}-a_{ij}^{\star}(1200)|$, where $a_{ij}\in{A}$ — the true system matrix, $a_{ij}(1200)^{\diamond}\in{A^{\diamond}(1200)}$, and $a_{ij}(1200)^{\star}\in{A^{\star}(1200)}$. Furthermore, Fig. 7 directly shows the difference between the component of the two estimated matrices which we can acquire online, i.e., $|a_{ij}^{\diamond}(1200)$-$a_{ij}^{\star}(1200)|$. It is evident from the figures that the proposed measurement-based method provides more accurate estimation than the model-based method. More importantly, the proposed method provides important insights regarding the location of the error. Since the largest discrepancies between the two estimated matrices (Fig. 7) occur in elements $a_{15,6}$, $a_{15,7}$, $a_{16,6}$, and $a_{16,7}$ corresponding to the ${\frac{\partial P_{e}}{\partial\delta}}$ of Generator 6 and 7, it indicates that the topology change is electrically close to Generator 6 and 7, which aligns exactly with the fact.

The results of the above two cases entail that the proposed recursive measurement-based method is able to provide a good estimation for the dynamic system state matrix in near real-time by utilizing the inherent statistical properties of the stochastic system. Moreover, the method may help detect the topological errors in the assumed network model, and more importantly, identify their sources.

III-C Numerical Example III: Identifying the Location of a Topology Change with Missing PMUs

It is assumed that all generators are monitored by PMUs when developing the methodology, which seems to be too strict. In this section, we relieve this assumption and consider a more realistic situation where some PMUs are missing.

Assuming the PMUs are only available for Generators $K$, where $K=\{m_{1},...,m_{p}\}$, then we can still accurately estimate the sample mean $\bar{x}_{i}$ where $i\in K$, and the sub-matrices $\hat{C}_{ij}$ and $\hat{G}_{ij}$, where $i,j\in K$, of the sample covariance matrix $\hat{C}$ and the sample correlation matrix $\hat{G}$, respectively. As a result, by (14), we can accurately achieve a sub-matrix $A_{ij}$, where $i,j\in K$, of the dynamic system state matrix $A$. More importantly, the estimated sub-matrix can still facilitate the detection and location of the topology change.

To show this, we reconsider the IEEE 39-bus 10-generator system in Section III-B. Assuming the PMUs at Gen 3 and Gen 9 are missing, i.e., $K=\{1,2,4,5,6,7,8\}$ and $\delta_{10}$ is the reference, we keep estimating the difference between the two sub-matrices recursively by the proposed measurement-based method and by the model-based method, respectively. The Frobenius norm of the difference with respect to the time is presented in Fig. 8, from which it can be concluded that the error in the network model occurring at 400s was not well detected. To locate the error, we further draw the 3D bar graph to spot the difference between each component of the two sub-matrices at a certain time step after the fault, say 1200s, as shown in Fig. 9. Note that only the submatrix ${-M^{-1}(\frac{\partial\bm{P_{e}}}{\partial\bm{\tilde{\delta}}})_{coi}}$ in (27) is presented since it possesses the largest discrepancies as revealed in Fig. 6-7. Although the estimation for $A_{3j}$, $A_{i3}$, $A_{9j}$ and $A_{i9}$ are missing, we still can observe that the largest discrepancies occur in Gen 6 and 7, indicating that the error is electrically close to these two generators which well aligns with the fact.

It is worth mentioning that pinpointing the exact location of the topology error for a large system may require additional information since the proposed method relies on the reduced Y-matrix for the detection, in which only the internal generator buses are retained. Nevertheless, the proposed method may help narrow the range of topology change and identify the generators that are electrically close to the topology change without any knowledge of the network model.

III-D Impact of Measurement Noise

Similar to all the measurement-based methods, the noise may affect the accuracy of the proposed method. To investigate this issue, we added a measurement noise with a standard deviation $10^{-3}$ to both $\bm{\delta}$ and $\bm{\omega}$ of Numerical Example II according to the IEEE Standards [33]. Likewise, Fig. 10(a) shows the recursive estimation of the normalized distance between the estimated dynamic system state matrix by the proposed method and the true one; Fig. 10(b) presents the normalized distance between the matrix estimated by the proposed measurement-based method and the one by the model-based method. For comparison purpose, the results with no measurement noise are also plotted. It can be found that the impact of measurement noise on the accuracy of the method is minor since only a slight increase in the error is observed, which does not affect the detection of the topology error. 3D bar graph similar to Fig. 6-7 can also be obtained, which provide essential information to locate the undetected topology change.

III-E Validation on a High-Order Model

Despite the fact that the methodology is developed using the classical second-order generator model, it is necessary to test the performance of the methodology in a real-life setting in which higher-order models controlled by exciters are typical. To this end, we revisit the WSCC 3-generator, 9-bus system in Numerical Example I. Specifically, all generators are modeled as third-order models equipped by automatic voltage regulators. Assuming an undetected line tripping between Bus 7 and Bus 8 occurs at 300s, we apply the proposed recursive algorithm to estimate the dynamic system state matrix, detect the topology change and locate the error. Simulation has been carried out in PSAT-2.1.10 [34].

Particularly, the initial window length is set to be 200s; $\beta=10$ and $w=0.2$ are chosen for calculating the adaptive smooth factor $\alpha$. As mentioned earlier and revealed by Fig. 6-7, the submatrix ${-M^{-1}(\frac{\partial\bm{P_{e}}}{\partial\bm{\tilde{\delta}}})_{coi}}$ typically possesses the largest difference between matrices, which is therefore a good index to demonstrate the estimation accuracy, detect and locate topological change. The normalized difference between the true submatrix and the one estimated by the proposed method is shown in Fig. 11(a), from which it can be seen that the proposed method can still provide accurate approximation regardless the implementation of third-order generator models with exciters. In addition, detecting topology change and identifying its location can still be realized by monitoring the difference between the submatrix estimated by the proposed measurement-based method and the one obtained from the model-based method, as illustrated in Fig. 11(b).

Similar results can be obtained using the fourth-order model equipped by both automatic voltage regulators and turbine governors as shown in Fig. 12. Although a small degradation in accuracy and convergence speed is noticed, the proposed method can still identify the topology change accurately and quickly, indicating the good feasibility of the proposed method in practice.

III-F Example IV: IEEE 145-bus System

In order to validate the proposed method on large test systems, we applied it to the IEEE 145-bus dynamic stability test case. This system contains 50 generators and 453 branches. The power flow and dynamic parameters of this system are available online: http://labs.ece.uw.edu/pstca/dyn50/pg_tcadd50.htm.

The undetected topology change occurs between bus 117 and bus 118. The recursive algorithm is implemented with a window length of 200s. $\beta=200$, and $\omega=0.3$. Since the errors in the transient is too large, the normalized distance in the steady-state after the topology change have been zoomed in, showing that the proposed method provides more accurate estimation ($\approx 8\%$) for the dynamic system state matrix than the model-based method ($\approx 25\%$) under the undetectable topology change.

IV Conclusions and Perspectives

In this paper, we have proposed a novel recursive measurement-based method for estimating the dynamic system state matrix in near real-time. Particularly, the proposed method leverages the statistical properties of the stochastic dynamical system to extract significant model information out of PMU data. A recursive estimation algorithm has also been developed to enable online implementation. The simulation study has shown that the proposed recursive method can provide an accurate approximation for the dynamic system state matrix, and more importantly, detect severe discrepancies in the assumed network model and further identify their locations. Various implementation conditions including measurement noise, missing PMUs, and detailed generator models have been tested to demonstrate the feasibility and the accuracy of the proposed method in real-life applications.

In the future, we plan to explore diverse applications of the estimated dynamic system state matrix in online oscillation analysis, model validation, emergency power re-dispatch, and load shedding, etc. Besides, we intend to use actual, rather than emulated, PMU measurement to validate the proposed methodology in practice.