Efficient Identification of Multiple Parameter Errors in Power Grids by Mixed-Effects Models and Generalized Least Squares

Accurate knowledge of static parameters forms the basis for nearly all applications in the energy management system. This paper proposes an efficient method for the simultaneous identification and correction of multiple network parameter errors based on a linear mixed-effects (LME) model and the generalized least squares (GLS) method. An LME model for parameter error identification is formulated using the residual equations of multiple-snapshot state estimation with equality constraints. The parameter errors are considered as the fixed effects and the measurement errors are considered as the random effects. Then, using the measurement error variances estimated by the LME model, the GLS method is used to estimate the parameter errors along with a hypothesis testing to infer whether each parameter error is zero. The semi-supervised adversarial autoencoder is used for bad data detection in the presence of erroneous parameters and limited labels such that only measurement snapshots that are free of any bad data are used. The proposed methodology is efficient in that the LME model is only used to estimate the variances of the measurement errors using a small number of measurement snapshots, therefore the huge computation burden needed to solve a large-scale LME model is avoided. In addition, the GLS only involves inversion of low-dimension matrices, which is very efficient even a large number of measurement snapshots are used. Thorough tests of the proposed methodology on a large number of scenarios are provided to show the effectiveness of the proposed methodology with promising results.


INTRODUCTION
Correct static parameters of transmission lines and transformers are essential for many applications in the energy management systems of all electric utilities, such as state estimation (SE), power flow, security assessment, etc. For example, it has been reported by the federal energy regulatory commission that finding a good solution for AC optimal power flow (OPF) could potentially save tens of billions of dollars annually (Cain et al., 2012). However, the OPF results may be hugely distorted and may even be infeasible due to potential errors in network parameters.
Traditional power system SE (Monticelli, 1999;Abur and Gómez-Expósito, 2004;Zhao et al., 2019) usually assumes that all of the parameters stored in databases are correct. However, it has been found that parameter deviation in power systems can be as large as 30% (Kusic and Garrison, 2004). Therefore, it is very important to identify and correct erroneous parameters in power networks.
The parameter error identification (PEI) and estimation problem has been studied for decades and several approaches have been proposed in the past (Zarco and Exposito, 2000). Relevant research in this field falls roughly into two categories according to different assumptions: constant parameters or time-varying parameters. Assuming constant model parameters during a short period, parameters can be periodically estimated and updated by operators. The contribution of this paper falls in this category. In (Debs, 1974), a recursive filtering-type algorithm is derived to correct inaccurate parameters by processing collected online data in an off-line mode. Parameters are augmented into the state vector, and existing parameters in the database are used as a-prior information. This method may suffer from modifying correct parameters to false values due to measurement noise. In (Liu et al., 1992;Liu and Lim, 1995), the authors propose to estimate parameters from measurement residuals of SE, which can be interpreted as a linear model linking measurement residuals to an unknown parameter error in the presence of noise. A suspicious parameter set is built based on normalized residuals to reduce computation cost. In (Castillo et al., 2011), a three-stage offline approach is presented to detect, identify, and correct series and shunt branch parameter errors. Suspicious parameters are flagged through a heuristic identification index, and are then estimated via an augmented state estimator, and are finally validated via a weighted least squares estimator.
Recently, a largest normalized Lagrange multiplier (NLM) method is proposed in (Zhu and Abur, 2006;Lin and Abur, 2018), which identifies an erroneous parameter to be the one with the largest NLM. The method is then extended to use multiple measurement snapshots (Zhang and Abur, 2013;Lin and Abur, 2017). The method has shown good performance for single or multiple non-interacting parameter errors, and is also able to distinguish erroneous parameters from bad data. However, the NLM method identifies parameter errors in a sequential manner, and thus may encounter the smearing and masking effect when multiple-interacting parameter errors exist. In (Zhao et al., 2018), the projection statistics is used for selection of suspicious parameters associated with leverage points, which is shown to be insensitive to the smearing and masking effects. In (Liang et al., 2021), a linear mixed-effects (LME)-based method is proposed for simultaneous identification of multiple network parameter errors. Suspicious parameter set is constructed by solving an LME model and hypothesis testing. An extended state augmentation (ESA) method is further proposed so that modifying correct parameters to false values can be avoided as much as possible. The method is shown to be promising for simultaneous identification of multiple network parameters, assuming all bad measurements have been removed. However, solving an LME model is challenging and expensive, especially when a large number of measurement snapshots are used. Due to this reason, only several typical scenarios are tested and a thorough testing on large number of scenarios is missing.
On the other side, assuming model parameters are timevarying, parameter tracking is also getting more and more attention. In (Slutsker et al., 1996;Bian et al., 2011), Kalman filter-based recursive parameter estimation methods are presented. The methods treat parameters as Gaussian random variables and can track parameters as they fluctuate due to changes in load and ambient conditions. In (Williams et al., 2016), the residual sensitivity analysis (RSA) method used in (Liu et al., 1992;Liu and Lim, 1995) is adopted for offline parameter tracking and the discovery of non-diurnal and nonseasonal changes of parameters in unbalanced distribution systems. The method leverages the increased deployment of distribution level measurement devices for a three-phase state estimator to estimate changes in impedance parameters over time. In (Ren et al., 2017), the authors propose to estimate the parameters of untransposed overhead transmission lines using synchronized measurements. SE and Kalman filter-based parameter tracking process are iteratively carried out to reduce the uncertainties that exist in the estimation function. All above references base their method on power system SE theory in a network-wide manner. Recently, non-SE type methods are also proposed for parameter tracking using local measurements, see (Wang et al., 2015;Dobakhshari et al., 2020) and references therein.
This paper proposes an efficient method for the simultaneous identification of multiple network parameter errors based on a LME model and the generalized least squares (GLS) method. The contributions are summarized as follows: 1) An LME model for PEI is formulated using the residual equations of multiple-snapshot SE with equality constraints. The parameter errors are considered as the fixed effects and the measurement errors are considered as the random effects. Using the measurement error variances (MEVs) estimated by solving the LME model, the GLS method is adopted to estimate the parameter errors along with a hypothesis testing to infer whether each parameter error is zero. The proposed methodology is efficient in that the LME model is only used to estimate the variances of the measurement errors using a small number of measurement snapshots, therefore the huge computation burden needed to solve a large-scale LME model is avoided. In addition, the GLS only involves inversion of low-dimension matrices, which is very efficient even a large number of measurement snapshots are used. 2) The semi-supervised adversarial autoencoder (semi-AAE) is used to detect existence of any bad data in each measurement snapshot such that only measurement snapshots that are free of any bad data are used. Thanks to the generalization ability of neutral networks, the method can achieve high bad data detection (BDD) accuracy even in the presence of erroneous parameters and limited labels. 3) Thorough tests on a large number of scenarios are provided to show the effectiveness of the proposed methodology.
The remainder of this paper is organized as follows. Power System State Estimation Section briefly introduces the power system SE theory. The LME Method for Estimation of MEVs Section presents the LME methodology. The Proposed LME-GLS Methodology Section presents the proposed LME-GLS methodology. Simulation Results Section describes the simulation results. Conclusions are drawn in Conclusion Section.

