Transmission Line Parameter Estimation Under Non-Gaussian Measurement Noise

A medium-length positive sequence transmission line model has been considered for the analysis conducted here. However, with appropriate modifications, the proposed formulation can also be applied to other line models. The line segment is modeled as a $\pi$-section (see Fig. 1), where $p$ and $q$ denote its sending and receiving ends, respectively.

Applying Kirchhoff’s circuit laws to the transmission line model shown in Fig. 1, the following equation can be written.

In (1), $I$ and $V$ denote the complex current and voltage measurements obtained from PMUs placed at the two ends of the line (represented by $p$ and $q$), $b_{pq}\in\mathbb{R}$ denotes the shunt admittance (susceptance) present at each end, and $y_{pq}\in\mathbb{C}$ denotes the series admittance, which is inverse of the series impedance, $z_{pq}\in\mathbb{C}$. The real and imaginary parts of $z_{pq}$ are the resistance, $r_{pq}\in\mathbb{R}$, and reactance, $x_{pq}\in\mathbb{R}$, of the line. The objective of this work is to estimate $r_{pq}$, $x_{pq}$, and $b_{pq}$ from the PMU measurements. Expressing the complex currents, voltages, and line parameters of (1) in their Cartesian form and rearranging the terms, we get

where, $y_{pq_{r}}\in\mathbb{R}$ and $y_{pq_{i}}\in\mathbb{R}$ denote the real and imaginary parts of $y_{pq}$. It can be observed from (2) that there is a summation of voltages in the right-hand side. This is a concern when non-Gaussian noise is present in the voltage measurements because it might be more difficult to approximate summation of two non-Gaussian noises by a GMM in comparison to approximating the individual noises (by GMMs). To alleviate this concern, linear algebra is used to express the voltages independently. This results in the following equation, where $t$ denotes the time instant at which the phasor measurements are taken.

In (3), $Y_{1}=y_{pq_{r}}$ $\in\mathbb{R}$, $Y_{2}=-(b_{pq}+y_{pq_{i}})$ $\in\mathbb{R}$, $Y_{3}=-y_{pq_{r}}$ $\in\mathbb{R}$, and $Y_{4}=y_{pq_{i}}$ $\in\mathbb{R}$, respectively. Once $\begin{bmatrix}Y_{1}&Y_{2}&Y_{3}&Y_{4}\end{bmatrix}^{T}$ are estimated, the line parameters, namely, $r_{pq}$, $x_{pq}$, and $b_{pq}$, can be calculated using (4).

However, in presence of noisy data, measurements obtained at a single time instant are not able to give accurate estimates of $\begin{bmatrix}Y_{1}&Y_{2}&Y_{3}&Y_{4}\end{bmatrix}^{T}$. Therefore, the goal is to estimate the optimal values of the line parameters from noisy measurements available from distinct time instants. Based on (3), the measurements from $s$ time instants can be concatenated to obtain a system of equation as

Essentially, (5) takes advantage of the fact that the line parameters change at a much slower rate than the speed at which the PMU measurements become available [24]. Equation (5) can now be expressed as a classical parameter estimation problem of the form shown below

where, $c\in\mathbb{R}^{n\times 1}$ denotes the dependent variables (current measurements), $D\in\mathbb{R}^{n\times p}$ denotes the independent variables (voltage measurements), and $x\in\mathbb{R}^{p\times 1}$ is the unknown parameter to be estimated (function of line parameters); the symbol ^∗ in the superscript indicates true values. It can be realized from (5) and (6), that $n=4s$ and $p$, the number of parameters to be estimated, is $4$.

In presence of Gaussian noise in PMU measurements, two scenarios arise. If the noise in the voltage measurements (independent variables) is ignored (i.e., Gaussian noise is only present in the current measurements (dependent variables)), the optimal estimate of the parameters can be obtained using LS. This scenario can be mathematically represented as

