Harmonic Source Localization Approach Based on Fast Kernel Entropy Optimization ICA and Minimum Conditional Entropy

Abstract

Based on the fast kernel entropy optimization independent component analysis and the minimum conditional entropy, this paper proposes a harmonic source localization method which aims at accurately estimating harmonic currents and identifying harmonic sources. The injected harmonic currents are estimated by the fast kernel entropy optimization independent component analysis (FKEO-ICA) in the absence of prior knowledge of harmonic impedances. Then, the minimum conditional entropy is applied to locate the harmonic sources based on the estimated harmonic currents. The proposed harmonic source localization method is validated on the IEEE 34-bus system. By applying the correlation coefficient and three error evaluation indicators, comparison has been made among the performances of the FKEO-ICA and three other ICA algorithms. The results show that the FKEO-ICA algorithm could achieve a significantly better accuracy of harmonic current estimation, while the minimum conditional entropy could determine the locations of harmonic sources precisely.

1. Introduction

The growing applications of power electronic apparatus and non-linear loads can result in serious harmonic pollution in electrical power systems. As a consequence, the harmonic estimation and localization of harmonic sources have drawn wide concern globally.

Harmonic State Estimation (HSE) [1,2] involves harmonic distribution and harmonic source identification [3,4,5,6,7,8,9,10]. With regard to the HSE problems, a multitude of methods have been employed, including least squares (LS) [1,2], singular value decomposition (SVD) [4], Kalman filter [5], neural networks (NN) [6], sparsity maximization [7], and particle swarm optimization (PSO) [8]. Additionally, it has been proved that the complete harmonic distribution can be derived based on selected measurement data [5]. However, in terms of the HSE methods, precise information on network parameters are required, which are rarely known in practice. Furthermore, sufficient measurements are also required to guarantee full observability, while considerable computational resources are demanded to ensure a relatively satisfactory processing speed, especially for a large distribution system. Consequently, the applications of HSE in the power system are restricted.

To estimate harmonic sources in the absence of prior information relating to the harmonic impedances, independent component analysis (ICA) [11,12] was conducted [13,14,15,16], using one of the blind source separation (BSS) approaches [17,18], which use measured bus voltages to estimate the magnitude of harmonic currents. In this process, it is assumed that the harmonic sources are statistically independent, non-Gaussian distributed, and linearly mixed [11]. The ICA algorithms applied include fast independent component analysis (Fast-ICA) [13], efficient variant Fast-ICA (EFICA) [14], supervised independent component analysis (SICA) [15], and single-channel independent component analysis (SCICA) [16]. For harmonic current estimation, different ICA algorithms display different levels of accuracy for different harmonic frequencies. Meanwhile, the estimation error at lower load levels is, accordingly, increased with the measurement noise [19]. However, the existing ICA-based harmonic current estimation methods often ignore measurement errors. Additionally, the orders of the components estimated by ICA algorithms are usually random, in which case it is difficult to determine the exact locations of harmonic sources. Thus, it appears necessary to propose a method which can precisely estimate harmonic currents and locate harmonic sources at the most concerned frequencies with the presence of measurement noises.

This paper proposes a method based on the fast kernel entropy optimization ICA (FKEO-ICA) aiming at estimating harmonic currents, and based on the minimum conditional entropy, which aims to identify the exact locations of harmonic sources. The main advantage of this method is that it only needs to measure the bus harmonic voltage, which is more accessible and more reliable compared with other harmonic measurements. There is no requirement for the data of system harmonic impedances. Section 2 discusses the relationship between the harmonic state estimation model and the ICA model. Additionally, the basic principle of the FKEO-ICA algorithm is elaborated in Section 3. Section 4 and Section 5 present the specific process of harmonic source localization with the employment of FKEO-ICA and minimum conditional entropy. Furthermore, the estimation and localization results for the IEEE 34-bus system are shown in Section 6. Lastly, Section 7 concludes the findings.

2. The Relationship between HSE Model and ICA Model

The model of HSE at the harmonic order $h$ without noise is as follows:

$$\begin{array}{cc}{\mathit{U}}_h(t_i)={\mathit{Z}}_h{\mathit{I}}_h(t_i),& i=1,2,\cdots ,T\end{array}$$

(1)