POWER SYSTEM STATE ESTIMATION
In this section, the basic theory of equality constrained power system state estimation is introduced, which forms the basis for the following sections on parameter error identification.
Consider a power network with the measurement model given by: where x is N s ×1 vector of the state variables consisting of voltage amplitudes and phase angles of all buses expect phase angle of the reference bus, N s = 2N b -1; z is the N m ×1 measurement vector, including real and reactive power injection measurements, real and reactive power flow measurements, and voltage magnitude measurements; h(·) is the N m ×1 measurement function vector; e is the N m ×1 measurement error vector, e ∼ N(0, R); R is the N m ×N m covariance matrix of the measurement error vector e; N b is the number of buses; N m is the number of measurements. We consider the following equality constrained weighted least squares (WLS) SE formulation: where r is the N m ×1 measurement residual vector; c(·) is the zeroinjection measurement function vector. The Lagrange function can be constructed as: where μ, λ are Lagrange multiplier vectors of the equality constraints.
Taking the derivative of the Lagrange function, we get the firstorder optimal condition: where H(·), C(·) are Jacobian matrices of the measurements and zero-injection measurements with respect to (w.r.t.) state variables evaluated at a suitable operating point.
Using a Taylor expansion at an initial operating point x 0 , the nonlinear equation can be linearized as follows: Starting from an initial state x 0 , the nonlinear equation can be solved iteratively. If we let the inverse of the coefficient matrix at the kth iteration be: (6) then at the last iteration before the iteration converges, starting from the true operating point x (x 0 = x), we get the following residual equation which relates the measurement residuals and measurement errors together: where matrix E 5 is evaluated at the estimated statex. If parameter errors exist, the linearization of residual at the true operating point x becomes r z-h(x)-H(x)Δx-H p (x)p e e-HΔx-H p p e , where p e is the N p ×1 parameter error vector and N p is the number of parameters. Then the residual equation becomes: where H p is the N m ×N p Jacobian matrix of the measurements w.r.t. parameters evaluated atx.
If there are no zero-injection measurements, the residual Eq. 8 reduces to the following widely known format: where S is the residual sensitivity matrix evaluated atx. As (8) and (9) have similar format, we will put an emphasis on Eq. 9 in the following. However, all the analysis are applicable if S is replaced by E 5 .
The residual equation provides an accurate linear approximation of the nonlinear measurement equations at the estimated operation point. Considering measurement error vector follows independent Gaussian distribution, i.e. e i ∼ N(0, R), then we have r i ∼ N(-S i H p,i p e , S i RS i T ) and it seems that p e can be estimated by GLS. However, S i is singular whose rank is equal to N m -N s and thus GLS cannot be directly used.