where, $c_{e}\in\mathbb{R}^{n\times 1}$ is the noise in $c$. For the linear system described by (6) and (7), the LS parameter estimate, $\hat{x}_{LS}$, is given by

In the second scenario, Gaussian noise is present in both the dependent as well as the independent variables, and the optimal parameter estimate is obtained using TLS. This scenario can be mathematically represented as

where, $D_{e}\in\mathbb{R}^{n\times p}$ is the noise in $D$. To obtain the TLS solution of the linear system described by (6) and (9), the concatenated noisy measurement matrix, $[D~{}~{}c]$, is factorized using singular value decomposition (SVD) into a matrix of singular values and left and right singular vectors. Then, using the Eckhart-Young Theorem, the TLS solution can be obtained as [25]

where, $v_{pq}$ is the vector of first $p$ elements and $v_{qq}$ is the $(p+1)^{th}$ element, respectively, of the $(p+1)^{th}$ column of the matrix of right singular vectors of the SVD of $[D~{}~{}c]$.

In this paper, we investigate both of the scenarios mentioned above under the condition that the PMU measurements have non-Gaussian noise. That is, in the analysis conducted henceforth, $c_{e}$ and $D_{e}$ will have non-Gaussian distributions. The proposed approach for solving the resulting problems is composed of two steps: a noise estimation step and a parameter estimation step. In the noise estimation step, the non-Gaussian noise is approximated by an appropriate GMM. The GMM parameters are obtained using expectation-maximization (EM). More details about this step are provided in Section III-C. The parameter estimation step is explained for the first scenario in Section III-D and for the second scenario in Section III-E.

We use $c_{e}$ to demonstrate the noise estimation step. If an $m$ component GMM is suitable for approximating $c_{e}$, then the following equation holds true.

In (11), $w_{g}\geq 0$, $\sum\limits_{g=1}^{m}w_{g}=1$, and $\mathcal{N}(c_{e};\mu,\Sigma)=\frac{1}{\sqrt{(2\pi)^{d}|\Sigma|}}\\ exp(-\frac{1}{2}(c_{e}-\mu)^{T}\Sigma^{-1}(c_{e}-\mu))$. Here, $w_{g},\mu_{g}$, and $\Sigma_{g}$ denotes the weight, mean vector, and covariance matrix of the $g^{th}$ Gaussian component of the GMM, respectively. These three variables are known as the GMM parameters. A set of these GMM parameters for the $g^{th}$ Gaussian component is denoted by $\phi_{g}\triangleq\{w_{g},\mu_{g},\Sigma_{g}\}$, and the entire set of GMM parameters to be estimated is denoted by $\theta=\{\phi_{g}\}_{g=1}^{m}$. Ideally, the optimal GMM parameters can be estimated by maximizing the log-likelihood of (11) with respect to the GMM parameters. However, the log-likelihood has summation inside the logarithm, which makes the problem of directly maximizing the log-likelihood hard. EM overcomes this problem by introducing a new variable, $z$, and maximizing the complete data log-likelihood, as explained below.

Let the $n$ independent and identically distributed samples of $c_{e}$ be denoted by $\mathcal{S}$, i.e., $\mathcal{S}=\{c_{e_{1}},c_{e_{2}},\dots,c_{e_{n}}\}$. The identity of the Gaussian component to which each of these data points belong (out of $m$ possible Gaussian components) is defined as the cluster membership of those data points. Let $z=[z_{1},z_{2},\dots,z_{m}]^{T}$ be an $m$-dimensional binary random variable. The $m$ possible values that $z$ can take can be denoted by $\{\zeta_{g}\}_{g=1}^{m}$, where $\zeta_{g}$ denotes an $m$-dimensional vector of zeros with one in the $g^{th}$ place. The probability of $z$ taking value $\zeta_{g}$ is $w_{g}$. In other words, $\mathbb{P}(z=\zeta_{g})=w_{g}$. Thus, the distribution of $z$ can be written as

The joint distribution of $c_{e}$ and $z$ can now be expressed as