where $t_i$ is the sampling time; $T$ is the number of samples; ${\mathit{U}}_h(t_i)$ are the measured harmonic voltage vectors which have been known previously; ${\mathit{I}}_h(t_i)$ are the unknown harmonic current vectors; and ${\mathit{Z}}_h$ is the harmonic impedance matrix, which can be obtained by ${\mathit{Z}}_h={\mathit{Y}}_h^{-1}$, where ${\mathit{Y}}_h$ is the nodal harmonic admittance matrix, as defined in [3].

BSS is to estimate P statistically-independent zero-mean signals, $\mathit{S}(N)={[S_1(N),S_2(N),\cdots ,S_P(N)]}^T$, from K measured signals, $\mathit{X}(N)={[X_1(N),X_2(N),\cdots ,X_K(N)]}^T$, where K ≥ P; N represents the number of signal samples. Generally, the model of ICA without noise [14] can be described as:

$$\mathit{X}=\mathit{M}\mathit{S}$$

(2)

where M is a $K\times P$ mixture matrix.

The ICA algorithm aims to obtain the separation matrix W so as to eliminate the mixing. That is, the estimation of the sources can be given as:

$$\mathit{Y}=\mathit{W}\mathit{X}$$

(3)

where Y is the estimation value of S.

In accordance with Equations (1) and (2), ${\mathit{U}}_h(t_i)$, ${\mathit{Z}}_h$, and ${\mathit{I}}_h(t_i)$ in the HSE model are associated with $\mathit{X}$, $\mathit{M}$, and $\mathit{S}$ in the ICA model, respectively. Therefore, the harmonic state estimation problem can be formulated as an ICA problem.

3. The Principle of the FKEO-ICA Algorithm

In view of the high dimension and large sample size of harmonic data involved in the electric power systems, the fast kernel optimization ICA method based on the FFT-based fast convolution, as presented in [20], is applied to estimate harmonic currents in this paper. The basic principle of the FKEO-ICA algorithm is described as follows.

The set of independent signal sources is denoted as $\{{\mathit{s}}_1,{\mathit{s}}_2,\cdots ,{\mathit{s}}_P\}$, where ${\mathit{s}}_p={[s_p(1),s_p(2),\cdots ,s_p(N)]}^T$ represents each individual source, $p=1,2,\cdots ,P$. With regard to the linear mixture of ${\mathit{x}}_1,{\mathit{x}}_2,\cdots ,{\mathit{x}}_P$, the set of measured signals is denoted as $\{{\mathit{x}}_1,{\mathit{x}}_2,\cdots ,{\mathit{x}}_P\}$. Additionally, $\{{\widehat{\mathit{s}}}_1,{\widehat{\mathit{s}}}_2,\cdots ,{\widehat{\mathit{s}}}_P\}$ is assumed to be the set of estimated sources.

According to Equation (3), ${[{\widehat{\mathit{s}}}_1,{\widehat{\mathit{s}}}_2,\cdots ,{\widehat{\mathit{s}}}_P]}^T=\mathit{W}{[{\mathit{x}}_1,{\mathit{x}}_2,\cdots ,{\mathit{x}}_P]}^T$ can be obtained. The aim of ICA is to calculate separation matrix W which can make the estimated sources independent. The common criterion is the mutual information of the estimated signal sources [14]. Furthermore, the mutual information of ${\widehat{\mathit{s}}}_1,{\widehat{\mathit{s}}}_2,\cdots ,{\widehat{\mathit{s}}}_P$ is expressed as follows:

$$MI({\widehat{\mathit{s}}}_1,{\widehat{\mathit{s}}}_2,\cdots ,{\widehat{\mathit{s}}}_P)={\displaystyle \sum _{p=1}^P{H}_{{\widehat{\mathit{s}}}_p}-\log \left|\det (\mathit{W})\right|}-H_C$$

(4)

where $H_{{\widehat{\mathit{s}}}_p}$ is the entropy of estimated signal ${\widehat{\mathit{s}}}_p$, and $H_C$ is a constant, generally, which can be ignored in the optimization process. Accordingly, the minimization problem to be solved becomes:

$$\underset{\mathit{W}}{\min }\left\{{\displaystyle \sum _{p=1}^P{H}_{{\widehat{\mathit{s}}}_p}-\log \left|\det (\mathit{W})\right|}+\mathsf{\beta }{\displaystyle \sum _{p=1}^P{(\Vert {\widehat{\mathit{s}}}_p-1\Vert )}^2}\right\}$$