THE LME METHOD FOR ESTIMATION OF MEVS
Based on the SE theory, this section presents the LME method (Liang et al., 2021), which is then improved in the next section to improve the computation efficiency. This LME method has shown to be very effective for simultaneous identification and correction of multiple network parameters. However, solving an LME model is challenging and expensive, especially when a large number of measurement snapshots are used. Theoretically, if MEVs can be accurately known, the LME method reduces to a much easier GLS method. However, directly solving the GLS model using empirical values of MEVs may lead to failing to identify the erroneous but less sensitive parameters. Fortunately, it is found that estimation of MEVs does not rely on a large number of measurement snapshots. In fact, a small number of measurement snapshots are enough to get acceptable estimates for MEVs, which greatly reduces the computation burden. Therefore, in this section we propose to use the LME model not for PEI, but only to estimate the MEVs using a small number of measurement snapshots, therefore the huge computation burden needed to solve a large-scale LME model is avoided. We briefly present the LME method in the following.
Formulation of the LME Model LME models extend linear models by incorporating random effects to account for correlation among observations. An LME model contains a population of individuals, and each individual has a group of observations. Fixed effects are associated with the entire population and random effects are associated with each individual. For the ith individual with n i observations, an LME model describing a response vector y i is as follows (Pinheiro and Bates, 2000): where y i is the n i ×1 response vector of the ith group; β is the p × 1 fixed effect vector; b i is the q × 1 random effect vector of the ith are n i ×p and n i ×q full-column-rank design matrices associated with the fixed effects and random effects, respectively. The response vectors are correlated within a group but independent among different groups. The equations of N groups can be stacked into one: For ease of simplicity, we can assume that N measurement snapshots share an equal number of measurements N m , then the ith residual equation is: By adding a virtual error term ε i , we can transform the residual Eq. 13 into an LME model: where the fixed effects β and random effects b i of the general LME model (10) are replaced by the deterministic parameter error vector p e and stochastic measurement error vector e i , respectively; the design matrices of the fixed effects and random effects can be constructed by Note that it is acceptable to introduce a virtual error term ε i , obeying a Gaussian distribution with an unknown but very small standard deviation σ that will automatically be estimated and does not affect the performance of the method. We set e i ∼ N (0, R) = N(0, σ 2 Ψ), where σ 2 is an extracted scale parameter equal to variance of the error term ε i and Ψ = R/σ 2 is the remaining term. At this time r i ∼ N(X i p e , σ 2 (I + Z i ΨZ i T )), and the covariance matrix σ 2 (I + Z i ΨZ i T ) is now nonsingular and can be viewed as a damped version of the original covariance matrix Z i RZ i T . However, the statistical properties of the measurement errors are still unknown as R −1 is usually casually set by experience and does not accurately reflect the realistic statistical properties. This motivates the usage of an LME model rather than directly usage of the GLS method.
Using N measurement snapshots, the LME model can be stacked as: Suppose all the measurements in each snapshot are divided into three groups (or random-effects terms) by measurement types: real power measurements, reactive power measurements, and voltage magnitude measurements. We Then, by proper splitting and reordering, (15) can be written as: where the three groups of measurements have standard deviations σ P , σ Q , and σ V , respectively. It should be noted that the measurements should be divided based on the types and precision levels of the instruments and are not necessarily divided into these three groups.