while the conditional probability of $z$ given $c_{e}$ is given by

Let $\mathcal{S}_{c}=\{(c_{e_{i}},z_{i})\}_{i=1}^{n}$. Now, EM maximizes $\mathcal{S}_{c}$ instead of $\mathcal{S}$, where the likelihood of $\mathcal{S}_{c}$ can be written as [26]

In (15), $z_{g}^{i}$ represents $g^{th}$element of $z$ for the $i^{th}$ data sample. The complete data log-likelihood can be obtained by taking logarithm on both sides of (15), as shown below.

From (16), EM can learn the GMM parameters for a given noise vector. EM has two main steps - expectation and maximization - which are repeated until convergence. The expectation step computes the conditional expectation of the complete data log-likelihood using given GMM parameters. In the maximization step, this conditional expectation of the complete data log-likelihood is maximized to obtain updated GMM parameters. After EM converges, the noise vector is clustered into different bins based on the cluster membership. These three processes are described below.

For the scenario in which non-Gaussian noise (expressed as a GMM) is present only in the dependent variables, the parameter estimation can be formulated as a minimization problem whose objective is to minimize the SSSE of the noise variable and the constraints are the error definitions. This is mathematically described by

Using a vector of Lagrange multipliers, $\lambda_{g}\in\mathbb{R}^{n_{g}\times 1}$, for each Gaussian component, (26) is expressed as an unconstrained optimization problem as shown below

The optimal parameter estimate is found by computing the $\hat{x}$ that minimizes (27). This is done by simultaneously solving $\partialderivative{\phi}{c_{e_{g}}}=0,g\in[m]$, $\partialderivative{\phi}{\lambda_{g}}=0,g\in[m]$, and $\partialderivative{\phi}{x}=0$, to yield the following equation for the parameter estimate:

Equation (28) can be directly used to compute the optimal parameter estimate when the noise characteristics and cluster membership of data points of the GMM is available. In the actual setting where the measurement noise characteristics are unknown, the noise estimation procedure explained in Section III-C must be used in conjunction with this parameter estimation step and iterated until convergence, leading to a joint estimation of noise and parameters. However, Section III-C assumes prior knowledge of $m$. For instance, Ahmad et al. have shown in [12] that PMU measurement noise (particularly the noise in the current phasors) can be represented using a three component GMM. However, as the power system is non-stationary, there is no guarantee that $m$ will always be equal to three. A strategy to circumvent the need for the a priori knowledge of $m$ is proposed below.

In absence of the prior knowledge of $m$, the proposed methodology for joint estimation of noise and parameters must be repeated for a range of values of $m$, say, $m=1$ to $m=m_{max}$. A reasonable value of $m_{max}$ could be ten. Then, an information theory-based model selection criterion called the Bayesian information criterion (BIC) can be used to identify the most appropriate value for $m$ [28]. The BIC has two components - one component rewards the goodness of fit of the model, while the other component penalizes the model complexity. Thus, the optimal number of GMM components will be the $m$ for which the lowest value of BIC is obtained. Note that by removing the need for the a priori knowledge of $m$, this strategy further improves the robustness of the proposed approach. The flowchart describing the overall procedure is presented in Fig. 2.

In Fig. 2, for each value of $m$, the expectation and maximization steps learn the characteristics of the GMM that approximates the noise vector. The cluster membership is used to cluster the data points in the noise vector into $m$ Gaussian components. These, then become inputs to the parameter estimation step, which is carried out using (28). For each value of $m$, noise and parameters estimates in presence of non-Gaussian noise in the dependent variables are obtained in an iterative manner. The optimal estimates of noise and parameters correspond to the value of $m$ that gives the lowest value of BIC. This concludes the flowchart.

The second scenario accounts for measurement noise in both the dependent as well as the independent variables, and are known as EIV problems [29]. The measurement model for this scenario is mathematically described by (6) and (9). The noise estimation methodology described in Section III-C is leveraged here to solve the EIV problem in presence of non-Gaussian noise in both dependent and independent variables. The $g^{th}$ Gaussian component of the standardized dependent variable noise vector is denoted as

