Robust Power System State Estimation Method Based on Generalized M-Estimation of Optimized Parameters Based on Sampling

Abstract

Robustness is an important performance index of power system state estimation, which is defined as the estimator’s capability to resist the interference. However, improving the robustness of state estimation often reduces the estimation accuracy. To solve this problem, this paper proposes a power system state estimation method for generalized M-estimation of optimized parameters based on sampling. Compared with the traditional robust state estimator, the generalized M-estimator based on projection statistics improves the robustness of state estimation, and the proposed optimized parameter determination method improves the overall accuracy of state estimation by appropriately adjusting its robustness. Considering different degrees of non-Gaussian distributed measurement noises and bad data, the estimation accuracy the proposed method is demonstrated to be up to 23% higher than the traditional generalized M-estimator through MATLAB simulations in IEEE 14, 118 bus test systems, and Polish 2736 bus system.

1. Introduction

Modern power systems need to grasp the real-time operation of the power system comprehensively and quickly. It can more accurately analyze and predict the system’s operation trends to ensure the power system’s economy and safety. Power system state estimation is an effective means to detect the real-time operation data of the power system which was first proposed by the Schweppe in 1970 [1]. Abundant research has been conducted on power system state estimation to optimize the real-time redundant measurement data of the power system, to realize the real-time and reliable monitoring of the power system, and to ensure the safe and reliable operation of the system [2,3,4,5,6,7,8].

The weighted least squares (WLS) method has been widely adopted since 1970, with the advantages of good convergence performance and accurate estimation results. Its disadvantage is that it is not sufficiently robust to measured noise that follows the Gaussian distribution without bad data interference. In the real power network, not all measurement noise can be represented strictly by a Gaussian distribution, and occasionally significant errors occur due to various known telemetry noises or faults. There may also be other types of outliers that strongly affect the estimated state but may or may not carry bad data. All of these factors can potentially affect the estimator’s estimation accuracy. A series of studies have been conducted to improve the robustness of state estimation, which is summarized in

Table 1. Reference summary.

Authors	Origin	Purpose	Advantages and Disadvantages
L. Mili, M. G. Cheniae, N. S. Vichare and P. J. Rousseeuw	USA	To describe a fast and robust method for identifying the leverage points [9].	The method is very fast and compatible with real-time applications, but it does not apply to all forms of lever points.
J. Zhao and L. Mili	USA	To develop a robust dynamic state estimator of a cyber-physical system [10].	The H-infinity filter is able to handle large system uncertainties as well as suppress outliers, but the estimation efficiency of this method is low.
M. B. Djukanovic, M. H. Khammash and V. Vittal	USA	To present a framework for robust stability assessment in multimachine power systems [11].	The proposed method significantly reduces computational complexity and at the same time preserves the accuracy in predicting stability robustness.
Z. Lyu, H. Wei, X. Bai, D. Xie, L. Zhang and P. Li	CHN	To propose an norm estimator [12].	The proposed estimator has high computational efficiency and robustness.
E. Kyriakides, S. Suryanarayanan and G. T. Heydt	USA	To demonstrate the Huber function technique in a power engineering application [13].	This technique reduces large residuals but not accuracy.
M. Göl and A. Abur	TR	To develop a PMU placement strategy [14].	This method can improve the stability and accuracy of estimation.
M. Netto, J. Zhao and L. Mili	USA	To develop a robust extended Kalman filter [15].	The robust extended Kalman filter exhibits good tracking capabilities under Gaussian process and observation noise while suppressing observation outliers, even in position of leverage. However, it presents poor performance under non-Gaussian noise.
I. Akingeneye, J. Wu and J. Yang	USA	To develop PMU placement algorithms to improve the power grid state estimation [16].	The performance of the low complexity algorithms approach that of the exhaustive search algorithm, but with a much lower complexity.
G. Wang, G. B. Giannakis and J. Chen	USA	To put forward a novel LAV estimator leveraging recent algorithmic advances in composite optimization [17].	The algorithm efficiently deals with the non-convexity and non-smoothness of LAV-based PSSE, but it relies on solving a sequence of convex quadratic subproblems.
M. Huang, Z. Wei, G. Sun and H. Zang	CHN	To propose a hybrid SE for distribution systems [18].	The estimator method provides more reliable estimation results with a limited number of SCADA measurements, while biased estimated results can exist if some buses are far away from the measuring points.
C. H. Ho, H. C. Wu, S. C. Chan and Y. Hou	CHN	To present a robust statistical approach [19].	The proposed approach outperforms conventional approaches using the ADMM with L1 outlier detection in state estimation accuracy and convergence speed.
J. Zhao, M. Netto and L. Mili	USA	To develops a robust iterated extended Kalman filter based on the generalized maximum likelihood approach [20].	GM-IEKF can achieve both robustness and statistical efficiency, but its vulnerability to system parameter and topology errors.