Maximum Likelihood Estimation of MEVs
Let the unknown Ψ be parameterized by a vector θ, i.e., Ψ becomes Ψ(θ). As Ψ denotes a covariance matrix and the measurements are assumed to be independent, Ψ should be a diagonal matrix, and θ should have a dimension of N m . Furthermore, assuming that the measurements in the same group share the same standard deviation, Ψ should be an isotropic matrix, and θ only has a dimension of three. By applying the Bayesian formulas, the likelihood of the observed residual r, given parameters p e , θ, and σ 2 , is (Pinheiro and Bates, 2000;Bates et al., 2015): where Frontiers in Energy Research | www.frontiersin.org February 2022 | Volume 10 | Article 840736 and |A| denotes the determinant of matrix A. Solving an LME model includes the following steps (refer to (Pinheiro and Bates, 2000;Bates et al., 2015) for more details): Step 1: Maximize P(r|θ,p e ,σ 2 ) with respect to p e and σ 2 for a given θ. The solutions p e (θ) and σ 2 (θ) are functions of θ.
It should be noted that the LME model can be solved only if matrix X = -SH p has full column rank. If the measurement redundancy is low and X does not have full column rank, then an identifiability analysis should be performed to find the unidentifiable parameters.

THE PROPOSED LME-GLS METHODOLOGY
This section presents the proposed LME-GLS methodology as a whole. The key lies in that the MEVs and parameter errors are no longer estimated simultaneously but in a sequential manner. Specifically, using the MEVs estimated by solving the LME model, the GLS method is adopted to estimate the parameter errors along with a hypothesis testing to infer whether each parameter error is zero. In addition, the semi-AAE is used to detect existence of any bad data in each measurement snapshot such that only measurement snapshots that are free of any bad data are used. The total flow chart of the proposed LME-GLS methodology is shown in Figure 1 and the details are explained in the following.
Step 1: Bad Data Detection by Semi-AAE Due to various disturbances in realistic power systems, it is necessary to detect bad data so as to prepare an accurate measurement dataset for subsequent PEI. In this paper, the BDD problem is viewed as a binary classification problem and solved by the semi-AAE (Makhzani et al., 2015). The semi-AAE model is designed and trained to accept each measurement snapshot as a data sample and output a detection indicator β. β = 1 indicates existence of bad data and β = 0 otherwise. The semi-supervised property makes it suitable for scenarios when only a small portion of measurement snapshots are labeled.
The AAE integrates the generative adversarial network (GAN) into the autoencoder (AE) framework. An encoder learns to convert the input data vector to a latent vector, while a decoder learns to reconstruct the original data vector from the latent vector. The aggregated posterior distribution of the latent data is defined as where X is the N m ×M input data; H is the q×M latent data; q(H|X) is the encoding distribution of H given X; p d (X) as the true distribution of X; M is the number of samples; q is the dimension of each latent vector. AAE introduces a GAN model to regularize the aggregated posterior distribution q(H) of the latent representation to an arbitrary prior distribution p(H), such as the standard Gaussian distribution. The encoder acts as the generator of the GAN model. It tries to learn an aggregated posterior distribution q(H) to fool the discriminator of the GAN and make it falsely believe a latent vector sample comes from q(H) but is actually from p(H). The discriminator D Gauss , on the other side, tries to classify the input samples correctly. Furthermore, to make good use of the label information of the input samples for classification use, the latent vector is augmented with an additional category variable vector using one-hot format along with an additional discriminator D Cat . Therefore, H consists of two parts: category latent data H C and continuous latent data H G , which are forced to approximate a binary categorical distribution Cat(2) and a standard Gaussian distribution N(0, I), respectively. In this situation, the encoder q(H C , H G |X) works as the generator of both GANs so that it can simultaneously predict fake category data H C and fake continuous data H G . The real data of the two discriminators, i.e. H C ' and H G ', can be randomly sampled from Cat (2) and N(0, I).
The AAE model can be trained in three stages: reconstruction stage, regularization stage and semi-supervised classification stage: 1) In the reconstruction stage, the AE updates its encoder q(H C , H G |X) and decoder to minimize the reconstruction loss where the input is the unlabeled data. The encoder maps the input data X to the representation H of the hidden layer: where h i is the ith column vector of H; x i is the ith column vector of X; W is the q×N m weight matrix; b is the q × 1 offset vector; s(·) is a nonlinear activation function. The model parameters of the encoder are expressed as θ = {W,b}. The decoder maps the latent data H to the reconstructed data X′ through the mapping function f′: where x i ′ is the ith column vector of X′; W′ is the N m ×q weight matrix; b′ is the N m ×1 offset vector. The model parameters of the decoder are expressed as θ' = {W′,b′}. The reconstruction error Loss R is minimized to obtain the optimized parameters: 2) In the regularization stage, the discriminator parameters of the GAN models are updated first, and then the generator parameters are updated. The discriminator should distinguish real samples from fake samples as much as possible, while the generator tries to generate fake data similar to the real samples to deceive the discriminator, resulting in a two-player minmax game formulated as where E is the expectations under a distribution; p(H C ), p(H G ) are the prior distributions of H C , H G ; G(·) is a function mapping samples from the prior distributions to the latent data space; D Cat (·) is a function which returns the probability that the input sample comes from the real categorical distribution Cat(2) (positive samples), rather than from our generative model (negative samples); D Gauss (·) is a function which returns the probability that the input sample comes from the real standard Gaussian distribution N(0, I) (positive samples), rather than from our generative model (negative samples).
3) Finally, in the semi-supervised classification stage, the encoder q(H C |X) is updated only by using the labeled data to minimize the cross-entropy loss, which is expressed as where q(H C ) is an aggregated posterior distribution of H C .
Step 2: MEV Estimation by Solving the LME Model Using the clean measurement dataset without outliers, the LME method presented in The LME Method for Estimation of MEVs Section is then used to estimate MEVs.
Step 3: Parameter Error Estimation by the GLS Method Using the clean measurement dataset without outliers, the LME method is then used to estimate MEVs. Once the statistical properties of the measurement errors are estimated, then we have r i ∼ (X i p e ,σ 2 (I + Z iΨ Z T i )). As follows from the Gauss-Markov theorem, p e can be estimated directly by the GLS method: As can be seen from (26), the GLS only involves inversion of low-dimension matrices, which is very efficient even a large number of measurement snapshots are used.