Similarly, the $g^{th}$ Gaussian component of the standardized independent variable noise vector can be denoted as

Note that $c_{e_{g_{N}}}^{T}c_{e_{g_{N}}}$ denotes the SSSE in $c_{e_{g}}$. However, as $D_{e_{g_{N}}}$ is a matrix, it has to be first converted to a column vector form before it can be converted to an SSSE form. Let $\text{vec}(D_{e_{g_{N}}})$ denote the vectorization operation of the matrix $D_{e_{g_{N}}}$, where the matrix is converted to a column vector by stacking the columns of the matrix on top of each other. Then, using the properties of the Kronecker product, $\otimes$, the resulting minimization problem for optimal parameter estimation for the EIV problem can be mathematically described by:

Using the method of Lagrange multipliers, (31) can be written as an unconstrained objective function as

The optimal parameter estimate for this objective function is the solution to the following non-linear system of equations:

where,

In (34), $j$ denotes the $j^{th}$ column for D and $j^{th}$ parameter for $x$. The optimal parameters are found by solving (33) using Newton’s method:

where, $\text{Jac}(f(x^{(k)}))$ denotes the Jacobian of $f(x^{(k)})$. The components of the optimal noise estimates of $c_{e}$ and $D_{e}$ can be obtained from the following equation,

The GMM noise characteristics and cluster membership can be obtained by applying EM to the noise estimates. Since the parameter estimation step uses the results from the noise estimation step and vice-versa, they are repeated until convergence. Lastly, the noise and parameter estimation is performed iteratively for $\forall m\in[1,m_{max}]$, and BIC is used to select the optimal $m$, and noise and parameter estimates. The overall process is described as an AlgorithmAlgorithm below. As this algorithm solves an EIV problem using a combination of GMMs, Lagrange multipliers, and EM, it is henceforth called the EGLE Algorithm.

The inputs for the EGLE AlgorithmAlgorithm are the noisy measurements, an initial guess for the parameters, the tolerances ($\epsilon_{0},\epsilon_{1},\epsilon_{2}$), and $m_{max}$. The tolerance values can be chosen based on the desired estimation accuracy and speed requirements. The outermost loop finds the optimal number of Gaussian components ($m^{*}$). The loop in the middle performs noise estimation using EM (see Section III-C) and provides inputs to the innermost loop. The innermost loop performs parameter estimation by solving (33) using Newton’s method. The parameter estimate obtained from the innermost loop serves as an input to the loop in the middle. The tolerance $\epsilon_{2}$ determines the termination of the parameter estimation loop, whereas $\epsilon_{1}$ determines when the joint noise and parameter estimation is completed for each value of $m$. The value of BIC corresponding to every value of $m$ is also calculated. Finally, the optimal value of $m$ is found using the index of the minimum value of $\mathrm{BIC}(m)$. The $m^{*}$ thus obtained gives the optimal number of GMM components that must be used to approximate the non-Gaussian noise, while the parameter estimate corresponding to that $m^{*}$ (namely, $x_{m}(m^{*})$) is the optimal parameter estimate, $x^{*}$.

Note that when $D$ is noise-free, then the EGLE AlgorithmAlgorithm simplifies to the flowchart shown in Fig. 2. As such, EGLE can be used when non-Gaussian noise is only present in the dependent variables as well as when non-Gaussian noise is present in both the dependent and the independent variables. The results obtained using EGLE for simulated and actual PMU data are described in the following sections.

