Location via proxy:   [ UP ]  
[Report a bug]   [Manage cookies]                
Next Article in Journal
Optimal Determination Method of the Transposition Steps of An Extra-High Voltage Power Transmission Line
Previous Article in Journal
Determinants of Managerial Competences Transformation in the Polish Energy Industry
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Hybrid Physics-Based Adaptive Kalman Filter State Estimation Framework

by
Arturo S. Bretas
1,*,
Newton G. Bretas
2,
Julio A. D. Massignan
2 and
João B. A. London Junior
2
1
Department of Electrical & Computer Engineering, University of Florida, Gainesville, FL 32611-6200, USA
2
Department of Electrical & Computer Engineering, University of Sao Paulo, Sao Carlos 13566-590, Brazil
*
Author to whom correspondence should be addressed.
Energies 2021, 14(20), 6787; https://doi.org/10.3390/en14206787
Submission received: 31 August 2021 / Revised: 7 October 2021 / Accepted: 13 October 2021 / Published: 18 October 2021

Abstract

:
State-of-the art physics-model based dynamic state estimation generally relies on the assumption that the system’s transition matrix is always correct, the one that relates the states in two different time instants, which might not hold always on real-life applications. Further, while making such assumptions, state-of-the-art dynamic state estimation models become unable to discriminate among different types of anomalies, as measurement gross errors and sudden load changes, and thus automatically leads the state estimator framework to inaccuracy. Towards the solution of this important challenge, in this work, a hybrid adaptive dynamic state estimator framework is presented. Based on the Kalman Filter formulation, measurement innovation analytical-based tests are presented and integrated into the state estimator framework. Gross measurement errors and sudden load changes are automatically detected, identified, and corrected, providing continuous updating of the state estimator. Towards such, the asymmetry index applied to the measurement innovation is introduced, as an anomaly discrimination method, which assesses the physics-model-based dynamic state estimation process in different piece-wise stationary levels. Comparative tests with the state-of-the-art are presented, considering the IEEE 14, IEEE 30, and IEEE 118 test systems. Easy-to-implement-model, without hard-to-design parameters, build-on the classical Kalman Filter solution, highlights potential aspects towards real-life applications.

1. Introduction

Electric power systems are critical infrastructures that rely on continuous and detailed monitoring of their assets, encompassed by power system state estimation (PSSE). Since its conception in the early 1970s, PSSE has been focused on processing the steady-state condition of the transmissions systems, under a static analysis framework [1]. With the penetration of renewable energy and new control equipment, the static perspective does not provide the modern energy management systems’ requirements. Thus, it becomes imperative to capture systems’ dynamics for advanced automated applications, in such a way to provide fast control actions, increase the reliability of the system, diagnosis of cyber-attacks, and optimize resources across the power grid [2,3,4].
Dynamic state estimation in power systems can be understood as to estimate the grid states, as voltages magnitude and phases, considering the time relation between them. There are many presented physics-model-based dynamic state estimation (DSE) solutions, being the Kalman Filter formulation and it is likewise a family of formulations most used [3,5,6]. The first formulations of DSE in power systems were proposed shortly after the consolidation of static state estimation in transmission system operation centers [7,8,9,10,11,12].
Such works established the theoretical foundations to capture the dynamic characteristic of the power system. They encompass a first-order state-space model to represent the temporal relationships among the state variables and the Kalman Filter technique application. However, the computational burden of the algorithm, and the lack of fast updating measurements in the Supervisory Control and Data Acquisition (SCADA) system hindered real-life applications [3,5,12]. Only in the late 2000s, dynamic state estimators regained significant interest, with the advent of faster computer architectures, by the allocation of new instrumentation and communication technologies with faster updating and sampling rates, and as well as the acute power system operational requirements, carried by a more volatile environment with renewables [5]. In the effort of surpassing many of the practical and theoretical challenges, different versions of the Kalman Filter were proposed, such as the Linear Kalman Filter [13], the Extended Kalman Filter [14], the Unscented Kalman Filter [15], the Cubature Kalman Filter [16], the Correntropy Kalman Filter [17] and the Ensemble Kalman Filter [18]. Nonetheless, all such approaches consider a single stationary scenario or that only known load changes occur on the system. This translates to the assumption that the system’s transition matrix is always correct, which might not hold in real-life applications. In a complementary perspective, robustness against measurement gross errors for dynamic state estimators gained a significant research effort. The work in [19] shows that incorporating the state-space model can benefit classical bad-data processing based on the weighted residual analysis. Also, recent works seek to improve the estimator robustness by changing its objective function, such as by the M-estimators [20], the generalized maximum likelihood estimation [21], the weighted least absolute value [22], the H-infinite [23] and hybrid approaches [24]. Altogether, robustness is often obtained at the cost of improving the estimator variance, computational burden, or model complexity. Conversely, such approaches fail to discriminate measurement gross errors from changes in the stationary load level, often treating them in the same theoretical perspective, or by modifying the residual-based tests. Amidst such a promising challenge, this work explores a family of Kalman Filters formulation for the DSE problem. The state space is modeled under a piece-wise stationary state space, and the change in the stationary load level must be detected and discriminated from measurement errors. Towards such a goal, a statistical asymmetry index, applied to the system innovation, is introduced. This metric is used in analytical-based tests towards continuous discrimination of anomalies type. In this work, it is further shown that considering the condition of changing the stationery load level, the asymmetry index will have a non-zero small value. In this case, the innovation will not have Normal distribution. Otherwise, in the measurement gross error condition, the innovation will present a large asymmetry index, due to the modeling structural discontinuity. Therefore, the presented analytical-based tests enable discrimination among different types of anomalies and automatically adapt the dynamic state estimator to maintain accuracy. The specific contributions of this work towards the state-of-the-art are:
  • The development of a hybrid state estimator framework, based on the Kalman Filter formulation, able to diagnosis measurement gross errors and sudden load changes automatically;
  • The exploration of the asymmetry index as an anomaly discrimination method, assessing the DSE process in different piecewise stationary levels, which can be used in any family of Kalman Filters-like dynamic state estimators framework.