Step 4: Suspicious Parameter Selection by Hypothesis Testing
Although p e can be estimated through GLS, the estimatedp e is not suitable for parameter correction asp e will have many nonzero elements, even elements corresponding to correct parameters may be relatively large due to measurement noise, and the optimal solution will show some amount of "overfitting". Therefore,p e is only used for hypothesis testing to test whether each parameter error is zero or not.
The covariance matrix ofp e can be calculated by The ith component of normalizedp e approximately follows a t-distribution with N T -N p degrees of freedom (Pinheiro and Bates, 2000): where γ i is the t-statistic of the ith parameter error; C is the N p ×N p covariance matrix ofp e ; C ii is the ith diagonal element of C; N T is the number of total measurements, N T = NN m .
The hypothesis testing identification (HTI) is then used to decide if a parameter is erroneous. The null hypothesis and alternative hypothesis are chosen as follows: The statistical significance of the ith parameter error can be tested by its two-sided p-value: where t α/2,d is the (1-α/2)th quantile of the t-distribution with d degrees of freedom, which indicates the probability of obtaining a value γ j more extreme than the critical value. p-values below the significance level provide strong evidence for rejecting the null hypothesis.
Step 5: Parameter Correction by the ESA Method The HTI may make two types of errors known as type I errors and type II errors. A type I error is the rejection of H 0 when it is indeed true. It is also referred to as a "false alarm" because correct parameters are misidentified as erroneous. To this end, the ESA method (Liang et al., 2021) is further used to move misidentified but correct parameters from the suspicious parameter set back to the correct parameter set while correcting the identified erroneous parameters.

Comparison With Existing Methods
The proposed method seems similar to the traditional RSA method (Liu et al., 1992;Liu and Lim, 1995). However, there are two differences between the LME-GLS method and the RSA method.
1) The RSA method is used to estimate parameter errors whereas the LME-GLS method is used to identify parameter errors. Due to different measurement configurations, measurement noise levels, and limited number of measurement snapshots used, it is possible to modify the existing values in the database, even when they are correct. On the contrary, in the LME-GLS method, only a limited number of parameters with strong evidence by data that they are undoubtedly erroneous are selected for estimation, which avoids modifying existing values to any false values as much as possible.
2) The RSA method formulates a linear model linking the residuals and the unknown parameter errors and assumes the measurement noise has a known statistic property. However, the proposed LME-GLS method formulates the residual equation as an LME model in which parameter errors and measurement noise are viewed as fixed effects and random effects, respectively. The parameter errors and variances in the measurement noise are simultaneously estimated from longitudinal residual data, followed by an HTI process.
The proposed method although cannot handle bad data directly, it makes use of a novel learning-based BDD method to guarantee only measurement snapshots that are free of any bad data are used. The comparison of the proposed learning-based bad data detection (BDD) method with classic BDD method, i.e. the Chi-squares test method using objective function values of state estimation, or the largest normalized residual (LNR) method, are summarized as follows: 1) Classic BDD methods are able to detect single bad data or multiple non-interacting bad data. However, they are incapable of detecting multiple interacting bad data whose residuals are masked by each other and the objective function value or their (normalized) residuals are smaller than the threshold even they exist. Note that the famously known false data injection attack is a special type of multiple interacting bad data, where hackers can circumvent the BDD module by tampering with the sensor measurements and injecting the elaborately designed false data without being detected. 2) Classic BDD methods depends on convergence of SE algorithms, which may be time-consuming or hard to converge, especially for large scale systems. On the other side, the learning-based BDD method is free of convergence problems and is very fast once it is well trained. In the meantime, it is very effective and can achieve very high accuracy, taking advantage of the strong feature extraction ability of neural networks.