In this sub-section, a two component GMM measurement noise of mean ($\begin{bmatrix}0&0.005\end{bmatrix}$), standard deviation ($\begin{bmatrix}0.0015&0.0015\end{bmatrix}$), and weight ($\begin{bmatrix}0.3&0.7\end{bmatrix}$) was added to the currents to create the noisy dependent variables, while the independent variables (voltages) were kept noise-free. The convergence of EGLE for the line which connects buses 38 and 65 ($L_{38-65}$) of the IEEE 118-bus system is shown in Fig. 3(a)-3(d). The figures also compare the performance of the proposed method with LS for all four $Y$ parameters. The colors magenta, red, and blue correspond to the true value, and estimates obtained from EGLE and LS, respectively. It can be observed from the figures that for each of the parameters, the estimate obtained using EGLE is: (a) very close to the true value, and (b) better than the LS estimate. It is also obvious that the superior performance in estimating the $Y$ parameters will lead to better performance in the estimation of the actual line parameters from the $Y$ values.

This sub-section corresponds to the most general setting in which both dependent and independent variables have non-Gaussian noise. The two component GMM noise that was defined in Section IV-A was also used here, except that this time it was added to both current and voltage measurements. The results corresponding to two 345 kV lines and two 138 kV lines of the IEEE 118-bus system are analyzed below. The 345 kV lines are between buses 38 and 65, and between buses 8 and 9. The 138 kV lines are between buses 47 and 69, and between buses 75 and 69. The comparison of the resistance, reactance, and susceptance estimates obtained using LS, TLS, and EGLE for these four lines is shown in Fig. 4, Fig. 5, and Fig. 6, respectively.

To draw reliable statistical inferences, the experiment was repeated 1,000 times by randomly generating noise based on the specified GMM characteristics and a random initial guess for the parameter estimates that was within $\pm 30\%$ of the true value. The $\mathrm{MARE}$ for LS, TLS, and EGLE for all the MC runs are shown as bar plots in Fig. 4-6. The $\mathrm{SDARE}$ are also calculated, a low value of which is an indication of the consistency of the results. The $\mathrm{SDARE}$ is displayed as black colored error bars in Fig. 4-6.

From Fig. 4-6, it can be observed that EGLE consistently performs better than LS and TLS when estimating line parameters in presence of GMM noise. This also translates to superior performance of EGLE in terms of the $\mathrm{MARE_{net}}$. For example, for $L_{38-65}$, the $\mathrm{MARE_{net}}$ for LS and TLS methods were $1.55\%$ and $1.33\%$, respectively, whereas the proposed method had a $\mathrm{MARE_{net}}$ of only $0.40\%$. The $\mathrm{SDARE_{net}}$ for the three methods were 0.06%, 0.06%, and 0.10% respectively, indicating consistency of the estimates. A similar trend was found for the other lines as well. It was also observed that the relative accuracies of the parameter estimates for the three algorithms across all four lines was in accordance with the conditioning of the $D^{T}D$ matrix.

The results obtained in the previous two sub-sections demonstrate the superior performance of EGLE over the conventional approaches in presence of non-Gaussian noise in the measurements. However, as TLPE is an offline process, an alternate strategy could be to first remove the non-Gaussian noise from the measurements, and then perform estimation using conventional approaches. This strategy, called the denoising-followed-by-(conventional)-estimation approach, is compared with EGLE in this sub-section.

The noisy data created in Section IV-B was used for this analysis. The denoising was done for both voltage and current phasor measurements using a moving window median absolute deviation filter that was developed in [12] to detect non-Gaussian noise in PMU measurements. Subsequently, the LS method was employed for TLPE. The results that were obtained when this approach was compared with EGLE are shown in Fig. 7. In Fig. 7, the purple bars indicate the $\mathrm{MARE}$ obtained using the denoising approach whereas the green bars show the $\mathrm{MARE}$ obtained using EGLE. It is clear from the figure that EGLE performs better than the denoising based approach, particularly for the resistance and susceptance estimates. One reason for this observation could be the sensitivity of the filter to the size of the window that was considered for analysis (the window size was set at 600 samples as mentioned in [12]). In summary, this analysis demonstrated the superiority of an integrated approach for noise and parameter estimation (EGLE) in comparison to a denoising-followed-by-(conventional)-estimation approach.