The remaining of the manuscript is structured as follows: Section 2 presents the theoretical background for DSE and the Kalman Filter (KF). Section 3 introduces the dynamic estimator and the innovation analysis for anomaly discrimination. Section 4 presents test results, followed by the conclusions, which are presented in Section 5.

2. Dynamic State Estimation

DSE in power systems can be understood to estimate the grid states, as voltages magnitude and phases, considering the time relation between them. DSE in power systems can be modeled as
x ̲ k + 1 = F k x ̲ k + G k w ̲ k + Γ k u ̲ k z ̲ k = H k x ̲ k + v ̲ k ,
With x ̲ k the state vector at time k, w ̲ k a known input and u ̲ k the state noise. Equivalently, z ̲ k the observation vector and v ̲ k the observation, or measurement, noise. F k is the state transition matrix, G k is the matrix relating the states to the input w ̲ k and H k is the matrix relating measurements to states at time k, all of them with the appropriate dimensions. Figure 1 illustrates a block diagram for the signal model with external known input that corresponds to the DSE problem.
Within the power systems realm, dynamic state estimation has been conceptualized under different perspectives, mainly related to the definition of state variables. The fundamental difference comprises the inclusion, or not, of internal variables from generators and loads (also referred to as dynamic state variables), besides the complex nodal voltages of the power network (also referred to as algebraic state variables) [3,5]. The first focuses on fast transient events that may cause instabilities, while the second on slower rates of change that affect the loading condition of the power networks. Nonetheless, both represent a system subject to a stochastic phenomenon and uncertainty. The difference lies in the time scale of events and observed response of the system. In this work, dynamic state estimation refers to the second definition, also referred to as tracking or forecasting-aided state estimation [4], focused on the power flow behavior over time, that is, the temporal relations among different operational conditions of the interconnected power grid. Hence the state variables comprise the complex nodal voltages in all nodes of the power systems x = { ( V i , θ i ) | i { 1 , 2 , , n b u s } } . The measurements comprise electrical quantities monitored in real-time by the SCADA systems. A detailed discussion about the employed models and physical aspects of the power system is presented in the Appendix A.

2.1. Statement of the Discrete Kalman Filtering Problem

Considering the linear, finite-dimensional, discrete system described by the previous equations, with v ̲ k and w ̲ k being independent, Gaussian white noise processes with
E [ v ̲ k v ̲ j T ] = R k a n d E [ w ̲ k w ̲ j T ] = Q k ,
where R k is the measurement noise covariance matrix and Q k is the state covariance matrix. Suppose x ̲ 0 , the initial value of the system state variables, being a Gaussian random variable. Assuming they have mean x ̲ ¯ 0 and covariance matrix 0 | 1 , the KF solution consists of estimating x ̲ k | k 1 and x ̲ k | k and their associated error covariance matrices k | k 1 and k | k . Here x ̲ k | k 1 means estimating the state vector at time step k with data until time step ( k 1 ). If the assumption of a Gaussian distribution for x ̲ 0 , υ ̲ k , and w ̲ k is dropped, but they remain described by their first and second-order statistical moments of their probability distribution functions, all the conclusions about unbiasedness and algorithm formulation are precisely the same. However, optimality, in the sense of the minimum variance estimator and minimum squared error, will not be the case. This is one of the major motivations to overcome gross errors in the measurement vector, that largely distort such assumption of Gaussian noise and negatively affects the estimation.

2.2. Solution of the Kalman Filtering problem

The KF model can be solved in two steps: the prediction and the estimation step.
(1) SKF Prediction Step:
The Kalman Filter prediction step is given by
x ̲ ^ k + 1 | k = F k x ̲ ^ k | k 1 + K k [ z ̲ k H k x ̲ ^ k | k 1 ] + Γ k u ̲ k ,
With x ̲ ^ 0 | 1 = x ̲ ¯ 0 . The Kalman gain is given by:
K k = F k k | k 1 H k T [ H k k | k 1 H k T + R k ] 1 ,
the inverse is guaranteed with R k being a positive definite matrix, while the conditional error covariance matrix k | k 1 is computed recursively by the Riccati equation:
k + 1 | k = F k [ k | k 1 H k T ( H k k | k 1 H k T + R k ) 1 H k k | k 1 ] F k T + G k Q k G k T ,
where, k + 1 | k is the updated error covariance matrix.
(2) KF Estimation Step:
The KF estimating solution is given by
x ̲ ^ k | k = x ̲ ^ k | k 1 + k | k 1 H k T ( H k k | k 1 H k T + R k ) 1 ( z ̲ k H k x ̲ ^ k | k 1 ) ,
assuming that the inverse exists, the error covariance matrix is given by
k | k = k | k 1 k | k 1 H k T ( H k k | k 1 H k T + R k ) 1 H k k | k 1 .

2.3. Bad Data Analysis

Considering DSE, besides measurement errors, one can have other anomalies as large load changes. In this case, the state transition matrix may not reflect the hull transition of states. This highlights the difficulty in modeling the system for DSE purposes. The state uncertainty can be quantified by white Gaussian noise with covariance Q k . This technique may incur some difficulties: if one assumes a high covariance so that the sudden load change is inside the modeling, the estimate may be significantly poor for the case of small load fluctuations. Otherwise, if the noise covariance is considered small, to represent well the small load fluctuations, the model will not simulate large load fluctuations, but, in this case, it could be used to identify the anomaly for large load variations [10].
To tackle the previous challenge, in this work, as the state-of-the-art, the state transition matrix is modeled through an identity matrix, with a small covariance matrix. This will represent a stationary system behavior. As a tradeoff, possible changes in the stationarity load level will be detected. However, measurement gross errors could be misidentified as a stationarity level change. Sudden large load changes, as well as measurement gross errors, will be named here as an “anomaly”. That said, they are anomalies of different characteristics and should be, in this way, automatically detected and identified. One should note that in the bad data condition, there exist equations “inconsistencies” or “a structural discontinuity” of the equations describing the DSE behavior [10,25]. Otherwise, considering the large load changes, there is no structural discontinuity in system equations [10,17,25]. In this work, the distinction between these anomalies is done through the analysis of the statistical distribution of the innovation process.
Due to the structural discontinuity in the equations, the innovation for the measurement gross error condition will present a large asymmetry characteristic. Otherwise, considering stationarity load level changes, the symmetry related to the innovation average value will keep its value, but dislocated from zero. In this case, the innovation will not have a Normal distribution [10]. These characteristics are illustrated in Figure 2.