SIMULATION RESULTS
In this section, thorough tests on a large number of scenarios are provided to show the effectiveness of the proposed methodology. Multiple measurement snapshots are generated by varying the loading condition between approximately 0.8 and 1.2 times their nominal values using MATPOWER (Zimmerman et al., 2011). The measurements are divided into five groups: real power flow measurements (PF), real power injection measurements (PI), reactive power flow measurements (QF), reactive power injection measurements (QI), voltage magnitude measurements (VM), and different levels of measurement noise are used for each group. The selected line parameters are added with an error value of +30%. The significance level for the p-values is set to 0.05.

BDD Results
In this section, the IEEE 14-bus system is used to test the effectiveness of the semi-AAE-based BDD method in the presence of erroneous parameters and limited labels. There are 20 series branches (transmission lines and transformers) and 40 series parameters (resistances and reactances) in this system. +30% relative error are added to parameters r 2-4 , r 2-5 , r 3-4 and r 4-5 . Full measurement configuration is assumed. Relative measurement errors are added to the true values of measurements by z i = z i,true ×(1 + e type ×σ) where z i , z i,true are measured and true value of measurement i, e type is percent relative error corresponding to measurement i's type, σ is a random variable obeying standard Gaussian distribution. Measurement errors are set as e PF = e QF = 0.5%, e PI = e QI = 1.0%, e VM = 0.2%. Monte Carlo simulation is used to generate 5000 normal measurement snapshots without bad data. Relatively large errors (10-20% of the normal values) are added to 8-12% of measurements in each measurement snapshot to generate another 5000 measurement snapshots with bad data. 80% of the whole 10,000 measurement snapshots are selected as the training dataset and the left 20% as the testing dataset. Only 10% of all snapshots are labeled β = 0 or 1, including 500 normal measurement snapshots with label β = 0 and 500 erroneous measurement snapshots with label β = 1, respectively. The performance of the semi-AAE method is evaluated by the following accuracy rate (A), which refers to the percentage that the detection results are correct among the prediction results of all samples: where TP means true positive, TN means true negative, FP means false positive, and FN means false negative. Their specific definitions are shown in Table 1.
The accuracy, detection time and training time are shown in Table 2. It can be seen that the detection accuracy on the training set and testing set reach 96.5 and 94.5%, respectively, showing that the algorithm is effective for BDD. At the same time, the training time is only 124.38 s, and the detection time is only 0.66 s, which shows that the AAE can meet the requirements of online BDD.
Different combinations of erroneous parameters are simulated to test the performance of the proposed methodology. If N e parameter errors are assumed to exist, then there will be C Ne N p different combinations of erroneous parameters, among which 100 combinations are selected for testing. In the LME step (Step 2 in Figure 1), N LME measurement snapshots are used to test the estimation accuracy of the MEVs and its impact on the whole success rate. In the GLS-HTI step (Step 3-4 in Figure 1), maximally 100 measurement snapshots (N GLS = 10, 20, . . . , 100) are used for parameter error estimation and hypothesis testing.
Simulation results show that three snapshots are enough to obtain acceptable MEVs. The estimated MEVs of the four cases when N LME = 3 are summarized in Table 3. If larger N LME is used, the estimated MEVs only slightly differ from the current values but much heavier computation burden is needed.
Then, based on the estimated MEVs in Table 3, the GLS is carried out to test the success rates of PEI. For each PEI test, the result is defined as "success" if all the erroneous parameters' p-values rank top k ones among the whole N p parameters' p-values sorted in ascending order. k = min(N max , N threshold ), where N max is a constant depending on N e (N max = 2N e in this paper), N threshold is the number of parameters whose p-values are smaller than a predefined significance level (0.05 through all simulations). In this way only N max parameters are selected as suspicious if a large number of parameters fail the hypothesis testing and thus the maximum number of suspicious parameters can be controlled.
The success rate results of the four cases are shown in Figures  2A-D respectively, where N e can be 2, three or 4. It can be seen that the success rates increase with increasing measurement snapshots. For Case 1 with the highest measurement accuracy and the largest measurement redundancy, Figure 2A shows that when the number of measurement snapshots N is 10, the minimum success rate is only about 60%. However, when the number of measurement snapshots is 20, the success rate reaches above 90%. By increasing N, the success rates can be higher than 95%, or even 100%. For Case 2 with the same redundancy and lower measurement accuracy, Figure 2B shows that the minimum success rate is only 33% when N is 10. As N increases, the success rate also gradually increases, but the final result is less than 100%. For Case 3 with high precision and low redundancy, Figure 2C shows that when N is 10, the minimum success rate is only 38%, but with the increase of N, the success rate gradually increases and can reach 100%. Finally for Case 4 with lowest measurement accuracy and redundancy, Figure 2D shows that the minimum success rate reduces to 10% when N is 10, and the highest success rate is only 82% when N increase to 100.
Comparing the results of Case 1-4 in Figure 2, it can be seen that the accuracy of measurements has a great impact on the FIGURE 3 | PEI success rate of the IEEE 14-bus system using empirical MEVs (Case 1). performance of the PEI method proposed in this paper, while the measurement redundancy has less impact on this method. In addition, the success rates gradually increase with the increase of the number of measurement snapshots used. Figure 3 shows the success rates of PEI using empirical values of MEVs, taking Case 1 as an example, i.e. e PF = e QF = 0.005 p.u., e PI = e QI = 0.010 p.u., e VM = 0.002 p.u. It can be seen that the maximum success rate is only 85%. Compared with Figure 2A, the maximum decline of success rate is more than 20%, and the success rate become lower when more parameter errors exist, which verifies the advantages of using the LME model to estimate the MEVs.