:

The problem of existing robust state estimators is that robustness of the state estimator is usually achieved at the expense of estimation accuracy. The reason is that the existing robust state estimators use fixed parameter settings, which cannot adjust the state estimator’s robustness and estimation accuracy in different scenarios. In [21], it is indicated that the measurement noise in PMU is likely to be non-Gaussian, which leads to the increased probability of outliers (measurements that significantly departs from their true values). Moreover, the types, parameters, and proportions of the non-Gaussian measurement distribution are different in different systems. Therefore, existing robust state estimators will be less robust when bad data and outliers are more frequent and will have insufficient accuracy when bad data and outliers are less frequent. To overcome this problem, this paper first evaluates the impact of bad data and non-Gaussian measurement noises on existing generalized M-estimation parameters based on extensive Monte Carlo simulations. Finally, this paper proposes a generalized M-estimator of optimized parameters based on sampling which can adaptively select the appropriate parameters according to the probability distribution of the system measurement noise, accurately identify the outliers, and significantly reduce the effects of non-Gaussian measurement noise and bad data.

The contributions of this paper are summarized as follows:

• Demonstrate that higher robustness does not necessarily improve the estimation accuracy of the state estimator, and the best accuracy can be achieved if the robustness is tuned at an appropriate level;
• Propose a new robust power system state estimation method that can adaptively tune its robustness according to different levels of non-Gaussian distributed measurement noise and bad data.

The paper is structured as follows. Section 2 introduces the existing robust state estimation methods. Section 3 describes the drawbacks of existing robust state estimation methods. Section 4 presents the proposed generalized M state estimator of optimal parameters based on sampling. Section 5 presents and discusses the simulation results. Finally, Section 6 concludes this paper and identifies opportunities for future developments.

2. Existing Robust State Estimators for Power System Estimation

This section presents two widely used robust state estimators for power system state estimation, M-estimator (including Huber estimator and LAV estimator) and generalized M-estimator, which is improved based on the M-estimator.

2.1. M-Estimation of the Static State Estimation Method

The concept of M-estimation was first used by Huber for robust estimation of distribution centers and subsequently generalized to regression [22]. In general, an M-estimator is a maximum-likelihood estimator. It minimizes an objective function that is expressed as a function ρ(r) for measuring the residue, according to the constraints given by the measurement equation:

Objective function:

$${\displaystyle \sum _{i=1}^m\rho (r_i)}$$

(1)

Constraint condition:

$$z=\mathrm{h}\left(x\right)+r$$

(2)

where ρ(r_i) is a selected function that can measure the residual r_i, z is the measurement vector, x is a state vector and, h(x) is a measurement function.

The Huber estimator objective function is expressed as:

$$J(x)={\displaystyle \sum _{i=1}^m{\omega }_i^2\rho (r_i)}$$

(3)

where the Huber function, ρ(r_i), can be defined as:

$$\rho (r_i)=\{\begin{array}{l}{r}_i{}^2/2\hspace{1em}\hspace{1em}\hspace{1em}\hspace{0.17em}\mathrm{for}\left|r_i\right|\le \beta \\ \beta \left|r_i\right|-{\beta }^2/2\hspace{1em}\mathrm{for}\left|r_i\right|>\beta \end{array}.$$

(4)

The first partial derivative ψ(r_i) of r_i is expressed as:

$$\psi \left(r_i\right)=\{\begin{array}{c}{r}_i,\hspace{1em}\hspace{1em}\hspace{1em}\left|r_i\right|<\beta \\ \beta \ast sign\left(r_i\right),else\end{array}.$$

(5)

The least absolute value (LAV) function can be defined as:

$$\rho \left(r_i\right)=\left|r_i\right|$$

(6)

where r_i is the normalized residual, the parameter β is a fixed value, generally set to 1.5. In fact, if β is set to infinite, the Huber estimator is equal to WLS estimator parameters, which has the highest accuracy. If β is set to 0, the Huber estimator turns into a LAV estimator, which has the highest robustness. Therefore, the Huber estimator can be regarded as an estimator which strikes a balance between the advantages of the WLS estimator and the LAV estimator.

The loss functions of the Huber and the LAV estimator are shown in Figure 1:

2.2. Generalized M-Estimator

The external interference of the system causes outliers. There are multiple types of external interference. For example, there can be errors in sampling, such as recorded bias, calculation errors, etc. It can also be caused by various accidental and abnormal factors. The presence of outliers affects the accuracy of state estimation. Although the Huber estimator and the LAV estimator have better robustness compared to the WLS estimator, it is still obviously affected by outliers. Therefore, early robust estimators suppress the effect of outliers by reducing their weights according to their normalized residuals.

The normalized residual is based on Mahalanobis distance, which is the weighted average distance between each point. The disadvantage is that the average point is significantly affected by outliers, leading to an inaccurate judgment of outliers. In contrast, the median is hardly affected by the outliers. However, the variance of each point to the median cannot be directly defined. It is more reasonable to use the projection statistics method. The idea of the projection statistics method is to study those coming from the coordinated median M, and through each data point direction v. Formally, we have v = l_i − M, where:

$$M={\left[\underset{j}{median}{l}_{j1},\cdots ,\underset{j}{median}{l}_m\right]}^T$$

(7)

where j is obtained from 1 to m. The resulting distance calculated with (7) will be referred to as the projection statistics, or the PS_i for short. Then, the i-th calculated projection statistic, PS_i, is compared to a given threshold. The labeled outliers are then degraded using the following weight function: ${\varpi }_i^2=\min (1,d^2/P{S}_i^2)$, where d = 1.5 is set to produce good statistical efficiency.

The generalized M-estimator combines outlier detection and the Huber loss function, whose objective function is expressed as:

$$J(x)={\displaystyle \sum _{i=1}^m{\varpi }_i^2\rho (r_s)}.$$

(8)

To solve the objective function, we calculated the first derivative about x:

$$\frac{\partial J(x)}{\partial x}={\displaystyle \sum _{i=1}^m-\frac{{\varpi }_i{a}_i}{s}\psi (r_{si})}=0$$

(9)

where s = 1.4826 median_i |r_i|.

It is clear that this is a set of equations which can be solved by iterated re-weighted least square (IRLS) algorithm [23].

3. Disadvantages of Existing Robust State Estimators for Power System State Estimation

Existing robust estimators mainly consider the suppression performance on bad data, while the assessment of their outlier suppression performance caused by non-Gaussian measurement noise is insufficient. This might result in a decreased estimation accuracy of the state estimator in the presence of non-Gaussian measurement noise. The two most common types of non-Gaussian measurement noise are bimodal Gaussian distribution and Laplace distribution, which are introduced in detail in Section 3.1. Section 3.2 discusses the impact of non-Gaussian measurement noise on the performance of existing robust estimators. Section 3.3 presents the purpose of this study.

3.1. Non-Gaussian Distributed Measurement Noises

3.1.1. Bimodal Gaussian Distribution

In practice, the errors of voltage and current measurements might follow bimodal Gaussian mixture (BGM) distribution [24]. The probability density function can be obtained by the superposition of two Gaussian probability density functions:

$$PD{F}_{BGM}={\displaystyle \prod _{i=1}^2{\omega }_i}{f}_{N({\widehat{\mu }}_{i,}{\widehat{\sigma }}_i^2)}(y)$$

(10)

where f_N(y) represents the Gaussian probability density function, N(µ,σ²) represents the normal distribution, and ω represents the weights corresponding to the combination of Gaussian components. The symbol “^” represents the estimated quantity, and the subscript i represents the i-th Gaussian component combination. The weight of each Gaussian component combination is the product of all the Gaussian component weights involved in that combination, satisfying:

$${\omega }_1+{\omega }_2=1.$$

(11)

For the bimodal mixed Gaussian distribution in the method, the mean and variance are as follows:

$${\mu }_{12}={\omega }_1{\mu }_1+{\omega }_2{\mu }_2$$

(12)

$${\sigma }_{12}^2={\omega }_1{\sigma }_1^2+{\omega }_2{\sigma }_2^2+{\omega }_1{\omega }_2{({\mu }_1-{\mu }_2)}^2$$

(13)

where ω represents the weight, μ represents the mean value, σ represents the variance, subscripts 1 and 2 indicate the components 1 and 2 of the BGM, respectively, and subscript 12 indicates the resulting BGM distribution. In the proposed method, μ₁₂ is set to 0, which results in

$$\frac{{\mu }_1}{{\mu }_2}=-\frac{{\omega }_2}{{\omega }_1}=k.$$

(14)

Combine (12) and assume σ₁ = σ₂ = σ, we have,

$${\mu }_1=-k{\omega }_2\sigma$$

(15)

$${\mu }_2=k{\omega }_1\sigma .$$

(16)

Therefore, the distribution of the BGM measurement noise can be changed by adjusting the parameter, k.

3.1.2. Laplace Distribution

The Laplace distribution is also called bi-exponential distribution because it can be seen as the combination of two exponential distributions at different positions. The probability density function of the Laplace distribution is:

$$f(x|\mu ,b)=\frac{1}{2b}\exp (-\frac{\left|x-\mu \right|}{b})=\frac{1}{2b}\{\begin{array}{c}\exp (-\frac{x-\mu }{b}),ifx\ge \mu \\ \exp (-\frac{\mu -x}{b}),ifx<\mu \hspace{1em}\end{array}$$

(17)

where μ is the position parameter, b > 0 is the scaling parameter. If μ = 0, then the positive half happens to be an exponential distribution of scale 1/2. The difference between Laplace distribution and Gaussian distribution is that the Gaussian distribution represents the square of the difference relative to the mean, while the Laplace distribution is represented by the absolute value relative to the difference [25]. Thus, the tail Laplace distribution is much flatter than that of the Gaussian distribution.

3.2. Effect of Non-Gaussian Measurement Noises on the Performance of Existing Robust State Estimators

State estimation aims to determine the most likely states in the system based on the measurements. One way to achieve this is through maximum likelihood estimation (MLE). Assuming that all measurement noises have known probability distributions, all measurements’ joint probability density function can be written with these unknown parameters, called the likelihood function, which will peak when the unknown parameter is chosen closest to its actual value. When the measurement noise follows the Gaussian distribution, the deduced optimal solution method is the WLS [26].

As shown in Figure 2, WLS is no longer the optimal solution when the measurement noise does not follow the Gaussian distribution. This is because, compared to the Gaussian distribution, the tail convergence rate of the bimodal Gaussian distribution and the Laplace distribution are faster than that of the Gaussian distribution tail convergence rate (so they are also known as a heavy tail distributions). Therefore, they have higher outlier probability than the Gaussian distribution, i.e., their absolute value in the white part is greater than that of the Gaussian distribution. If WLS is still used in the existence of non-Gaussian distributed noises, outliers will be assigned with large weights, resulting in reduced estimation accuracy.

In discussing the effect of non-Gaussian measurement noise on the performance of existing robust state estimators, this subsection only discusses M-estimators since the generalized M-estimator is still developed based on the M-estimators. As shown in (5), the value of the target function of LAV is proportional to the absolute value of the deviation value, so when the deviation value is large when the outlier appears, the weight is much less than the WLS, so it has high robustness for the outliers. Although there may be significant proportion of non-Gaussian measurement noises, the Gaussian measurement noise still constitutes the main part for real systems. Therefore, the LAV estimator might not achieve the highest estimation accuracy.

The Huber estimator combines the advantages of both WLS and LAV according to (4). The size of the parameter β determines Huber’s emphasis on WLS and LAV: when β→∞, Huber estimator is equivalent to WLS; when β = 0, it is equal to LAV. The parameter β is generally set to be the fixed value of 1.5. Therefore, the existing Huber estimators and generalized M-estimators based on fixed β values may have insufficient robustness or accuracy in the existence of the non-Gaussian distribution measurement noise with different proportions and parameters. Therefore, it is urgent to propose a generalized M-estimator (or Huber estimator) that optimizes the value according to the measurement noises.

3.3. Research Purpose of Generalized M State Estimation of Optimization Parameters Proposed

As shown in the above sub-sections, non-Gaussian distributed measurement noises are more likely to cause outliers than the Gaussian distributed measurement noises. Therefore, different levels of robustness of the estimator are required to ensure the best estimation accuracy under different types and proportions of non-Gaussian distributed measurement noises. Traditional generalized M-estimation based on Huber estimation with fixed parameter β can ensure high robustness against bad data and high estimation accuracy under Gaussian distributed noise. Its estimation accuracy is decreased if the measurement noises are non-Gaussian distributed. So, the purpose of this study is to propose a method that can adaptively change the value of β so that the generalized M-estimator can have high estimation accuracy under different non-Gaussian distributed measurement noises while ensuring good robustness against bad data.

4. Generalized M State Estimator of Optimized Parameters Based on Sampling

To overcome the drawback of the traditional generalized M-estimator, a novel generalized M state estimator of optimized parameters based on sampling is proposed in this section, where the random sampling method is described in Section 4.1, and the algorithm process of sampling-based optimization parameter is introduced in Section 4.2.

4.1. Optimized Parameter Selection Method Based on the Random Sampling Method

In Section 3.2, it is explained that non-Gaussian noise has a higher outlier probability than Gaussian noise. Hence, the value of β needs to be reduced to increase the estimator’s robustness. Because the optimal value is affected by many uncertain factors, such as the configuration of measurement noise, the probability distribution of measurement noise, and the probability of bad data occurrence, it is almost impossible to obtain the exact optimal value of β. Therefore, this paper proposes the sagging optimization selection method to find a value of β close to its optimal value according to the probability ratio of the outliers:

$$\beta =1.5-a\frac{{\displaystyle \sum _{i=1}^N{P}_{Truthi}}(x\ge 1.5)}{{\displaystyle \sum _{i=1}^N{P}_{Gaussi}}(x\ge 1.5)}.$$

(18)

Define $n_1={\displaystyle \sum _{i=1}^N{P}_{Truthi}}(x\ge 1.5)$, which is the number of times when the actual standardization error ((measurement value-mean)/standard deviation) of all the measurements in the system is greater than 1.5, is the quantity to be determined. Define $n_2={\displaystyle \sum _{i=1}^N{P}_{Gaussi}}(x\ge 1.5)$, which is the number of times when all the measurements obey the Gaussian distribution, is the known quantity. The variable a is the droop coefficient and the quantity to be determined. This paper takes the following steps to determine n₁ and a:

Step 1: for a system with m measurements, randomly select N measurements;

Step 2: For the N measurements, record L groups of data from their corresponding devices (L is large enough), and record the number of normalized errors greater than 1.5, n₁.

Step 3: Calculate the droop coefficient a by (19):

$$a=\left(1.5-{\beta }_{\min }\right)\frac{n_2}{n_{1,\max }}$$

(19)

where β_min is the set minimum value. In order to ensure the convergence ability of the estimator, this paper takes β_min = 0.1, which is a conservative value. n₁,max is the number of all measurements whose normalized errors are greater than 1.5 that are calculated from the probability density function at the highest non-Gaussian degree.

If most of the non-Gaussian distribution measurements in the system obey a bimodal Gaussian distribution, then the formula for the optimization values can be approximately reduced to:

$$\beta =1.5-Akh$$

(20)

where A is the droop coefficient to be sought, and k is the average bimodal Gaussian error coefficient, which is defined in Equation (14), as the non-Gaussian measurement ratio. This paper takes the following methods to determine A, k, and h, as shown in Figure 3.

4.2. The Proposed Generalized M State Estimation Algorithm of Optimized Parameters Based on Sampling

Based on the optimization parameter selection method proposed in the previous section, this paper proposes an improved generalized M state estimation algorithm. The steps are shown in Figure 4.

5. Simulated Examples

To verify the effectiveness of the proposed method, the IEEE 14 and 118 bus test systems are selected to simulate the performance generalized M-estimator with different parameters. All simulations were conducted in MATLAB, using an Intel Core i7-9750CPU (@2.6 Hz) 16 GB memory computer. Section 5.1 describes the impact of bad data and non-Gaussian measurement noise on generalized M-estimators with different values of β by conducting exhausted Monte Carlo simulations in the IEEE 14 bus test system. Section 5.2 demonstrates the effectiveness of the proposed generalized M-state estimator of the optimized parameter in a large system, the IEEE 118 bus test system.

5.1. Effect of Bad Data and Measurement Noise on the Performance of the Generalized M-Estimator

To ensure the generality of the results, all of the results shown in this section are obtained by averaging the results over the 500 Monte Carlo simulations. To ensure the validity of the generalized M-state estimation method, a set of redundant conventional measurements and PMU measurements were selected in IEEE 14 bus test system as shown in Figure 5.

For Gaussian errors, the standard deviation is set to a percentage of the measurement according to the type of measurement as follows:

SCADA: P, Q measured value (0.02), voltage (0.002)

PMU: Voltage amplitude (0.002), Voltage phase angle (0.01)

The performance of the generalized M-estimator can be adjusted by varying values of β. In this simulation test, the least value of β is 10⁻². With this configuration, the generalized M-estimator can be considered to be an LAV estimator. The largest value of β is set to be 10⁴. The corresponding estimator can be considered to be a WLS estimator. In the range of 10⁻² to 10⁴, β increases at the fixed ratio of 10^1/30, all estimators with different values of β are tested in the IEEE 14 bus test system. The root mean square error (RMSE) is used as the performance index of the estimators.

Case 1: Existence of bad data with fixed percentage errors. We compared the robustness of the state estimators by introducing fixed percentage errors to the voltage magnitude measurement at bus 1 and obtaining the RMSE of the estimators when the noise probability density function follows the standard Gaussian distribution. The voltage magnitude measurement at bus 1 is changed to 0, 0.2, 0.4, 0.6, 0.8 of its true value (100%, 80%, 60%, 40%, 20% errors, respectively), and tested separately. The simulation results are shown in Figure 6.

The following conclusions can be obtained:

(1) After introducing fixed errors of different magnitudes, the trend of RMSE of each curve is roughly the same: the general trend is that RMSE increases as a larger β is selected. It will not change after a certain value;

(2) There is an optimal value of β for each fixed percentage error, and the value decreases with the percentage of the error parameter, as summarized in

Table 2. Optimal β values of each estimator with different fixed error parameters.

Percentage of Fixed-Error Parameters	100%	80%	60%	40%	20%
The Optimal β Value	0.3415	0.5412	0.6813	0.8577	2.3263

.

Case 2: Existence of bad data with random percentage of errors. The robustness of generalized M state estimators with different values of β is compared by five individual simulations with three random bad data points in the measurements and observing their performances. The test results are shown in Figure 7.

The following conclusions can be obtained:

(1) When there are no bad data in the measurements, the trend of RMSE with the change of β is decreasing in general; the value of β which achieves the lowest RMSE is 1.5, which is the exact value of β for normal generalized M-estimator;

(2) After the introduction of three random errors, the trend of the RMSE of each simulation with the change of β is roughly the same: The general trend is increasing. It will not change after a certain value;

(3) There is an optimal β whose corresponding estimator achieves the lowest RMSE for each simulation, which varies between 10⁻¹ and 10¹ with different bad data.

Case 3: Existence of measurements with bimodal Gaussian distribution noise. When the bimodal interval size k (as introduced in Section 3.1) is between 0~5, the trend of the RMSE with estimators with different values of β is shown in Figure 8.

The following conclusions can be obtained:

(1) When k ≤ 1, RMSE decreases with the increase of β, and when k > 1, RMSE decreases with the increase of β;

(2) The optimal β whose corresponding estimator achieves the lowest RMSE decreases with the increase of k;

(3) In the presence of non-Gaussian measurement noise, the value range of the optimal β is between 10⁻¹ and 10¹, which can roughly cover the value range of the optimal β in the presence of bad data;

(4) Since the occurrence of bad data is low and the value range of the optimal β in the bad data case can be roughly covered by that in the non-Gaussian measurement noise case, the optimal β can be determined by only considering the influence of non-Gaussian measurement noise.

The optimal βs obtained by simulation in the presence of measurement noises in bimodal Gaussian distribution with different values of k are summarized in

Table 3. Optimal β values in the presence of different bi-modal Gaussian distributed measurement noises.

The Bimodal Interval Size (k)	0	1	2	3	4	5
The Optimal β	1.5000	0.4642	0.2712	0.1849	0.1324	0.1000

and plotted in Figure 9. The optimized βs calculated according to the proposed optimized parameter selection method (20) are also plotted in Figure 9. It can be seen that with the increase of k, the calculated optimized β decreases, and are equal to the optimal β at k = 0 and k ≈ 5. Although the optimized β is larger than the optimal β for 0 < k < 5, the decreasing trend guarantees significantly better estimation accuracy compared to the estimator with a fixed c value while k increases.

5.2. Simulation Examples of the IEEE118 Bus Test System

To verify the proposed generalized M-estimation method of optimized parameters based on sampling, simulations are conducted in the IEEE 118 bus test system with the network diagram and measurement configuration shown in Figure 10.

Among the total 186 measurements in the system, 28 of them have bimodal Gaussian distributed noises with the bimodal interval size (k) randomly set between 0~5. According to the proposed optimized parameter selection method, the optimized β, β opd, is determined as follows:

Step 1: Randomly select 80 measurements;

Step 2: For each measurement selected, test the measurement noises and record 10,000 groups of data; 12 non-Gaussian measurement noises are counted according to the statistics, and the proportion of non-Gaussian measurement, h, is 15%;

Step 3: Calculate the sagging coefficient A:

$$A=\frac{1.5-{\beta }_{\min }}{k_{\max }{h}_{\max }}=\frac{1.5-0.1}{5\times 0.5}=0.56$$

where k_max is the maximum bimodal Gaussian error coefficient possible, and h_max is the maximum proportion of the bimodal Gaussian distribution measurement noise possible.

Step 4: Find the average k according to the statistic bar chart of all bimodal Gaussian distributed measurements noises. To test the proposed method in different scenarios, 5 different simulations with average k ≈ 1, 2, 3, 4, 5 are conducted with the statistical bar charts shown in Figure 11.

To verify the effectiveness of the proposed generalized M-estimator of optimized parameters, the RMSEs of the estimators obtained from Monte Carlo simulations for different ks and against β are shown in Figure 12. The optimized β, β_opd, can be calculated according to (20) with the obtained A and h. The optimal β, β_opl, and the RMSEs at β_opl, β_opd, β = 1.5 for different ks can be obtained from Figure 12. The results are summarized in

Table 4. Average standard residuals of different estimators in IEEE118 systems at different non-Gaussian noise ratios.

k	5	4	3	2	1
βopd	0.1005	0.3804	0.6603	0.9402	1.2201
βopl	0.1001	0.1467	0.2326	0.3687	0.9261
RMSE at βopd	1.55 × 10−3	1.49 × 10−3	1.45 × 10−3	1.17 × 10−3	1.15 × 10−3
RMSE at βopl	1.52 × 10−3	1.43 × 10−3	1.26 × 10−3	1.06 × 10−3	1.11 × 10−3
RMSE at β = 1.5	1.82 × 10−3	1.73 × 10−3	1.63 × 10−3	1.32 × 10−3	1.18 × 10−3
((RMSE at β = 1.5) − (RMSE at βopd))/RMSE at β = 1.5 × 100%	14.84%	13.87%	11.04%	11.36%	2.54%

.

The analysis performed in Figure 12 and Table 4 leads to the following conclusions:

(1) β_opd is larger than β_opl, but the RMSE is not much different from the actual RMSE;

(2) The RMSE at β_opd is significantly lower than the RMSE at β = 1.5.

The conclusions above demonstrate the effectiveness of the proposed method. The underlying reason for these results is that the proposed method reduces the value of β to increase the estimator’s robustness when the level of non-Gaussian measurement noise increases. The linear drooping characteristic of β gives the proposed estimator significantly better estimator accuracy than the traditional generalized M-estimator, but slightly lower estimator accuracy than the estimator with β_opl, which is acceptable considering the difficulty to obtain β_opl.

5.3. Simulation Examples of the Polish 2736 Bus System

To verify that which is proposed, simulations are conducted in the Polish 2736 bus test system [27].

Among the total 3269 measurements in the system, 280 of them have bimodal Gaussian distributed noises with the bimodal interval size (k) randomly set between 0~5. According to the proposed optimized parameter selection method, the optimized β, β opd, is determined as follows:

Step 1: Randomly select 800 measurements;

Step 2: For each measurement selected, test the measurement noises and record 10,000 groups of data; 120 non-Gaussian measurement noises are counted according to the statistics, and the proportion of non-Gaussian measurement, h is 15%;

Step 3: Calculate the sagging coefficient A;

Step 4: Find the average k according to the statistic bar chart of all bimodal Gaussian distributed measurements noises. To verify the effectiveness of the proposed generalized M-estimator of optimized parameters, the RMSEs of the estimators obtained from Monte Carlo simulations for different ks and against β are shown in

Table 5. Average standard residuals of different estimators in Polish 2736 systems at different non-Gaussian noise ratios.

k	5	4	3	2	1
βopd	0.6751	0.8330	0.9876	1.167	1.340
βopl	0.1079	0.1467	0.2326	0.3687	0.9261
RMSE at βopd	1.69 × 10−3	1.54 × 10−3	1.39 × 10−3	1.22 × 10−3	1.10 × 10−3
RMSE at βopl	1.61 × 10−3	1.48 × 10−3	1.34 × 10−3	1.16 × 10−3	1.05 × 10−3
RMSE at β = 1.5	2.03 × 10−3	1.96 × 10−3	1.78 × 10−3	1.57 × 10−3	1.43 × 10−3
((RMSE at β = 1.5) − (RMSE at βopd))/RMSE at β = 1.5 × 100%	16.74%	21.43%	21.91%	22.29%	23.08%

.

Analysis performed in Table 5 leads to the following conclusions:

(1) The β value calculated by sagging coefficient is slightly different from the actual β value, but its RMSE is not much different from the actual RMSE;

(2) The RMSE corresponding to the calculated β value is significantly higher than that when β = 1.5, and the RMSE is significantly higher when the k value is larger.

These conclusions are almost the same with those in Section 5.2, demonstrating the adaptability of the proposed even in very large networks.

6. Conclusions and Future Works

This paper presents a generalized M-estimation method for power system state estimation with sampling-based optimization parameters. The sampling statistics method can set the parameter according to the different non-Gaussian measurement noise parameters and the proportion to adjust its robustness, which increases its overall estimation accuracy. The generalized M-estimation defines standardized residuals based on the median, thus being more robust than ordinary M-estimation methods. Exhaustive simulation results in the IEEE14 bus test system demonstrate the impacts of non-Gaussian measurement noise and bad data on estimation accuracy, which are summarized as follows:

(1) Small β can increase the robustness of the estimator but decrease the estimation accuracy in the normal operation case where no bad data occur, and vice versa. There is an optimal β that can achieve the highest overall estimation accuracy;

(2) The optimal β decreases with the increase of the magnitude of bad data and the non-Gaussian degree of measurement noise;

(3) The selection of the optimal β value mainly depends on the non-Gaussian degree and proportion of the measured noise.

Simulations in the IEEE 118 bus test system and Polish 2736 bus system verify that the proposed method can improve the estimation accuracy by up to 23% compared to the traditional generalized M state estimator in the presence of different degrees of non-Gaussian measurement noise.

Although this paper has reached the expected goal, further research can be carried out by considering more newly developed estimators, adopting more complex grid models, and considering more complex situations.