3. Dynamic State Estimation Formulation for Gross Error Analysis

The power system state can be defined by an n-dimensional x ̲ vector consisting of the voltages magnitudes and their phase angles, except one of them, taken as reference. The system state can be modeled as
x ̲ k + 1 = x ̲ k + w ̲ k .
With w ̲ k being a white Gaussian, noise vector with
E [ w ̲ k w ̲ j ] = Q k j a n d Q = α 2 d i a g ( q i 2 ) = α 2 Γ .
The previous q i is the maximum ratio of the state i obtained from historical data of bus i. To represent power flows, the measurement model in power systems dynamic state estimation consists of the following set of nonlinear equations ( h ( x k ) ) given by
z ̲ k = h ( x ̲ k ) + υ ̲ k .
and
E [ υ ̲ k υ ̲ j T ] = [ δ k j R k j ] , E [ υ ̲ k ] = 0 ̲ .
The estimate x ̲ ^ of x ̲ is obtained using the Kalman filter theory:
k + 1 | k = k | k + Q k
K k = k | k H k T [ H k k | k 1 H k T + R k ] 1
x ̲ ^ k | k = x ̲ ^ k | k 1 + K k [ z k h ( x ̲ ^ k | k 1 ) ]
k | k = [ I G k H k ] k | k 1
With H k the Jacobian matrix of h ̲ computed at x ̲ = x ̲ ^ .
The hull differentiation between the anomalies, gross error, and sudden load change, will be based on the deformation of the innovation distribution function [10]. In the following, the Innovation process is introduced. While using an analytical model, the implicit hypothesis is that system behavior is considered stationary by parts.

3.1. The Innovation in the Dynamic State Estimation Process

The system innovation is modeled by
υ ̲ k = z ̲ k h ̲ ( x ̲ ^ k | k 1 )
As one knows, υ k is a white Gaussian process, so that:
E [ υ ̲ k ] = 0 , E [ υ ̲ k υ ̲ j T ] = H k k | k 1 H k T + R k .
Normalizing the innovation process, one obtains:
λ k , i = λ k , i / ρ k , i
ρ k , i 2 = H k k | k 1 H k T + r k , i 2
where υ ̲ k , i is the i t h component of υ ̲ k and H k , i is the i t h line of H k . In this way, one can define the random variable Λ k (the set of all individual innovations υ ̲ k , i ) whose sample is the normalized innovation process, presents a normal distribution with average zero and unitary covariance. Considering the bad data condition, this would not be in fact an innovation, since the innovation probability function would change. Otherwise, considering the stationary level change, that would be the case, since the probability function would be dislocated from zero, but still a normal distribution [10].

3.2. Anomaly Detection

In this work, two types of anomalies are considered, that is, the sudden change of system states and measurement gross errors. Considering the measurement value z ̲ k , in a measurement gross error condition, the predicted state value will present significant errors in the innovation Equation (16) [10,25].
When this anomaly occurs, the innovation distribution ( Λ ) in the normalized form is considerably different from the Λ distribution for the case of stationarity level change, as illustrated in Figure 2, where f( λ ) is the probability distribution function for λ [10]. If those measurements are deleted or corrected, considering the measurement gross error condition, the distribution for Λ is equivalent to the stationary case situation. Otherwise, for the abrupt stationary level changes case, a significant number of elements of the innovation process are affected, and the Gaussian behavior will stay valid, but with an average magnitude, not zero, that is, the symmetry related to the non-zero average value will stay valid, which does not occur for the measurement gross error condition. The detection of sudden load changes is mainly related to the sensitivity of the state variables to the active and reactive power injections of each bus of the system. Thus, the detection of anomalies in a specific node may be more difficult if the system is less sensitive to it, for instance, buses with the smallest loads. On the contrary, it will be more sensitive to variations in buses with more loads or generation. The main factor in such detection feature is the tuning of the state covariance matrix Q k of the Kalman Filter. The elements of the covariance matrix should reflect the impact on the state variables caused by the expected load and generation variations. Considering the previous, in this work anomalies will be detected by the test:
| λ k , i | < α , i = 1 , 2 , n
With α a chosen threshold value (in this work assumed equal 3.0). This value is chosen considering the desired significance level and degrees of freedom. If at least one element of the innovation process is out of the pre-specified region, an anomaly condition is characterized.

3.3. Anomaly Discrimination and Adaption

The anomaly discrimination is in this work done by estimating the deformation of the innovation distribution. This deformation is estimated by the asymmetry index of the innovation distribution related to the average [11]:
γ k = M 3 , k / σ k 3 ,
With M 3 , k the third moment of the distribution, at time k, and σ k 3 the standard deviation of that distribution. The third moment is given by
M 3 , k = E [ Λ k 3 ] 3 μ k E [ Λ k 2 ] + 2 μ k 3
With
σ k 3 = E [ Λ k 2 ] μ k 2 , μ k = E [ Λ k ]
When the innovation deformation tends to zero, the distribution is symmetric and the constant γ k tends to zero, as in the normal distribution. That said, for the measurement bad data condition, the distribution deformation is large, and then γ k becomes large as well, as symmetry is lost. Otherwise, for an abrupt change in the stationary level, the symmetry is not lost and γ k will have a low value. Thus, considering the γ k calculation for the normalized innovation, the anomaly will be discriminated by
  • γ k > | γ m a x | Measurement Gross error situation
  • γ k < | γ m a x | Suddenly state variable change