(5)

where $\mathsf{\beta }$ is a constant weight, $\mathsf{\beta }{\displaystyle \sum _{p=1}^P{(\Vert {\widehat{\mathit{s}}}_p-1\Vert )}^2}$ is a weighted sum which is applied to solve the ambiguities issue from the scale invariance of mutual information, and its impact for the calculation accuracy is negligible [20].

Based on Parzen-windows (PW), the estimation of $H_{{\widehat{\mathit{s}}}_p}$ is capable of differentiating estimated entropies, while an approximate expression is required for the differential entropy gradients [20]. For a given r, the PW estimator for the probability density function is:

$$\widehat{p}(r|{\widehat{\mathit{s}}}_p)\equiv (1/N){\displaystyle \sum _{n=1}^N\zeta [r-{\widehat{\mathit{s}}}_p(n)]}$$

(6)

where ${\widehat{\mathit{s}}}_p(n)$ is sampled from estimated signal ${\widehat{\mathit{s}}}_p$ and $\zeta [\cdot ]$ represents the smoothing kernel function.

The PW estimator for $H_{{\widehat{\mathit{s}}}_p}$ is:

$${\widehat{H}}_{{\widehat{\mathit{s}}}_p}=-\frac{1}{N}{\displaystyle \sum _{t=1}^N\log \widehat{p}[{\widehat{\mathit{s}}}_p(t)|{\widehat{\mathit{s}}}_p]}=-\frac{1}{N}{\displaystyle \sum _{t=1}^N\log \left\{\frac{1}{N}{\displaystyle \sum _{n=1}^N\zeta [{\widehat{\mathit{s}}}_p(t)|{\widehat{\mathit{s}}}_p(n)]}\right\}}$$

(7)

The gradient of $\log \left|\det (\mathit{W})\right|$ is ${({\mathit{W}}^{-1})}^T$ [16]. Thus, the gradient of $MI({\widehat{\mathit{s}}}_1,{\widehat{\mathit{s}}}_2,\cdots ,{\widehat{\mathit{s}}}_P)$ can be expressed as:

$${\nabla }_{\mathit{W}}MI({\widehat{\mathit{s}}}_1,{\widehat{\mathit{s}}}_2,\cdots ,{\widehat{\mathit{s}}}_P)={\nabla }_{\mathit{W}}{\displaystyle \sum _{p=1}^P{H}_{{\widehat{\mathit{s}}}_p}-{({\mathit{W}}^{-1})}^T}$$

(8)

where ${\nabla }_{\mathit{W}}{\displaystyle \sum _{p=1}^P{H}_{{\widehat{\mathit{s}}}_p}}$ can be calculated from:

$$\begin{array}{l}{\nabla }_{\mathit{W}}{\displaystyle \sum _{p=1}^P{H}_{{\widehat{\mathit{s}}}_p}}={\left[{\nabla }_{{\widehat{\mathit{s}}}_1}{H}_{{\widehat{\mathit{s}}}_1},{\nabla }_{{\widehat{\mathit{s}}}_2}{H}_{{\widehat{\mathit{s}}}_2},\cdots ,{\nabla }_{{\widehat{\mathit{s}}}_P}{H}_{{\widehat{\mathit{s}}}_P}\right]}^T\left[{\mathit{x}}_1,{\mathit{x}}_2,\cdots ,{\mathit{x}}_P\right]\\ {\nabla }_{{\widehat{\mathit{s}}}_p}={\left[\frac{\partial H_{{\widehat{\mathit{s}}}_p}}{\partial {\widehat{\mathit{s}}}_p(1)},\frac{\partial H_{{\widehat{\mathit{s}}}_p}}{\partial {\widehat{\mathit{s}}}_p(2)},\cdots ,\frac{\partial H_{{\widehat{\mathit{s}}}_p}}{\partial {\widehat{\mathit{s}}}_p(N)}\right]}^T\end{array}$$

(9)

The derivatives of Equation (7) are:

$$\begin{array}{l}\frac{\partial H_{{\widehat{\mathit{s}}}_p}}{\partial {\widehat{\mathit{s}}}_p(l)}=-\frac{1}{N}{\displaystyle \sum _{t=1}^N\frac{(1/N){\displaystyle \sum _{n=1}^N{\zeta }^{\prime }[{\widehat{\mathit{s}}}_p(t)-{\widehat{\mathit{s}}}_p(n)]}[{\delta }_{tl}-{\delta }_{nl}]}{(1/N){\displaystyle \sum _{n=1}^N\zeta [{\widehat{\mathit{s}}}_p(t)-{\widehat{\mathit{s}}}_p(n)]}}}\\ \begin{array}{ccc}& & =\end{array}-\frac{1}{N}\frac{(1/N){\displaystyle \sum _{n=1}^N{\zeta }^{\prime }[{\widehat{\mathit{s}}}_p(l)-{\widehat{\mathit{s}}}_p(n)]}}{\widehat{p}[{\widehat{\mathit{s}}}_p(l)|{\widehat{\mathit{s}}}_p]}+\frac{1}{N}{\displaystyle \sum _{t=1}^N\frac{(1/N){\zeta }^{\prime }[{\widehat{\mathit{s}}}_p(t)-{\widehat{\mathit{s}}}_p(l)]}{\widehat{p}[{\widehat{\mathit{s}}}_p(t)|{\widehat{\mathit{s}}}_p]}}\end{array}$$

(10)