The IEEE 30-Bus System
In this section, the IEEE 30-bus system is used to test the effectiveness of the proposed PEI method. There are 41 series branches and 82 series parameters in this system. As the true values of the loads and power flow are very small for this system, small absolute errors for measurements are added. The following cases are designed: Case 1: e PF = e QF = 0.3%, e PI = e QI = 0.5%, e VM = 0.2%, η = 4.31. Case 2: e PF = e QF = 0.3%, e PI = e QI = 0.5%, e VM = 0.2%, η = 2.92. Similarly, 100 different error parameter combinations are selected to test the performance of the proposed method. The estimated MEVs of the two cases are summarized in Table 4. Then, based on the estimated MEVs in Table 4, the GLS is carried out for PEI. For Case 1, one scenario (N e = 6) is randomly selected from 100 combinations and the p-values of all 82 line parameters using 10-100 snapshots are shown in Figure 4. It can be seen that all erroneous parameters except for r 3-4 are identified. The erroneous parameter r 3-4 is missed because the resistance value of this line is too small and the objective function has low sensitivity to it.
The success rate results of the two cases using at most 100 snapshots are shown in Figure 5, where N e can be 2, 3, 4, 5 or 6. For Case 1, Figure 5A shows that the success rate increases if more measurement snapshots are used from the overall perspective. Specifically, when 100 snapshots are used, all the success rates of PEI are above 95%. For Case 2 with high precision but low redundancy, Figure 5B shows the success rate is greatly affected. When 10 measurement snapshots are used, the minimum success rate is only 11% (N e = 6), and the highest success rate is only 48% (N e = 2). In addition, the success rate does not show a monotonous increase as N increases. When N is 100, mostly success rates concentrated at 70%, which is far worse than Case 1.
However, the PEI results show much better performance if the definition of "success" is relaxed, i.e. for each PEI test, the result is defined as "success" if all the erroneous parameters' p-values "but one" rank top k ones among the whole N p parameters' p-values sorted in ascending order. In other words, missing of one erroneous parameter is permitted and still viewed as "success". The success rates under this definition for Case 2 are shown in Figure 6. We can see much higher success rates can be achieved than in Figure 5B. For example, the success rate rises from 48 to 91% when N = 10, N e = 6, which means that in most tests only one erroneous parameter is missed by the proposed method due to this missed parameter's low sensitivity. Figure 7 shows the PEI results using empirical MEVs, taking Case 1 as an example. It can be seen that the success rate of PEI is less than 30% when 10 snapshots is used, and the success rate increases with the increase of N. Even when N is 100, the maximum success rate is only 65%. Compared with Figure 5A, the success rate is reduced by 30%, which proves the superiority of using LME method for MEV estimation.  In this section, the proposed method is evaluated on the IEEE 118-bus system with 186 series branches and 372 series parameters. The following cases are designed for testing: Case 1: e PF = e QF = 0.5%, e PI = e QI = 1.0%, e VM = 0.2%, η = 4.67. Case 2: e PF = e QF = 1.0%, e PI = e QI = 1.5%, e VM = 0.5%, η = 4.67. Case 3: e PF = e QF = 0.5%, e PI = e QI = 1.0%, e VM = 0.2%, η = 3.09. Case 4: e PF = e QF = 1.0%, e PI = e QI = 1.5%, e VM = 0.5%, η = 3.09.
Similarly, in Case 1 and Case 2, full measurement redundancy is assumed, resulting in a high η value. In Case 3 and Case 4, we delete all the real/reactive power flow measurements at the end of each branch, resulting in a lower η value. Three measurement snapshots are used for MEV estimation, which is then used for the GLS-based PEI. The estimated MEVs of the four cases are summarized in Table 5.  For Case 1, one scenario (N e = 10) is randomly selected from 100 combinations and the p-values of all 372 line parameters are shown in Figure 8. It can be seen that all the erroneous parameters rank first and are correctly identified. Note that the p-value curve of parameter x 100-103 is too small to be displayed.
The success rate results of the four cases are shown in Figure 9, where N e can be 2-10. For Case 1 and Case 3, Figures 9A,C show that even if a large number of error parameters exist, the success rates can still reach a high level as long as the measurement accuracy is high, and the success rate will gradually increase with the increase of N. For Case 2 and Case 4 with lower measurement accuracy, Figures  9B,D show that the number of successful PEI is greatly reduced, especially for Case 4 with both lower measurement accuracy and lower measurement redundancy.
Then, if we relax the definition of "success", i.e. missing of one or more erroneous parameter is permitted and still viewed as "success". Figures 10, 11 show that much better performance can be achieved, compared with Figures 9B,D. This means that in most tests only one or two erroneous parameters are missed while most erroneous parameters are actually correctly identified by the proposed method. Figure 12 shows the success rates of PEI using empirical MEVs taking Case 1 as an example. By comparing the results of Figure 9A and Figure 12, it can be concluded that the MEVs estimated by solving the LME model result in much better PEI performance than empirical ones.