where, γ m a x is a threshold obtained by experience and at the same time by simulations (in this work assumed equal 2.0). This value is well-chosen considering a level of confidence and degrees of freedom. Toward surpassing the effects of anomalies in the state estimation, an adaptive analytical model is presented. The end goal is to restore the accuracy of the estimation process. The adaptive analytical model consists of two tests:
  • If a gross error is detected (if | γ k , i | α and γ k | γ m a x | ) , the measurement with the largest normalized innovation is not considered, and the estimation process is repeated without this measurement;
  • If a gross error is detected (if | γ k , i | α and γ k < | γ m a x | ) , the estimation considers only the current snapshot of measurements, thus not considering the time relations between states, and the estimation process is repeated using a static estimator:

4. Presentation of Performed Tests and Analysis of Results

This section describes the performed tests and analyses for the presented anomaly detection model. The algorithm behavior for anomalies identification, through the asymmetry index of the distribution curve of the innovation, is presented. The simulation methodology includes random noise in the results of a reference scenario to emulate the measurement set. A load flow calculation obtains the reference scenario, and the measured values are given by (24) [26]. The power system’s temporal aspect is simulated by a sequence of power flow scenarios, each instant with its corresponding loading/generation scenario. Each sample represents a different discretized instant k and the respective measured values. Thus, the simulation consists of a sequence of quasi-stationary conditions. Gross errors are included in the measurements by including larger error values in a specific set of measurements. Anomalies in the power system are simulated as sudden changes in the reference load flow scenario, changing the load/generation values and considering slow dynamic effects between successive instants. Measurements are considered with a sampling rate of one sample per second, a typical period of practical operation and supervisory systems.
z ̲ k , i = z k , i r e f + μ i p r i | z k , i r e f | 3
where z k , i r e f is the reference value for the i t h measurement at time k obtained by the load flow solution, p r i is the metering device precision (assumed as 2% for active and reactive power and 1% for voltage magnitude), and μ i is a random value obtained from the standard Gaussian distribution. A similar equation also provides the diagonal elements of the measurement noise covariance matrix R i i = p r i | z ̲ i , k | / 3 . The process noise covariance matrix was calculated in a similar manner, by associating a percentage of admissible variation for each state variable in two successive instants ( α ), the maximum state ratio in (9), according to a percentage of the previous state value ( x i , k 1 ) , that is
Q i i = α 2 q i = ( α x i , k 1 ) 2

4.1. Conceptual Example and Asymptotic Performance of DSE

This section presents the initial simulation with a 3-bus test system, illustrated on Figure 3, as a conceptual example. The network data and reference loading scenario are provided in Appendix B. This simulation intends to show two main aspects of DSE in power systems, and the importance of supporting the algebraic state variables with the state-space model.
The state vector and measurement vector are defined by the following electrical quantities of the network, with the voltage phase angle at bus #1 as reference ( θ 1 = 0 0 ) :
x = [ V 1 , V 2 , V 3 , θ 2 , θ 3 ] T
z = [ P 1 , P 2 , P 3 , Q 1 , Q 2 , Q 3 , P 1 2 , P 1 3 , P 2 1 , P 2 3 , Q 1 2 , Q 1 3 , Q 2 1 , Q 2 3 , V 1 , V 2 ] T
The measurement model h ( x ) in (9) is built according to the respective monitored electrical quantities. The temporal relations are included in the model as described in (8), that is, with an Identity matrix as a transition matrix. The measurement covariance matrix is a (16 × 16) diagonal matrix with elements obtained according to the respective measurement precision. While the state process covariance matrix is a (5 × 5) diagonal matrix with each element calculated according to (25), using a state rate equal to 0.01%.
In this test, the simulation consists of a stationary load condition under 10 min, with samples acquired at each second, to illustrate the asymptotical behavior of the Kalman Filter in comparison with traditional static state estimation, that is neglecting temporal relations among state variables. Figure 4 illustrates the state variables at bus #2 during the simulation. It illustrates the major conceptual enhancement of dynamic state estimation algorithms, an asymptotical improvement of the estimated values over time. Such advantage is valid under stationary conditions, as shown in the next simulations.

4.2. Performance of the Adaptive Dynamic State Estimator Algorithm