where ${\zeta }^{\prime }[\cdot ]$ represents the derivatives of $\zeta [\cdot ]$, and ${\delta }_{tl}$ represents the kroneker delta, ${\delta }_{tl}=\{\begin{array}{l}\begin{array}{ccc}0& if& t\ne l\end{array}\\ \begin{array}{ccc}1& if& t=l\end{array}\end{array}$.

Based on the above-mentioned derivation process, the gradient equations have been converted into accurate analytical formulae, and the separation matrix W can, thus, be obtained. Further acceleration based on the FFT-based fast convolution is described in [20], which has proved that the estimation error is reasonably small. Moreover, through the FKEO-ICA algorithm, both high computational accuracy and low computational complexity can be achieved.

The harmonic current estimation model using FKEO-ICA is illustrated in Figure 1.

4. Identification of Harmonic Source Localization Using Minimum Conditional Entropy

Similar to other ICA methods, FKEO-ICA is not capable of determining the exact bus of harmonic sources. The relationship between the harmonic currents and the bus voltages is $\begin{array}{cc}{\mathit{U}}_h(t_i)={\mathit{Z}}_h{\mathit{I}}_h(t_i),& i=1,2,\cdots ,T\end{array}$, which indicates that the injected harmonic current ${\mathit{I}}_h$ produced by the harmonic source at the injected bus has a certain correlation with the bus harmonic voltage ${\mathit{U}}_h$. As a consequence, the bus harmonic voltage ${\mathit{U}}_h$ and the injected harmonic current ${\mathit{I}}_h$ can be considered as two random variables with associated distributions. In accordance with the information theory [21,22,23], the conditional entropy has been employed to measure the correlation between two random variables in case one of them is known. The pair-wise conditional entropy between the estimated harmonic current ${\tilde{I}}_{h\_est}$ and the measured bus voltage $U_h$ at each frequency can be expressed as follows:

$$H({\tilde{I}}_{h\_est}|U_h)={\displaystyle \sum _{t_i=1}^{T}p(U_h(t_i),{\tilde{I}}_{h\_est}(t_i))\log }[p_r(U_h(t_i))/p(U_h(t_i),{\tilde{I}}_{h\_est}(t_i))]$$

(11)

where $p(U_h(t_i),{\tilde{I}}_{h\_est}(t_i))$ is the joint probability distribution of the measured bus voltage and the estimated harmonic current; $p(U_h(t_i),{\tilde{I}}_{h\_est}(t_i))=p({\tilde{I}}_{h\_est}(t_i)|U_h(t_i))p(U_h(t_i))$, where $p({\tilde{I}}_{h\_est}(t_i)|U_h(t_i))$ is the probability of the estimated harmonic current under the given measured bus voltage; $p_r(U_h(t_i))$ is the marginal probability distribution of the measured bus voltage averaging over information about the estimated harmonic current, and it can be calculated by summing the joint probability distribution over the estimated harmonic current.

The stronger the correlation between the two random variables, the smaller their conditional entropy is. Owing to the current division effect between branches, the injected harmonic current has minimum conditional entropy with the voltage of its injected bus, compared with the conditional entropies with other bus voltages. Thus, The minimum pair-wise conditional entropy $H({\tilde{\mathit{I}}}_{h\_est}|{\mathit{U}}_h)$ can be used to determine the harmonic source location, where ${\tilde{\mathit{I}}}_{h\_est}$ is the normalized estimated harmonic current, and ${\mathit{U}}_h$ is the normalized measured bus voltage. The principle of harmonic source localization based on the minimum pair-wise conditional entropy is shown in Figure 2.

5. Harmonic Source Localization Using FKEO-ICA and Minimum Conditional Entropy

The process of the harmonic source localization based on FKEO-ICA and minimum conditional entropy is as follows:

(1) Measurement of bus harmonic voltages

The harmonic voltages at all buses are measured and represented as vector ${\mathit{U}}_h$.

(2) Centralization of the harmonic voltage data

The first pre-processing procedure is centralization, which aims to transform vector ${\mathit{U}}_h$ into the zero-mean vector ${\mathit{U}}_{hC}$ by subtracting the mean of ${\mathit{U}}_h$.

(3) Whitening of the centralized data

The second pre-processing procedure is to whiten the centralized data. This means that the centralized vector ${\mathit{U}}_{hC}$ can be linearly transformed into ${\tilde{\mathit{U}}}_{hC}$. In the process, the components of ${\tilde{\mathit{U}}}_{hC}$ are independent, while the variances of ${\tilde{\mathit{U}}}_{hC}$ are normalized into a unity. The linear whitening transformation can be defined as:

$${\tilde{\mathit{U}}}_{hC}=(\mathit{O}{\mathit{\Lambda }}^{-1/2}{\mathit{O}}^T){\mathit{U}}_{hC}$$

(12)

where $\mathit{O}=(o_1,o_2,\cdots ,o_n)$ indicates the orthogonal matrix with vectors ${(o_1,o_2,\cdots ,o_n)}^T$ as the unit norm eigenvectors of the covariance matrix relating to ${\mathit{U}}_{hC}$, and $\mathit{\Lambda }=diag({\mathsf{\Lambda }}_1,{\mathsf{\Lambda }}_2,\cdots ,{\mathsf{\Lambda }}_n)$ is the diagonal matrix of the eigenvalues concerned with the covariance matrix of ${\mathit{U}}_{hC}$.

(4) Determination of the number of independent components

Prior to the estimation of harmonic currents, the number of independent components needs to be identified through the principal component analysis [24].

(5) Estimation of harmonic currents

The FKEO-ICA algorithm elaborated in Section 3 can be used to obtain the separation matrix W. Then the harmonic currents ${\mathit{I}}_{h\_est}$ can be estimated based on Equation (3).

(6) Localization of the harmonic sources

The last procedure is to calculate the conditional entropies $H({\tilde{\mathit{I}}}_{h\_est}|{\mathit{U}}_h)$ between the normalized estimated currents and the normalized voltages based on Equation (11), and find the precise location of harmonic sources based on the minimum conditional entropy.

In summary, the process of the proposed algorithm for harmonic source localization can be illustrated in Figure 3.

6. Example Test

The IEEE 34-bus power system [25] as shown in Figure 4 is selected as the test system, with the system parameters listed in the Appendix.

In the simulation, the fundamental frequency of the system is 50 Hz. The harmonic loads are modeled as constant power loads in the fundamental frequency power flow calculations. Then the power flow solution, together with the spectrum of the six-pulse HVDC as shown in

Table 1. The harmonic spectrum of the six-pulse HVDC.

Harmonic Order	Magnitude (p.u.)	Angle (deg.)
1	1.0000	−49.56
5	0.1941	−67.77
7	0.1309	11.90
11	0.0758	−7.13
13	0.0586	68.57
17	0.0379	46.53

, are used to obtain the harmonic current source models. In our work, the 5th, 7th, and 11th harmonics are selected for the test, as the amplitude of them is larger than that of the higher order harmonics. Then, the three harmonic sources containing the 5th, 7th, and 11th harmonic orders are located at bus 12, 30, and 34, respectively. To simulate altered operating situations, all loads are multiplied with the random variables, which are statistically independent and obey a Laplace distribution with the variance of 0.002. Furthermore, the Backward/Forward Sweep-based algorithm [26] has been applied to calculate harmonic power flow and generate harmonic voltages. To take measurement noises into account, the white Gaussian noise with SNR = 60 dB is added to the harmonic power flow data using the “AWGN” function in the Communications System Toolbox of MATLAB. In this study, 600 samples of RMS harmonic bus voltages and injected harmonic currents for the 34 buses were created.

Comparison has been made among the estimation performances of FKEO-ICA, Fast-ICA [17], equivalent robust ICA (ERICA) [27], and unbiased qNewton ICA (UNICA) [27]. All algorithms are implemented on MATLAB R2014a in an Intel Pentium Dual Core, 2.0 GHz, 1 GB RAM computer. The correlation coefficients (CC) between the actual and estimated currents are calculated as shown in

Table 2. CC between actual and estimated currents.

Harmonic Order	Bus Number	FKEO-ICA	Fast-ICA	ERICA	UNICA
5th harmonic	12	0.9983	0.9961	0.9200	0.8690
	30	0.9975	0.9929	0.9821	0.9471
	34	0.9990	0.9980	0.9665	0.8942
7th harmonic	12	0.9984	0.9962	0.9193	0.8671
	30	0.9973	0.9927	0.9809	0.9451
	34	0.9991	0.9977	0.9651	0.8943
11th harmonic	12	0.9984	0.9966	0.9186	0.8653
	30	0.9976	0.9927	0.9795	0.9431
	34	0.9990	0.9978	0.9662	0.8942

. Meanwhile, the errors are measured by:

$$MAE=\frac{1}{n}{\displaystyle \sum _{t=1}^T\left|I_{h\_est}(t)-I_{h\_act}(t)\right|}$$

(13)

$$MSE=\frac{1}{n-1}{\displaystyle \sum _{t=1}^T{(I_{h\_est}(t)-I_{h\_act}(t))}^2}$$

(14)

$$MAPE=\frac{1}{n}{\displaystyle \sum _{t=1}^T(\left|I_{h\_est}(t)-I_{h\_act}(t)\right|/I_{h\_act}(t))}$$

(15)

where $I_{h\_est}(t)$ and $I_{h\_act}(t)$ represent the estimated and actual current values under the same conditions of time and frequency, respectively; MAE indicates the mean absolute error; MSE indicates the mean squared error; and MAPE indicates the mean absolute percentage error. The errors of four ICA algorithms are shown in

Table 3. MAE between actual and estimated currents.

Harmonic Order	Bus Number	FKEO-ICA	Fast-ICA	ERICA	UNICA
5th harmonic	12	7.5904 × 10−3	1.5510 × 10−2	7.3730 × 10−2	1.0639 × 10−1
	30	9.2604 × 10−3	1.5313 × 10−2	3.6112 × 10−2	7.5313 × 10−2
	34	5.1122 × 10−3	6.8862 × 10−3	4.4565 × 10−2	8.1041 × 10−2
7th harmonic	12	6.8337 × 10−3	1.6352 × 10−2	7.5081 × 10−2	1.0790 × 10−1
	30	9.9669 × 10−3	1.6092 × 10−2	3.5680 × 10−2	7.5697 × 10−2
	34	4.7057 × 10−3	7.3396 × 10−3	4.6096 × 10−2	8.1386 × 10−2
11th harmonic	12	7.1430 × 10−3	1.4016 × 10−2	7.5491 × 10−2	1.0900 × 10−1
	30	8.7369 × 10−3	1.5055 × 10−2	3.5280 × 10−2	7.5760 × 10−2
	34	5.4324 × 10−3	7.3482 × 10−3	4.1690 × 10−2	7.8008 × 10−2

,

Table 4. MSE between actual and estimated currents.

Harmonic Order	Bus Number	FKEO-ICA	Fast-ICA	ERICA	UNICA
5th harmonic	12	9.4362 × 10−5	3.7155 × 10−4	8.1654 × 10−3	1.653 × 10−2
	30	1.4417 × 10−4	3.9838 × 10−4	1.7862 × 10−3	7.4094 × 10−3
	34	5.0132 × 10−5	8.4152 × 10−5	2.7657 × 10−3	9.0424 × 10−3
7th harmonic	12	9.3467 × 10−5	4.0118 × 10−4	8.4285 × 10−3	1.6989 × 10−2
	30	1.6414 × 10−4	4.3235 × 10−4	1.7705 × 10−3	7.5194 × 10−3
	34	4.4444 × 10−5	9.2396 × 10−5	2.9502 × 10−3	9.1054 × 10−3
11th harmonic	12	1.0107 × 10−4	3.0905 × 10−4	8.5180 × 10−3	1.7331 × 10−2
	30	1.5514 × 10−4	3.8765 × 10−4	1.7543 × 10−3	7.5693 × 10−3
	34	5.9595 × 10−5	9.0453 × 10−5	2.4558 × 10−3	8.4541 × 10−3

and

Table 5. MAPE between actual and estimated currents.

Harmonic Order	Bus Number	FKEO-ICA	Fast-ICA	ERICA	UNICA
5th harmonic	12	3.0513 × 10−5	6.0647 × 10−5	2.8736 × 10−4	4.1395 × 10−4
	30	3.6046 × 10−5	5.8307 × 10−5	1.4993 × 10−4	3.0977 × 10−4
	34	2.5417 × 10−5	2.8354 × 10−5	1.6142 × 10−4	2.9292 × 10−4
7th harmonic	12	3.2119 × 10−5	6.4991 × 10−5	2.9288 × 10−4	4.1988 × 10−4
	30	3.8687 × 10−5	6.1150 × 10−5	1.4829 × 10−4	3.1145 × 10−4
	34	2.3089 × 10−5	3.0296 × 10−5	1.6676 × 10−4	2.9412 × 10−4
11th harmonic	12	3.4319 × 10−5	5.4946 × 10−5	2.9424 × 10−4	4.2399 × 10−4
	30	4.0944 × 10−5	5.7622 × 10−5	1.4622 × 10−4	3.1124 × 10−4
	34	2.6597 × 10−5	3.0091 × 10−5	1.5252 × 10−4	2.8334 × 10−4

.

The tables above show that, compared with other two ICA algorithms, the FKEO-ICA and the Fast-ICA algorithm can estimate harmonic currents more accurately. For all of the concerned frequencies, the correlation coefficients of the FKEO-ICA algorithm are closer to 1 compared with that of other ICA algorithms. Meanwhile, the error of the FKEO-ICA algorithm is relatively lower at all concerned frequencies, thus proving better accuracy of the FKEO-ICA algorithm compared with other ICA algorithms.

In order to test the running time of different ICA algorithms concerned with harmonic current estimation, simulations are conducted repeatedly 100 times. It can be shown in

Table 6. The running time statistics of different methods in 100 simulation tests.

Running Times	FKEO-ICA	Fast-ICA	ERICA	UNICA
Maximum value (s)	3.2792	3.0620	3.5602	2.6382
Minimum value (s)	3.0463	2.3236	2.4081	2.4381
Mean value (s)	3.1246	2.3714	2.5059	2.4913
Standard deviation	0.0474	0.0726	0.1141	0.0392

that the running time of all algorithms is around three seconds, while the proposed method can reach a higher estimation precision with a slightly longer time in comparison with other methods.

Based on the FKEO-ICA algorithm, the comparisons of normalized actual and estimated harmonic currents of bus 12, 20, and 34 at h = 5 are shown in Figure 5, Figure 6 and Figure 7, which indicate that the estimated harmonic currents are highly consistent with the actual harmonic currents.

To find the exact bus position of harmonic sources, the pair-wise conditional entropy $H({{\displaystyle \tilde{\mathit{I}}}}_{h\_est}|{\mathit{V}}_h)$ is calculated for each concerned frequency. The results are shown in Figure 8, Figure 9 and Figure 10, where the base of the algorithm in Equation (11) is assumed as 2.

As indicated in Figure 8, Figure 9 and Figure 10, although the orders of the estimated independent components (i.e., estimated harmonic currents) are random, the harmonic injection bus can be accurately identified when its voltage exhibits the minimum conditional entropy with the estimated harmonic currents.

7. Conclusions

This paper presents a method to identify harmonic sources based on FKEO-ICA and the minimum conditional entropy. In this study, the injected harmonic currents are estimated through the FKEO-ICA algorithm, in which the processed harmonic bus voltages are used as inputs. The advantage of the FKEO-ICA algorithm lies in that it does not require the prior information about system parameters. Additionally, the minimum conditional entropy is applied to locate the bus position of harmonic sources. The results indicate that the FKEO-ICA could achieve a better accuracy of harmonic current estimation, while the minimum conditional entropy can precisely locate the harmonic sources. In the simulations, only weak measurement noises are taken into consideration. In practice, the harmonic measurement quality and the detection technology can be further improved so as to avoid the increase of estimation errors caused by strong measurement noises.