Computation Time
All the tests above are performed on a personal computer with a 3.70 GHz i7-8700K CPU and 32 GB RAM. Table 6 shows the computation times of different methods under different scenarios. It can be seen that the proposed LME-GLS method greatly saves the calculation time than the original LME-based PEI method. The main reason is that this method only uses a small number of measurement snapshots to estimate the MEVs, which avoids the computational burden of solving a very largescale LME model. In addition, the GLS only involves inversion of low-dimension matrices thus is very fast without losing the PEI accuracy.

CONCLUSION
This paper proposes an efficient method for the simultaneous identification of multiple network parameter errors based on an LME model and the GLS method. The residual equations from multiple-snapshot state estimation are used to formulate the LME model in which the parameter errors are considered as the fixed effects and the measurement errors are considered as the random effects. Then, using the measurement error variances estimated by solving the LME model, the GLS is used to estimate the parameter errors along with a hypothesis testing to infer whether each parameter error is zero. The key advantage of the proposed methodology lies in FIGURE 12 | PEI success rate of the IEEE 118-bus system using empirical MEVs (Case 1). that the LME model is only used to estimate the variances of the measurement errors using a small number of measurement snapshots, therefore the huge computation burden needed to solve a large-scale LME model is avoided. In addition, the GLS only involves inversion of low-dimension matrices, which is very efficient even a large number of measurement snapshots are used. The proposed method is shown to be promising for the simultaneous identification of multiple parameter errors. It can be used periodically by system operators in an offline manner using multiple measurement snapshots, to identify and correct multiple erroneous parameters once large residuals are reported by state estimator and existence of parameter errors are suspected. However, parameter errors are treated as unknown constants in this paper rather than random variables. Considering continuous varying of parameters during various operating conditions and load levels, future work will be focused on dynamic parameter tracking using phasor measurements.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusion of this article will be made available by the authors, without undue reservation.

AUTHOR CONTRIBUTIONS
DL contributed toward supervision, conceptualization, and writing-review and editing. YC contributed toward methodology, software, data curation, and writing-original draft. SS contributed toward software and data curation. XW contributed toward writing-review and editing. LZ contributed toward writing-review and editing.