This section presents the application in the IEEE 14 test system [27] with the measurement set illustrated in Figure 5. The simulation consists of an initial stationary load condition under the first 20 s. A load increase at a constant rate of 0.2% per second follows from t = 20 s to t = 60 s, that is, a load ramp for all nodes, followed by more than 10 s of stationary condition. Then a sudden load increase of 5.0% occurs from t = 70 s to t = 90 s at four specific nodes (nodes #3, #4, #13, and #14). The state ratio in this test case was assumed as 0.03% of the previous state variable value.
Figure 6 presents the estimated state variables for node #14 in comparison to the reference load flow during this simulation. Figure 7 presents the estimation error boxplot between the estimated state variables and the reference values from the load flow calculation.
As one can see from the figures, when anomalies occur in the system, the estimation process deteriorates, especially for the phase angle. This is though expected, as, during such events, the hypothesis that the state transition matrix is equal to an identity matrix becomes false, that is, there is a change in the stationarity level of the system. Moreover, the proposed adaptation strategy successively identifies and adapts the estimator to each change on the stationarity level and in the presence of gross errors.
The adaptation strategy is triggered during the load ramp, at time steps #26s, #35s, #44s, and #54s. This is because the load ramp is a persistent anomaly, at that time during 40 s, which continuously changes the stationarity level of the system behavior during that interval. That said, the sudden load change (load increase and load decrease) are sudden changes that occur in a single instant, at time step #70s and #90s. In this case, the anomaly needs to be suppressed only once, in the instant of the load change, and the dynamic state estimator will return to the correct and new stationary level in a normal way. Figure 8 illustrates maximum absolute innovation value from all measurements and the calculated asymmetry index, in each instant.
One can see that, after the inclusion of measurement gross error, the asymmetry index increases. Besides, when of the sudden load increase of 5%, the asymmetry index keeps its value below the threshold, as expected.

4.3. Effect of Process Noise and Anomalies Size

This case study considered the reference load flow scenario of the IEEE 30 bus test system [1]. A random load variation was included with a uniform variation of 0.5% in all loads in this case. Still, a sudden load increase of 5.0% occurs from t = 40 s to t = 90 s at six specific nodes (#3, #4, #15, #19, #21, and #30). The state ratio in this test case was also assumed as 0.03% of the previous state variable value. Initially, different levels and locations of gross errors were inserted in the simulations. Different levels of gross errors were added in different measurements during the simulations. Figure 9 illustrates the normalized innovation and the asymmetry index for this simulation. It shows that anomaly discrimination, based on the asymmetry index and the innovation level, enables an automatic suppression of different levels and locations of gross errors. An estimation accuracy like the previous section was obtained.
This second test also evaluated different sizes of sudden load changes and the detection by the asymmetry index of such change on the stationarity level. Different values of additional loads were included to the same six nodes (#3, #4, #15, #19, #21, and #30), that comprise 20% of the total net load. Table 1 presents the anomaly discrimination indexes, the maximum normalized innovation, and the asymmetry index, for each level of sudden load change, and different percentages for the state ratio. The method’s ability to detect such changes on the stationarity level depends on the assumed state covariance values. For instance, the state covariance matrix increase pushes the performance of the dynamic state estimator towards the static state estimator. Thereby, it becomes insensitive to changes on the stationarity level and neglects the additional information from the state-space model. It is noteworthy that the asymmetry index kept inside the accepted variation during such sudden load increases, i.e., it kept its value below the threshold value of 2.0.
Finally, different levels of random load variation were also evaluated. Table 2 shows the mean absolute error, among all state variables, for different load variation levels and different percentages of state ratios
As can be seen, the DSE increases its performance under stationary conditions. Such performance is directly related to proper covariance matrix tuning. Besides, by increasing such covariance, the estimation reaches an error plateau given by a static state estimator’s performance. It is noteworthy that load variations of 2.0% are very large for the time span of the SCADA sampling rate, a few seconds. Such size of load change is more prone to be associated with sudden changes type of anomaly.

4.4. Effect of Different System Anomalies

To further explore and illustrate the asymmetry index, this simulation considers the IEEE 118 test system, a medium-size transmission system [27]. The simulation comprises one hundred sequential seconds, with the inclusion of 0.1% of random load variation in all nodes. This test scenario also considers a loss of generation at node #100, reducing 252 MW (6.1% of net load) at t = 40 s. The state ratio in this test case was assumed as 1.0% of the previous state variable value. Figure 10 illustrates the normalized innovation and the asymmetry index, demonstrating the automatic suppression of both gross errors and sudden changes while keeping accuracy throughout the estimation.
Figure 11 illustrates the cumulative probability plot of the normalized innovation for this test during different anomalies. It clearly shows the asymmetry present in the normalized innovation sample in the presence of gross errors, captured by the proposed anomaly discrimination strategy. Despite the results during the initial stationary event, as well as during the loss of generation, show a heavy-tailed distribution, it has more symmetry when compared to the presence of gross errors.
To demonstrate the effect of tuning and the influence of the sensitivity of the state variables, different levels of sudden load changes in individual buses were included across the power system. Two different approaches for tuning the state covariance matrix were evaluated, the first with a constant percentage of possible state variation, as in the previous simulations, and the second by choosing individual values of the state ratio parameter for each state variable. This later choice seeks to capture a more sensitive detection for anomalies that may occur in different parts of the system according to the respective state sensitivity.
Figure 12 illustrates a map with the detection level of sudden load change required to trigger the anomaly detection method in these two tuning approaches. The results are also present in Figure 13 with the obtained load change detection levels in each bus of the system. As it can be seen, the sensitivity of buses with less load requires larger levels of sudden changes, sometimes much more than 40%, to be captured by the anomaly detection method. However, as shown in the second simulation, by changing the tuning strategy for the state covariance was possible to increase such sensitivity accurately detect individual load changes less than 10% in most of the buses. Thus, showing the covariance matrix tuning as an important aspect of dynamic state estimators and their response against different anomalies.

4.5. Computational Aspects

Finally, the computational aspects of the previous simulations are presented in Table 3. The tests were performed using a microcomputer with a Core i7 3.60 GHz, 16GB RAM with C programming language. The computational times only indicate the algorithm feasibility, and there is a large space for improvements, such as exploring sparsity treatments, improving data structures, and employing parallel computing techniques.

5. Conclusions

This paper presented a family of Kalman filter modeling for the power system dynamic state estimation problem. The presented formulation encompasses different piecewise stationary levels by discriminating different anomalies in the power system, such as gross errors or sudden load changes. The presented method automatically adapts itself to surpass the harmful effects of such anomalies while maintaining accuracy. Simulation results in IEEE transmission systems demonstrate the performance under dynamic anomalies, like different scenarios of gross errors, load ramps, sudden load changes, and generator contingency. Incorporating a dynamic perspective into the state estimation enables a sensible improvement on estimation accuracy when the system keeps a stationary level. Besides, by incorporating the asymmetry index into the innovation analysis, it is possible to track significant changes in the system’s dynamics accurately and adapt the estimator response. The state covariance tuning significantly impacts the Kalman filter performance, pushing its results between the static estimator and a smooth function. Future works comprise new Kalman filter tuning strategies based on power flow sensitivity, exploring the effects of bias into the state-space, and enhancing computational performance.

Author Contributions

Conceptualization, A.S.B.; methodology, A.S.B.; software, A.S.B. and J.A.D.M.; validation, A.S.B. and J.A.D.M.; formal analysis, N.G.B.; investigation, A.S.B. and J.B.A.L.J.; resources, A.S.B.; data curation, J.A.D.M.; writing—original draft preparation, A.S.B., N.G.B. and J.A.D.M.; writing—review and editing, A.S.B. and N.G.B.; visualization, J.A.D.M. and J.B.A.L.J.; supervision, N.G.B.; project administration, A.S.B.; funding acquisition, A.S.B. All authors have read and agreed to the published version of the manuscript.

Funding

This work was supported by NSF grant ECCS-1809739.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

Not applicable.

Conflicts of Interest

The authors declare no conflict of interest.

Appendix A. Power Systems Nonlinear Measurement Model

The measurement model for power system dynamic state estimation, with m measurements and n state variables, comprises nonlinear equations that relate the measured electrical quantities to the state variables. The state vector (n × 1) consists of the complex nodal voltages, on polar coordinates, in all buses of the network also referred to as the algebraic state of the electric power network. The measurement vector (m × 1) comprises electrical quantities monitored across the power system, typically by SCADA systems, which provide voltage magnitudes, active and reactive power flows, and injections, at different locations.
x = [ V , θ ] T
z = [ P i , Q i , P i j , Q i j , V i ] T
where V and θ are the voltage magnitude and voltage phase vectors in all nodes of the power system, P i and Q i are the active and reactive power injections measured at the respective node i, and P i j and Q i j are the active and reactive power flows measured at the branch between nodes i and j. Taking a power system with nbus the number of nodes, the total state variables is n = 2nbus − 1, since one of the phase angles is taken as a reference. The branch-bus model represents the components of the power grids and relates the monitored electrical quantities to the state variables, as illustrated in Figure A1. Hence, power flow equations represent most of the measured values according to the following:
P i j = V i 2 g i j a i j V i V j ( g i j cos ( θ i j ϕ i j ) + b i j sin ( θ i j ϕ i j ) )
Q i j = V i 2 ( b i j + b i j s h ) a i j V i V j ( g i j s i n ( θ i j ϕ i j ) b i j c o s ( θ i j ϕ i j ) ) .
where g i j and b i j are the series conductance and susceptance of the branch component (obtained by the series resistance and reactance r i j and x i j ), b i j s h is the shunt susceptance of the branch component (obtained by its nominal reactive loading), a i j is the off-nominal transformer relation and ϕ i j its phase-shifting.
Figure A1. Branch-bus model for power system components connected between two nodes i-j. The model can represent transmission lines and transformers.
Figure A1. Branch-bus model for power system components connected between two nodes i-j. The model can represent transmission lines and transformers.
Energies 14 06787 g0a1
For power injections, the summation of adjacent power flows provides the mathematical model, as described below.
P i = j Ω i P i j
Q i = V i 2 b i s h + j Ω i Q i j
where b i s h is the shunt susceptance for elements connected directly in bus i, and Ω i is the set of adjacent nodes to bus i. Finally, voltage magnitude measurements are directly related to the corresponding state variable of the bus being monitored.
The Jacobian matrix elements are obtained by the derivatives of the above equations and respective structure of the measurement vectors, as described in (). The equations for the derivatives are directly obtained by differentiating the above equations, and further details can be refereed in [1].
H ( x ) = h ( x ) x = P i ( x ) V P i ( x ) θ Q i ( x ) V Q i ( x ) θ P i j ( x ) V P i j ( x ) θ Q i j ( x ) V Q i j ( x ) θ V i ( x ) V 0 ,

Appendix B. Three-Bus Test System Data

The data for the 3-bus test system is provided below. A power flow calculation with the described loading scenario provides the reference case for the simulations.
Table A1. Network data for nodes on the reference case on the 3-bus test system.
Table A1. Network data for nodes on the reference case on the 3-bus test system.
VoltageActiveReactiveActiveReactive
Bus IDTypeSetpointGenerationGenerationLoadLoad
(p.u)(MW)(MVAr)(MW)(MVAr)
Bus 1 V θ 1.02040.615.6--
Bus 2PV1.00057.014.4--
Bus 3PQ---95.019.0
Table A2. Network data for branches on the reference case on the 3-bus test system (base power 100 MVA and base voltage 69 kV).
Table A2. Network data for branches on the reference case on the 3-bus test system (base power 100 MVA and base voltage 69 kV).
Bus FromBus ToVoltage
Ratio α ij
Resistance
(%)
Reactance
(%)
Nominal Reactive
Loading (MVAr)
Bus 1Bus 211.955.905.30
Bus 1Bus 215.4022.304.90
Bus 2Bus 314.7019.804.40

References

  1. Bretas, A.S.; Bretas, N.G.; London, J.B.; Carvalho, B.E. Cyber-Physical Power Systems State Estimation; Elsevier: Amsterdam, The Netherlands, 2021. [Google Scholar]
  2. Bretas, A.S.; Bretas, N.G.; Carvalho, B.; Baeyens, E.; Khargonekar, P.P. Smart grids cyber-physical security as a malicious data attack: An innovation approach. Electr. Power Syst. Res. 2017, 149, 210–219. [Google Scholar] [CrossRef]
  3. Zhao, J.; Netto, M.; Huang, Z.; Yu, S.S.; Gomez-Exposito, A.; Wang, S.; Kamwa, I.; Akhlaghi, S.; Mili, L.; Terzija, V.; et al. Roles of dynamic state estimation in power system modeling, monitoring and operation. IEEE Trans. Power Syst. 2020, 36, 2462–2472. [Google Scholar] [CrossRef]
  4. Rodriguez Paz, M.; Ferraz, R.; Bretas, A.S.; Leborgne, R. System unbalance and fault impedance effect on faulted distribution networks. Comput. Math. Appl. 2010, 60, 1105–1114. [Google Scholar] [CrossRef] [Green Version]
  5. Zhao, J.; Gómez-Expósito, A.; Netto, M.; Mili, L.; Abur, A.; Terzija, V.; Kamwa, I.; Pal, B.; Singh, A.K.; Qi, J.; et al. Power system dynamic state estimation: Motivations, definitions, methodologies, and future work. IEEE Trans. Power Syst. 2019, 34, 3188–3198. [Google Scholar] [CrossRef]
  6. Kalman, R.E. A new approach to linear filtering and prediction problems. Trans. ASME—J. Basic Eng. 1960, 35–45. [Google Scholar] [CrossRef] [Green Version]
  7. Debs, A.S.; Larson, R.E. A dynamic estimator for tracking the state of a power system. IEEE Trans. Power Appar. Syst. 1970, 89, 1670–1678. [Google Scholar] [CrossRef]
  8. Schweppe, F.; Masiello, R. A tracking static state estimator. IEEE Trans. Power Appar. Syst. 1971, 90, 1025–1033. [Google Scholar] [CrossRef]
  9. Nishiya, K.; Hasegawa, J.; Koike, T. Dynamic state estimation including anomaly detection and identification for power systems. Proc. IEEE 1982, 129, 192–198. [Google Scholar] [CrossRef]
  10. Falcao, D.; Cooke, P.; Brameller, A. Power system tracking state estimation and bad data processing. IEEE Trans. Power Appar. Syst. 1982, PAS-101, 325–333. [Google Scholar] [CrossRef]
  11. Bretas, N. An iterative dynamic state estimation and bad data processing. Int. J. Electr. Power Energy Syst. 1989, 11, 70–74. [Google Scholar] [CrossRef]
  12. Rousseaux, P.; Van Cutsem, T.; Liacco, T.D. Whither dynamic state estimation? Int. J. Electr. Power Energy Syst. 1990, 12, 104–116. [Google Scholar] [CrossRef]
  13. Sarri, S.; Zanni, L.; Popovic, M.; Le Boudec, J.Y.; Paolone, M. Performance assessment of linear state estimators using synchrophasor measurements. IEEE Trans. Instrum. Meas. 2016, 65, 535–548. [Google Scholar] [CrossRef]
  14. Fan, L.; Wehbe, Y. Extended Kalman filtering based real-time dynamic state and parameter estimation using PMU data. Electr. Power Syst. Res. 2013, 103, 168–177. [Google Scholar] [CrossRef]
  15. Valverde, G.; Terzija, V. Unscented Kalman filter for power system dynamic state estimation. IET Gener. Transm. Distrib. 2010, 5, 29–37. [Google Scholar] [CrossRef]
  16. Sharma, A.; Srivastava, S.C.; Chakrabarti, S. A cubature Kalman filter based power system dynamic state estimator. IEEE Trans. Instrum. Meas. 2017, 66, 2036–2045. [Google Scholar] [CrossRef]
  17. Massignan, J.A.; London, J.B.; Miranda, V. Tracking Power System State Evolution with Maximum-correntropy-based Extended Kalman Filter. J. Mod. Power Syst. Clean Energy 2020, 8, 616–626. [Google Scholar] [CrossRef]
  18. Zhou, N.; Meng, D.; Huang, Z.; Welch, G. Dynamic state estimation of a synchronous machine using PMU data: A comparative study. IEEE Trans. Smart Grid 2014, 6, 450–460. [Google Scholar] [CrossRef]
  19. Do Coutto Filho, M.B.; De Souza, J.C.S.; Guimaraens, M.A.R. Enhanced bad data processing by phasor-aided state estimation. IEEE Trans. Power Syst. 2014, 29, 2200–2209. [Google Scholar] [CrossRef]
  20. Durgaprasad, G.; Thakur, S. Robust dynamic state estimation of power systems based on M-estimation and realistic modeling of system dynamics. IEEE Trans. Power Syst. 1998, 13, 1331–1336. [Google Scholar] [CrossRef]
  21. Zhao, J.; Mili, L. A framework for robust hybrid state estimation with unknown measurement noise statistics. IEEE Trans. Ind. Inform. 2017, 14, 1866–1875. [Google Scholar] [CrossRef]
  22. Rouhani, A.; Abur, A. Linear phasor estimator assisted dynamic state estimation. IEEE Trans. Smart Grid 2016, 9, 211–219. [Google Scholar] [CrossRef]
  23. Zhao, J. Dynamic State Estimation With Model Uncertainties Using H Extended Kalman Filter. IEEE Trans. Power Syst. 2017, 33, 1099–1100. [Google Scholar] [CrossRef]
  24. Wang, Y.; Sun, Y.; Dinavahi, V.; Cao, S.; Hou, D. Adaptive robust cubature Kalman filter for power system dynamic state estimation against outliers. IEEE Access 2019, 7, 105872–105881. [Google Scholar] [CrossRef]
  25. Do Coutto Filho, M.B.; de Souza, J.C.S. Forecasting-aided state estimation—Part I: Panorama. IEEE Trans. Power Syst. 2009, 24, 1667–1677. [Google Scholar] [CrossRef]
  26. Castillo, M.R.; London, J.B.; Bretas, N.G.; Lefebvre, S.; Prévost, J.; Lambert, B. Offline detection, identification, and correction of branch parameter errors based on several measurement snapshots. IEEE Trans. Power Syst. 2010, 26, 870–877. [Google Scholar] [CrossRef]
  27. IEEE Transmission System Test Cases: IEEE14, IEEE30 and IEEE118. IEEE Test Systems. Available online: https://labs.ece.uw.edu/pstca/ (accessed on 15 August 2021).
Figure 1. Signal model with external input known.
Figure 1. Signal model with external input known.
Energies 14 06787 g001
Figure 2. Distribution function of the New Information: (a) Dashed line for Normal Distribution and Solid line for Sudden load change situation; (b) Measurement Gross Error case.
Figure 2. Distribution function of the New Information: (a) Dashed line for Normal Distribution and Solid line for Sudden load change situation; (b) Measurement Gross Error case.
Energies 14 06787 g002
Figure 3. The 3-bus transmission system and respective measurement set for the simulation.
Figure 3. The 3-bus transmission system and respective measurement set for the simulation.
Energies 14 06787 g003
Figure 4. Example of estimated state variables (voltage magnitude and phase angle) at bus #2 during the simulations, with the dynamic state estimator and with a static state estimator.
Figure 4. Example of estimated state variables (voltage magnitude and phase angle) at bus #2 during the simulations, with the dynamic state estimator and with a static state estimator.
Energies 14 06787 g004
Figure 5. The IEEE 14 bus transmission system and respective measurement set for the simulation.
Figure 5. The IEEE 14 bus transmission system and respective measurement set for the simulation.
Energies 14 06787 g005
Figure 6. Example of estimated state variables (voltage magnitude and phase angle) at bus 14 during the simulations with the adaptive dynamic state estimator.
Figure 6. Example of estimated state variables (voltage magnitude and phase angle) at bus 14 during the simulations with the adaptive dynamic state estimator.
Energies 14 06787 g006
Figure 7. Boxplot of estimation error in all state variables at each instant of the simulation: (a) conventional Kalman Filter; (b) proposed Adaptive Dynamic State Estimator.
Figure 7. Boxplot of estimation error in all state variables at each instant of the simulation: (a) conventional Kalman Filter; (b) proposed Adaptive Dynamic State Estimator.
Energies 14 06787 g007
Figure 8. Maximum normalized innovation and asymmetry index for the IEEE14 in the simulations.
Figure 8. Maximum normalized innovation and asymmetry index for the IEEE14 in the simulations.
Energies 14 06787 g008
Figure 9. Maximum normalized innovation and asymmetry index for the IEEE30 in the simulations.
Figure 9. Maximum normalized innovation and asymmetry index for the IEEE30 in the simulations.
Energies 14 06787 g009
Figure 10. Maximum normalized innovation and asymmetry index for the IEEE118 in the simulation.
Figure 10. Maximum normalized innovation and asymmetry index for the IEEE118 in the simulation.
Energies 14 06787 g010
Figure 11. Probability plots for the normalized innovation in different instants: (left) in a stationary condition; (middle) with multiple gross errors; (right) during a loss of generation event.
Figure 11. Probability plots for the normalized innovation in different instants: (left) in a stationary condition; (middle) with multiple gross errors; (right) during a loss of generation event.
Energies 14 06787 g011
Figure 12. Locational aspect of sudden load change detection with two different state covariances tuning: the first with a fixed state ratio for all buses and the second with an individual level of state ratio for each bus.
Figure 12. Locational aspect of sudden load change detection with two different state covariances tuning: the first with a fixed state ratio for all buses and the second with an individual level of state ratio for each bus.
Energies 14 06787 g012
Figure 13. Amount of load change, in the percentage of the nominal value, in each individual node, detected as an anomaly for two different state covariances tuning.
Figure 13. Amount of load change, in the percentage of the nominal value, in each individual node, detected as an anomaly for two different state covariances tuning.
Energies 14 06787 g013
Table 1. Dynamic State Estimator anomaly discrimination during different levels of sudden load change for the IEEE30 system with different state ratios for the state covariance matrix.
Table 1. Dynamic State Estimator anomaly discrimination during different levels of sudden load change for the IEEE30 system with different state ratios for the state covariance matrix.
State RatioDiscrimination IndexSudden Load at Nodes: #3, #4, #15, #19 and #30
( α ) 1.00%2.00%3.00%5.00%10.00%
α = 0.01%Max.Norm.Innovation2.23952.58902.96697.064814.7944
Asymmetry Index−0.75640.43540.24640.13180.1314
α = 0.1%Max.Norm.Innovation2.24721.69491.68573.13966.1013
Asymmetry Index−1.55510.52880.44430.64520.4482
α = 1.0%Max.Norm.Innovation2.21651.34831.33672.08013.2864
Asymmetry Index−1.60050.5130−0.06420.59331.0896
Table 2. Mean absolute error of the estimation for different levels of load variation for the IEEE30 system with different state ratios for the state covariance matrix.
Table 2. Mean absolute error of the estimation for different levels of load variation for the IEEE30 system with different state ratios for the state covariance matrix.
Load VariationMaximum State Ratio in the State Covariance Matrix
α = 0.00% α = 0.01% α = 0.1% α = 1.0% α = 10.0%
0.00%8.9252 × 10−62.6840 × 10−57.9222 × 10−51.8912 × 10−42.3610 × 10−4
0.50%4.4572 × 10−42.9600 × 10−42.1010 × 10−41.9453 × 10−42.3684 × 10−4
1.00%12.1285 × 10−45.8336 × 10−43.7315 × 10−42.0537 × 10−42.3761 × 10−4
2.00%61.1972 × 10−411.6417 × 10−47.1472 × 10−42.3882 × 10−42.3922 × 10−4
Table 3. Mean processing time for different stages of the proposed DSE.
Table 3. Mean processing time for different stages of the proposed DSE.
Test SystemPrediction StepEstimation StepAnomaly Discrimination
IEEE14<1.0 ms<1.0 ms1.9 ms
IEEE30<1.0 ms16.0 ms9.7 ms
IEEE11815.2 ms360.0 ms18.4 ms
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Bretas, A.S.; Bretas, N.G.; Massignan, J.A.D.; London Junior, J.B.A. Hybrid Physics-Based Adaptive Kalman Filter State Estimation Framework. Energies 2021, 14, 6787. https://doi.org/10.3390/en14206787

AMA Style

Bretas AS, Bretas NG, Massignan JAD, London Junior JBA. Hybrid Physics-Based Adaptive Kalman Filter State Estimation Framework. Energies. 2021; 14(20):6787. https://doi.org/10.3390/en14206787

Chicago/Turabian Style

Bretas, Arturo S., Newton G. Bretas, Julio A. D. Massignan, and João B. A. London Junior. 2021. "Hybrid Physics-Based Adaptive Kalman Filter State Estimation Framework" Energies 14, no. 20: 6787. https://doi.org/10.3390/en14206787

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop