Independent Component Analysis-Based Harmonic Transfer Impedance Estimation for Networks with Multiple Harmonic Sources ^†

Abstract

This paper presents a novel methodology to estimate the harmonic transfer impedances in electric power systems with multiple harmonic sources (HSs). The purpose is to determine the responsibility of each HS for the total harmonic distortion at a specific bus within the system, addressing a critical issue in the power quality field. To achieve this objective, it is necessary to estimate not only the individual HS, but also the transfer impedances between each source and the bus under analysis (BUA). Most methods for solving this problem are based on proper network modeling or restrict variations in harmonic sources to a single source at a time. The proposed methodology has overcome this limitation. For this, synchronized current and voltage phasors are measured at the BUA. Once the measurements are gathered, the Independent Component Analysis (ICA) method is applied to estimate the Norton equivalent. The harmonic transfer impedance (HTI) is then determined using the information provided by the ICA. To enhance the accuracy of HTI estimation, three procedures are employed for data mining the parameters provided by ICA over time to generate a well-conditioned system. Once the HTI is satisfactorily determined, the individual harmonic contributions (IHCs), i.e., the harmonic responsibility, can be estimated accurately. The effectiveness and performance of the method are demonstrated based on computational simulations using distribution and transmission systems. Additionally, the methodology is validated with real data collected from a Brazilian transmission system monitored by synchronized power quality measurement units. Simulated results show that the Total Vector Error (TVE) is less than 0.4%, and for the field data test, the TVE is less than 2%.

1. Introduction

The ever increasing use of nonlinear loads and power electronic devices results in power quality disturbances that must be analyzed during power system operation. One of the most significant problems is the proliferation of harmonic currents that are associated with several negative consequences such as additional power losses, overheating in cables, the premature aging of capacitor banks, the malfunction of electronic equipment and protection devices [1,2,3]. Due to the importance of this topic, it is utterly important to identify the major harmonic sources and to determine their individual contributions and responsibility associated with power quality degradation [4,5,6].

The assessment of harmonic distortion responsibility is a difficult task due to a reduced number of sensors installed in power grids. Current studies have predominantly focused on single-point analysis with the objective of determining the harmonic contributions from utility and customer sides connected at the same Point of Common Coupling (PCC) [7,8,9]. Within this context, it is extremely important to estimate the corresponding Norton equivalent at the PCC. There are some procedures that can be used to accurately estimate the impedances values including intrusive [10,11] and non-intrusive methods [12].

Intrusive methods are proposed based on the analysis of time-varying voltage and current signals assuming intentional disturbances introduced into the electric network caused by switched capacitor banks or passive filters, for example. However, these techniques can potentially affect the Power Quality (PQ) of the power grid. Some examples of the intrusive approaches include harmonic and inter-harmonic current injections at the PCC [13], the introduction of subtle signals [11], or the controlled activation of electric devices such as switched capacitor banks [14].

On the other hand, non-intrusive methods consider inherent variations in voltage and current measurements to estimate the network harmonic impedance. However, a negative aspect of these techniques is their susceptibility to provide inaccurate results particularly when harmonic components with lower magnitude and energy spectrum density are collected from the grid. Among the non-invasive methods, there are different approaches, including the linear regression method [15], the fluctuation method [16], methods based on data correlation [17,18], methods based on random vectors [19], methodologies based on Independent Component Analysis (ICA) [12,20,21,22], and techniques using data mining [23,24].

Since harmonic contributions assessment represents a contemporary issue for power quality, the exploration of multi-point methods emerges due to the ever increasing proliferation of harmonic sources distributed along the electric systems [25]. In a typical scenario with multiple harmonic sources, a consumer may report severe harmonic distortion at a specific point in the system, known as the Bus Under Analysis (BUA). The power utility recognizes that there are several harmonic sources within the system suspected of causing the problem or partially contribute to the overall distortion [26]. The goal of multi-point techniques is to estimate the harmonic impact of each suspected harmonic source on the harmonic distortion at the specified point, making it an essential approach for power quality management in complex electric systems [27].

Several multi-point techniques can be found in the literature. For example, in [28], a method based on Harmonic State Estimation (HSE) is presented. However, HSE-based methods require the proper modeling of the power system at the frequency of interest, which represents a challenge for real-time applications once a frequency domain approach is often used. To overcome this problem, data-driven methods have been developed. In [29], a method based on multiple linear regression analysis was proposed, while in [19], the least squares correlation method was applied to determine the harmonic contribution of each source. In this case, it is necessary to find a dataset where only one of the harmonic sources varies while the others remain approximately constant, making the process challenging.

These challenging characteristics have been addressed by multiple researchers, and new methods have been developed to overcome this problem, such as multiple linear regression [30], the constrained recursive least squares algorithm [31], the Bayesian partial least squares method [32], and the widely linear complex partial least squares method [25]. Despite the satisfactory performance of these methods, they were developed considering accurate network models and precise impedance parameters values, which can be extremely difficult to be obtained in real scenarios due to the scarcity of measurement units installed along power grids.

Methods based on ICA for multiple harmonic sources have been successfully used to identify the individual contribution and estimate the transfer matrix, considering a different approach in which the network impedances can be calculated during the estimation process. Reference [33] presents a modified ICA, called sparse component analysis for solving undetermined cases when the number of harmonic sources are higher than the number of measurement units installed in the power grid. In [34], an improved version of Complex ICA (CICA) was used considering non-GPS synchronized measurements. Even using the ICA algorithm to determine the impedance transfer matrix, the appropriate choice of harmonic measurement data impacts the estimated results. It may compromise the estimation of the harmonic contributions yielding unsatisfactory results, being extremely important to apply additional strategies for guaranteeing the observability of the system considering the selected set of measurements.

In summary, the methods for multiple harmonic sources in the current literature have one of the following limitations: (i) they require an accurate mathematical model of the power system; (ii) they depend on precise impedance parameters; (iii) they impose variations in harmonic sources to a single source at a time; or (iv) they must deal with the high-order problem of ICA.

This paper is an expanded version of the paper presented in [35], where a novel methodology was introduced for estimating harmonic responsibility in systems with multiple harmonic sources, without relying on detailed harmonic models or the localization of all harmonic sources. The proposed approach begins by estimating the Norton model at each monitored PCC using a single-point ICA approach, leveraging synchronized harmonic phasor data collected from Harmonic Phasor Estimation Units (HPMUs) or power quality monitors installed at the PCCs. In this way, the proposed method uses second-order ICA, i.e., two sources, two measurements, and is based solely on synchronized harmonic phasors. The Norton model is estimated dynamically over time, allowing us to continuously track the harmonic current injection and other parameters of the Norton model. With this synchronized information, it is possible to calculate the Harmonic Transfer Impedance (HTI) using data mining techniques, resulting in a well-conditioned system. The main contributions of the present paper are highlighted as follows:

• The proposed method consists of two steps: first, estimating the Norton model at each monitored PCC using the ICA method, and second, estimating the HTI.
• The method does not depend on detailed harmonic models of the power system, ensuring accurate HTI and Individual Harmonic Contribution (IHC) estimation.
• HTI and IHC are dynamically estimated using synchronized harmonic phasors from a monitored system.
• If a harmonic source is not monitored, its contribution is included in the supplier’s harmonic responsibility.
• A data mining methodology is employed to improve estimation accuracy.

This paper is structured as follows, considering this first Introduction section. In Section 2, a review of ICA is provided, and the process of obtaining the Norton model using single-point ICA is presented. In Section 3, the proposed methodology for the analysis of harmonic contribution with multiple harmonic sources is presented, and in Section 4, the results of the proposed method are presented, compared, and discussed. Finally, conclusions are highlighted in Section 5.

2. Norton Model Estimation Based on ICA

As mentioned in the introduction, the proposed method follows a two-step approach. In the first step, the Norton equivalent model is estimated at each monitored PCC to determine the primary harmonic currents of the sources. This estimation is performed using Independent Component Analysis ICA.

This section presents a brief review about ICA and shows how to use it to estimate the Nothon equivalent impedances and current sources in a single-point problem.

2.1. Independent Component Analysis (ICA)

ICA is a powerful statistical signal processing technique designed to uncover latent variables or mixed source signals in situations where the details of their mixing are unknown [36]. The basic ICA mathematical model can be expressed as

$$\mathbf{x}\left(t\right)=\mathbf{As}\left(t\right)$$

(1)

where $\mathbf{x}\left(t\right)={[x_1\left(t\right),\dots ,x_M\left(t\right)]}^T$ is the vector of observed variables, $\mathbf{s}\left(t\right)={[s_1\left(t\right),\dots ,s_N\left(t\right)]}^T$ is the vector of original sources, and $\mathbf{A}$ is an unknown constant mixture matrix. Note that the observed variables and the desired sources are discrete time sequences that are mixed with a constant matrix $\mathbf{A}$. ICA can be used for the separation of complex valued signal being a valuable tool in the context of power systems, facilitating thorough analysis within the frequency domain.

To take advantage of ICA’s ability, certain conditions regarding the original source signals and the measured signals must be met:

(C1:)

The number of measured signals, denoted as M, must match or exceed the number of sources, represented as N;

(C2:)

Among the original source signals, only one can exhibit a Gaussian distribution at most;

(C3:)

The original sources need to exhibit statistical independence from each other.

The purpose of the ICA is to estimate the separation matrix $\mathbf{H}={\mathbf{A}}^{-1}$ and then estimate the sources as follows:

$$\widehat{\mathbf{s}}\left(t\right)=\mathbf{H}\mathbf{x}\left(t\right)$$

(2)

where $\widehat{\mathbf{s}}\left(t\right)$ is the vector containing the estimation of the original sources.

There is a computationally efficient and robust approach to solve the ICA problem, called the Fast Fixed-Point Algorithm (FastICA). The FastICA is known for its rapid convergence, making it significantly faster than many other ICA algorithms; also it is robust to variations in the data and can handle noise effectively. Finally, the algorithm allows the choice of different non-entropy functions and whitening strategies, which can be adjusted as needed for different applications. These advantages make FastICA a popular choice for blind source separation problems and signal analysis. Furthermore, as the measures used in the methodology are harmonic phasors, this work uses the complex version of FastICA [37], which manipulates complex variables.

2.2. Single-Point ICA

Figure 1 shows the Norton equivalent for load and supply sides, seen from the PCC, for a specific harmonic. The index of the harmonic order is omitted in this figure and in the following equations for the sake of readability. In this figure, the impedances $Z_S$ and $Z_L$ of the Norton model compose the elements of the mixing matrix $\mathbf{A}$ in the ICA model, the measuring variables are the current phasor $I_{PCC}\left(t\right)$ and the voltage phasor $V_{PCC}\left(t\right)$; the desired sources are the injected harmonic source from the utility side $I_S\left(t\right)$ and load (customer) side, $I_L\left(t\right)$. The variable t in the phasor representation emphasizes that the proposed methodology explores the time varying behavior of these variables for accurate estimation. Only the Norton impedances are assumed constant during the estimation process. This ICA model is known in the literature as single-point ICA [12].

The restrictions presented on (C1) to (C3) are quite plausible in single-point ICA model. Note from Figure 1 that there are two sources to be determined (${\dot{I}}_S\left(t\right)$ and ${\dot{I}}_L\left(t\right)$) and two measurement variables, the current (${\dot{I}}_{pcc}\left(t\right)$) and voltage (${\dot{V}}_{pcc}\left(t\right)$). Therefore, (C1) is verified. Since there are only two sources in the single-point ICA problem, it can be assumed that they are independent. As the source of the utility side is a combination of several sources, except the source of the customer side at the considered PCC, from the Central Limit Theorem, it may be modeled as Gaussian. On the other hand, the customer source can be assumed to be non-Gaussian. Then, statements (C2) and (C3) are very likely to occur.

It is well known that the ICA may not preserve either the magnitude or the order of the sources. The way to overcome these limitation and to perform the correct estimation of the Norton parameters is presented in [23].

2.3. FastICA

The FastICA algorithm was initially developed for real-valued signals. However, there is an extension for complex signals. This extension for complex values retains the property of the cubic global convergence of its real counterpart. Furthermore, unlike gradient-based algorithms, FastICA does not require a prior decision on any learning rate parameter, making it simpler to use.

In this paper, the FastICA algorithm described in [20,38] is used. Initially, a step of centering and whitening is performed via singular value decomposition of the collected phasors, aiming to reduce the complexity of optimization and improve the accuracy of the results. FastICA can be utilized through negentropy or kurtosis, as it is an algorithm based on maximizing the non-Gaussianity of the sources. In this work, the approach via negentropy will be addressed.

Negentropy is a measure of non-Gaussianity in which, given a set whose random variables have the same variance, the Gaussian distribution is the variable with the highest entropy. Therefore, it is suggested to minimize Gaussianity by minimizing the marginal entropies of the estimates. Given a random variable y, its negentropy can be obtained as follows:

$$J_{Neg}\left(y\right)=H\left(y_{\mathrm{gs}}\right)-H\left(y\right),$$

(3)

where $H(·)$ is the entropy, and $y_{\mathrm{gs}}$ corresponds to a Gaussian-distributed random variable with the same mean and variance as y. Furthermore, maximizing negentropy maximizes non-Gaussianity. Due to the mathematical complexity of its determination, approximations such as the non-polynomial as defined in (4) is used:

$$J_{Neg}\left(y\right)\approx \alpha {\left(E\left\{G\right(y\left)\right\}-E\left\{G\right(v\left)\right\}\right)}^2,$$

(4)

where $\alpha$ is a constant, v is a zero-mean, unit-variance Gaussian random variable, and $G(·)$ is a nonlinear and non-quadratic function. In this work, a cubic function is used.

$$G\left(y\right)=y^3$$

(5)

This function was used due to its simplicity and computational efficiency. Since FastICA is based on the above approximation, it is assumed that it is desirable to estimate only one component, or, equivalently, to adjust a single line of $\mathbf{W}$ (equivalent to ${\mathbf{w}}_i^T$). Thus, an estimate will correspond to $y_i={\mathbf{w}}_i^T\mathbf{x}$. The algorithm can be summarized as follows:

• Randomly initialize ${\mathbf{w}}_i$;
• Execute the iterations using

  $${\mathbf{w}}_i\leftarrow E\left\{\mathbf{x}{G}^{\prime }\left({\mathbf{w}}_i^T\mathbf{x}\right)\right\}-E\left\{G^{″}\left({\mathbf{w}}_i^T\mathbf{x}\right)\right\}{\mathbf{w}}_i;$$

  (6)
• Apply the normalization:

  $${\mathbf{w}}_i\leftarrow {\displaystyle \frac{{\mathbf{w}}_i}{\parallel {\mathbf{w}}_i\parallel }};$$

  (7)
• Repeat step 2 and step 3 until convergence is achieved.

3. Proposed Method for Harmonic Contribution Analysis for Multiple Sources

For estimating the harmonic distortion responsibilities from the utility and load side at the PCC, it is crucial to estimate the Norton equivalent, as presented in the previous section. When there are more than two Harmonic Sources (HSs) in the network, one faces the issue of multiple HSs. In this case, it is important to estimate the IHC of each HS at a specific point, referred as the BUA. To accurately estimate the IHC, the HTI from every point where the HS is connected to the BUA must first be estimated. The HTI relates the voltage at the BUA to the current injected by the individual HS. This section will explain the methodology to obtain the HTI and the corresponding IHC. Figure 2 illustrates a power system with multiple harmonic sources.

In situations like the one shown in Figure 2, it may be interesting to determine the contribution of each source on bus k, for harmonic h. This problem may be written according to (8).

$$V_k\left(t\right)=\sum _{i=1}^M{z}_{ki}{I}_{si}\left(t\right)+V_0$$

(8)

where $V_k\left(t\right)$ is the voltage harmonic of order h on bus k, $I_{si}\left(t\right)$ is the ith current source for harmonic h, $z_{ki}$ is the transfer impedance between the current source to the voltage at bus k, $V_0$ is the background harmonic, and M is the number of harmonic sources. The current sources were previously obtained by single-point ICA. At this point it is important to highlight that all measuring must be synchronized. The challenge, however, is to estimate the HTI.

Suppose the voltage on bus k and the individual current sources are known at time instants, $t=t_1,t_2,\dots t_G$. Also, suppose the network HTI remain constant in the interval $t\in \{t_1\dots t_G\}$. So, it is possible to write (8) as

$$\left[\begin{array}{c}{V}_k\left(t_1\right)\\ ⋮\\ V_k\left(t_G\right)\end{array}\right]=\left[\begin{array}{cccc}1& I_{s1}\left(t_1\right)& \dots & I_{sM}\left(t_1\right)\\ ⋮& ⋮& \ddots & ⋮\\ 1& I_{s1}\left(t_M\right)& \dots & I_{sM}\left(t_G\right)\end{array}\right]\left[\begin{array}{c}{V}_0\\ z_{k1}\\ ⋮\\ z_{kM}\end{array}\right]$$

(9)

In matrix form, Equation (9) can be rewritten as

$${\mathbf{V}}_{\mathbf{k}}={\mathbf{I}}_{\mathbf{s}}\times {\mathbf{z}}_{\mathbf{k}}$$

(10)

The elements of matrix ${\mathbf{I}}_{\mathbf{s}}$ are the primary harmonic sources, obtained using single-point ICA, for each source presented in the system and for each chosen time instant. In practice, the number of observations G is greater than the number of harmonic sources M. So, there is a need to scrutinize ${\mathbf{I}}_{\mathbf{s}}$ in order to find the best invertible matrix, named ${\mathbf{I}}_{\mathbf{ss}}$. The proposed methods used to scrutinize ${\mathbf{I}}_{\mathbf{s}}$ are described in Section 3.1. The HTI is calculated as

$${\mathbf{z}}_k={{\mathbf{I}}_{\mathbf{ss}}}^{-1}\times {\mathbf{V}}_{\mathbf{k}}$$

(11)

With the elements of the ${\mathbf{z}}_k$ vector calculated, it is possible to estimate the harmonic contributions generated by source n on the bus k according to

$$V_{kn}=z_{kn}\times I_{sn}$$

(12)

where $V_{kn}$ is the harmonic voltage, of order h, generated by the hth harmonic source on the bus k.

The IHC can be evaluated in relation to the projection $V_{kn,proj}$ on the total harmonic voltage $V_k$, as shown in (13). Figure 3 shows the relationship between the voltage vectors and the projection for calculating the IHC.

$$IHC(\%)={\displaystyle \frac{|V_{kn}|·cos({\varphi }_k-{\varphi }_{kn})}{|V_k|}}\times 100\%$$

(13)

where ${\varphi }_k$ and ${\varphi }_{kn}$ represent the phase angle of $V_k$ and $V_{kn}$, respectively.

It is worth highlighting that the transfer impedances must remain constant during the time interval $t\in \{t_1\dots t_G\}$ used to perform the transfer impedances and the $IHC(\%)$ estimation. In order to manage this, the Norton equivalent impedances, estimated by the ICA, are monitored along time. The general flowchart of the proposed technique is presented in Figure 4.

In each segment in which $IHC(\%)$ will be estimated, the methodology starts with acquiring ${\dot{V}}_{pcc}$ and ${\dot{I}}_{pcc}$ at each access point, to estimate the Norton equivalent parameters through ICA. With these parameters in hand, the Norton impedance is evaluated. If the value of the impedance has changed, the system is reinitialized, since the impedances should remain constant in the estimation interval. Otherwise, the ${\mathbf{I}}_{\mathbf{s}}$ matrix and ${\mathbf{V}}_{\mathbf{k}}$ vector are built with the concatenation of the measured voltage and estimated load currents for each time instant, until the number of lines of ${\mathbf{I}}_{\mathbf{s}}$ is equal to G. Since $G>M$, some strategy might be used to solve (11), by either conditioning the matrix or using its pseudo inverse. Some strategies to solve this issue will be discussed in the next subsection. Then, the transfer impedances are obtained, and the IHC is estimated. The process repeats over time, so it is possible to estimate the variations in IHC during the day, for example.

It is worth highlighting that the methodology takes into account load fluctuations, as it uses time-varying harmonics as a basis, so load fluctuation does not affect the stability of the algorithm. On the other hand, grid reconfiguration alters the Norton impedance, and the algorithm to estimate HTI requires that the Norton impedance remain constant during the application of Equation (11). Since grid reconfiguration is a sporadic event during the day, this fact will not affect algorithm stability.

3.1. Scrutinize ${\mathbf{I}}_{\mathbf{s}}$ Matrix

To solve Equation (10), it is essential to find an invertible ${\mathbf{I}}_{\mathbf{s}}$ matrix. In this paper, three approaches are proposed to overcome the ill-conditioning problem of matrix ${\mathbf{I}}_{\mathbf{s}}$. The three methods presented will be compared in terms of accuracy and computational effort, in order to help end users make the most appropriate choice for their application.

3.1.1. Conditioning Number Approach

Since, typically, the number of observations G for every current source during the period in which the network impedance is constant exceeds the number of sources M, there are $N_p$ possible square matrices ${\mathbf{I}}_{\mathbf{s}}$, $M\times M$. The value of $N_p$ is determined by the equation

$$N_p={\displaystyle \frac{G!}{M!(G-M)!}}$$

(14)

The selection of the appropriate matrix is based on system conditioning. The conditioning of a linear system reflects how errors propagate from input to output data. A system is considered ill conditioned when the solution of the linear system is sensitive to small changes in the coefficients.

The mathematical parameter that expresses the conditioning of a linear system is the conditioning factor $\kappa$. For two matrices ${\mathbf{I}}_{s1}$ and ${\mathbf{I}}_{s2}$, the better conditioned system is the one with the smaller conditioning factor. For example, the system defined by ${\mathbf{I}}_{s1}$ is better than ${\mathbf{I}}_{s2}$ if $\kappa \left({\mathbf{I}}_{s1}\right)<\kappa \left({\mathbf{I}}_{s2}\right)$. The conditioning factor $\kappa \left({\mathbf{I}}_s\right)$ may be estimated using the following relationship:

$$\kappa \left(I_s\right)={\displaystyle \frac{{\lambda }_{max}}{{\lambda }_{min}}}$$

(15)

where, ${\lambda }_{max}$ and ${\lambda }_{min}$ are, respectively, the maximum and minimum eigenvalues of ${\mathbf{I}}_s$.

Then, if there are $N_p$ possible matrices to solve (10), the better one is the one with $\kappa \left(I_s\right)$ close to the unity value. However, depending on the number of sources (M) and the measurements (G) the number of tests to find the better well-conditioned solution may be prohibitive. For example, if there are $M=7$ harmonic sources and we obtain $G=3\times M=21$ measurements, the possible combination given by (14) will be $N_p$ = 116,280, which represents a prohibitive number of tests for online operations. The strategy to be used in this work, when the number of combinations exceeds 50, will be to choose, by random draw of the matrices’ rows, among the G available rows, fifty $I_s$ matrices to be tested. The number of possible tests will depend on the application and the available hardware to estimate the harmonic contribution.

3.1.2. Least Squares Approach

An alternative to solving Equation (10) is the application of the Least Squares (LS) method. This method is useful for estimating an approximate solution for overdetermined systems of linear equations, that is, systems, described by matrix ${\mathbf{I}}_s$, that have more equations than unknowns ($G>M$), with $rank\left({\mathbf{I}}_s\right)=M$.

The LS solution is obtained using the pseudo-inverse matrix, and is given by

$$\widehat{{\mathbf{z}}_k}={\left({\mathbf{I}}_s^{†}{\mathbf{I}}_s\right)}^{-1}{\mathbf{I}}_s^{†}{\mathbf{V}}_k$$

(16)

where the operator † denotes an Hermitian operator (the conjugate transpose).

In Equation (16), all measurements are used to form the matrix ${\mathbf{I}}_s$ of dimensions $G\times M$. However, many of the included rows may have temporal dependence, since the harmonic sources may not vary significantly in a short time period. Therefore, the pseudo-inverse can produce an ill-conditioned solution due to the collinearity of the measurements [30]. To improve the conditioning of the overdetermined matrix ${\mathbf{I}}_s$, some type of data mining should be incorporated to the algorithm. The next subsection will describe a simple procedure to improve this solution.

3.1.3. Reduced Least Squares Approach

In this work, we introduce a vector analysis algorithm with the purpose of selecting the most pertinent data from an available set. A straightforward approach would be to implement a mechanism that selects a vector $\mathbf{x}$ only if its similarity score meets certain criteria. The chosen similarity function is expressed by Equation (17), where ${\mathbf{x}}_s$ and ${\mathbf{x}}_i$ are two possible rows to be included in matrix ${\mathbf{I}}_s$. If the similarity factor is higher than $\tau$, row ${\mathbf{x}}_i$ is not included in the matrix; else, it can be included. This is presented in (18) where the limit value $\tau$ is empirically established. In this equation, ${\mathbf{x}}_{\mathbf{s}}$ represents the stored vectors, and ${\mathbf{x}}_{\mathbf{i}}$ is the input vector for analysis.

$$cos({\mathbf{x}}_s,{\mathbf{x}}_i)={\displaystyle \frac{{\mathbf{x}}_s^T.{\mathbf{x}}_i}{\left|\right|{\mathbf{x}}_s\left|\right|.\left|\right|{\mathbf{x}}_i\left|\right|}}$$

(17)

$$\left\{\begin{array}{cc}\mathrm{if}cos({\mathbf{x}}_s,{\mathbf{x}}_i)\le \tau ,& {\mathbf{x}}_i\mathrm{is}\mathrm{included}\\ \mathrm{else},& {\mathbf{x}}_i\mathrm{is}\mathrm{not}\mathrm{included}\end{array}\right.$$

(18)

In this context, the similarity function has the purpose of searching for vectors that tend to be orthogonal with the reference vector [39].

4. Results

The proposed methodology was evaluated in three distinct scenarios, all of them processed offline using Matlab/Simulink software (version 2020a). The first one used a reference system simulated at Simulink; the second one used real transmission system topology and parameters to generate simulated data; finally, the third one used filed data, acquired in a real transmission system.

To compare the estimated transfer impedance with the reference values, the Total Vector Error (TVE) is used as defined by Equation (19):

$$TVE(\%)=\sqrt{{\displaystyle \frac{{(Z_R-{\widehat{Z}}_R)}^2+{(Z_{Im}-{\widehat{Z}}_{Im})}^2}{{(Z_R-Z_{Im})}^2}}}\times 100\%$$

(19)

where ${\widehat{Z}}_R$ and ${\widehat{Z}}_{Im}$ are the real and imaginary parts of the estimated impedance. Analogously, $Z_R$ and $Z_{Im}$ denote the corresponding real and imaginary parts of the reference value.

The TVE is a valuable metric in phasor measurements and impedance analysis, as it simultaneously takes into account the magnitude and phase of the complex measurements.

4.1. Reference System Results

A benchmark test system for evaluating methods of harmonic responsibility was proposed in [40]. Figure 5 presents a schematic of this system, which has been slightly modified for use in the present study. The original benchmark is a three-phase symmetrical system consisting of an equivalent 110 kV high-voltage (HV) network with a rated short-circuit power of 3200 MVA and an X/R ratio of 10. The system includes a 110 kV/20 kV (YNyn) transformer with a transformation ratio of 110/21. Three loads are supplied from the medium-voltage (MV) busbars via separate MV lines. The customers’ PCCs are located on the MV side of their respective MV-to-low-voltage (LV) transformers. Customer 1, connected at PCC1, has a three-phase thyristor rectifier as a source of harmonic distortion. Customer 2, connected at PCC2, also has a three-phase thyristor rectifier in addition to linear loads. Customer 3, connected at PCC3, only has linear loads in the original configuration.

For the purposes of this study, we modified the network as follows: the nonlinear loads were modeled using a current source equivalent in parallel with an impedance. Additionally, the linear load for Customer 3 was changed to a nonlinear load. These modifications allow for better control of the experiment, as the harmonic phasors need to exhibit time variability for the methodology to function properly. Furthermore, the addition of a new nonlinear load creates another branch for estimating harmonic impedance transfer. More details about the original power system can be found in [40].

To analyze the performance of the proposed method, three cases were studied. In the first case (Case A), the MV compensator is disconnected from the system, and buses MV, PCC1, PCC2, and PCC3 are monitored. In the second case (Case B), the same buses are monitored; however, the MV compensator is connected to the system, changing the system impedance. The third case (Case C) is the same as Case A except that PCC3 is not monitored. In all cases, transfer impedances and harmonic contributions are analyzed at the MV bus.

To simulate the transfer impedance at harmonic frequencies up to the 50th order, the assumption that the sources produce all these harmonics was used. Clearly, this is not a practical situation, but it was assumed here to show that the method is capable of estimating the transfer impedances for any harmonic component that is present in the network. The sources that inject harmonics into the system are HV-GRID, Load 1, Load 2, and Load 4. The harmonic variation profile in the sources in this work is represented by a slow variation, that represents the variations along the day, added with a fast variation, which represents the instantaneous variations. The fast variation was obtained by adding a random component with Laplacian distribution as suggested in [41].

4.1.1. Case A—MV Compensator Disconnected from the System

In the first case, the MV compensator (“COMP. MV” in Figure 5) is disconnected from the system. After four sets of linearly independent voltage and current values, the transfer impedance estimation was performed, and the results are shown in Figure 6. Each plot represents the transfer impedance from a specific source to the bus MV. The top plots show, respectively, from left to right, the transfer impedance from the HV GRID bus to MV bus and from Load 1 to the MV bus. The bottom plots show, respectively, from left to right, the transfer impedance from Load 2 to the MV bus and from Load 4 to the MV bus.

Simulations were performed using the three conditioning algorithms explained in Section 3 as well as the Elastic Net Regression Method, presented in [30]. For simplicity, the transfer impedance results shown in Figure 6 are only for the “Conditioning Number Approach”. Note that the estimated results are very close to the reference ones given by Simulink. Then, in order to compare the accuracy of each method, the TVE for each one is portrayed in Figure 7.

From Figure 7, it can be seen that all tested conditioning techniques were able to accurately estimate the transfer impedance values, with the largest TVE being lower than $0.4\%$. However, the “Conditioning Number Approach” presented the best results.

With the use of the estimated HTI and the estimated harmonic injection of individual sources, it is possible to estimate the individual harmonic contribution of each source on the MV bus using (13). The IHC results are shown in

Table 1. Individual harmonic contribution generated by each system source in Case A.

Harmonic Order	HV GRID IHC (%)	LOAD 1 IHC (%)	LOAD 2 IHC (%)	LOAD 4 IHC (%)
1st (50 Hz)	100	0	0	0
5th (250 Hz)	59.78	18.35	8.60	13.26
7th (350 Hz)	35.59	39.89	9.40	15.12
11th (550 Hz)	51.47	−26.59	18.75	56.37
13th (650 Hz)	25.33	−8.11	10.48	72.29

.

From Table 1, each line corresponds to the IHC for harmonic h on bus MV. Note from this table that Load 1 behaves as a harmonic sink for the 11th and 13th harmonics. This means that Load 1 causes a damping effect attenuating the harmonic propagation in the grid. It is important to mention that the proposed method allows for computing the IHC over time, for example, every 10 min. This makes it possible to dynamically update the HTI, estimate the harmonics injected by each harmonic source, and determine the IHC at the specified point.

The method proposed in this paper does not analyze a single point in the network with the objective of determining contributions from both utility and customer sides. Differently from this traditional approach, the proposed method analyzes the contributions from each load to the whole power system in a multi-point and systemic approach. In this case, note that, according to Table 1, the results can be interpreted as follows:

• For the fundamental frequency (50 Hz), 100% of the contribution is attributed to the main substation, indicating that this is the only energy source for the distribution network. The other buses do not contribute, in this case, since there is no distributed generation or other source of energy in the system;
• For the fifth harmonic order (250 Hz), there are harmonic currents being injected into the system by all the load buses of the network. In this case, all points of the network individually contribute to the harmonic distortion of the power system. In this case, 59.78% of the contribution is associated with the HV grid, indicating that it represents the major contribution to the harmonic distortion, while the other loads have a total contribution of 40.22%. Differently from other methods in the literature, the proposed method indicates not only the total contrition of the loads but also calculates the individual harmonic contrition from each harmonic source. In this case, 18.35%, 8.60%, and 13.26% are the individual contributions from Loads 1, 2, and 4, respectively;
• For the eleventh harmonic order, there are harmonic injections from HV grid, Load 2 and Load 4. However, in this specific case study, no harmonic current is injected by the Load 1 in the 11th harmonic order. In this case, Load 1 absorbs the harmonic currents attenuating the propagation of the harmonic currents throughout the grid. This is commonly called the damping effect of the loads, as discussed in reference [42]. Note that the proposed method is able to identify this important and practical case which can occur in real scenarios. A given linear load may be affected by other distributed harmonic sources in the distribution system even if it does not inject harmonic current in the grid. In this case, the linear load serves as a “filter” to the harmonic distortions in the network. In practice, it yields several problems to the consumer (even if it does not injects harmonic currents) such as overload, power losses, premature ageing of electric devices, temperature rising in cables damaging its equipments and machines.

4.1.2. Case B—MV Compensator Connected to the System

In this case, the MV compensator is connected to the system, changing the system impedance. As shown in the flowchart of the proposed technique in Figure 4, if there is an impedance variation, the transfer impedance calculation returns to the initial state. Similarly to the previous case, the technique stores four sets of linearly independent voltage and current values to estimate the transfer impedances, which are shown in Figure 8.

In analyzing the results of Figure 8, it can be seen that the proposed technique was able to adapt to the system impedance variation, and it was capable to correctly obtain the new HTI, very close to the expected values. The TVE values, are presented in Figure 9 and, as in the previous case, the “Condition Number ”approach performs better than the other condition techniques.

Similarly to the previous case, the IHC for some sources were estimated and are presented in

Table 2. Individual harmonic contribution generated by each system source in Case B.

Harmonic Order	HV GRID IHC (%)	LOAD 1 IHC (%)	LOAD 2 IHC (%)	LOAD 4 IHC (%)
1st (50 Hz)	100	0	0	0
5th (250 Hz)	67.29	14.89	7.00	10.81
7th (350 Hz)	55.31	27.57	6.57	10.54
11th (550 Hz)	73.50	−16.84	11.11	32.22
13th (650 Hz)	52.37	−5.47	6.96	46.14

. Again, Load 1 behaved as a sink for the 11th and 13th harmonics.

In order to help users to choose the better method to estimate the HIT,

Table 3. Results about the method used to scrutinize ${\mathbf{I}}_s$.

	Conditioning Number	Least Squares	Least Squares Reduced	Elastic Net
TVE (%)	0.0527	0.0942	0.1383	0.0742
Simulation Time (s)	3.2045	2.0585	2.4913	2.52775

summarizes the performance of the four methods used in this paper, in terms of accuracy and simulation time. As can be seen the Conditioning Number method is the most accurate; however, it has the highest computation effort. It is worth highlighting that the typical time interval for power quality parameter estimation is 10 min, so the simulation time required by all algorithms is shorter than that and may be implemented in an online application.

4.1.3. Case C—Partially Monitored System

This test is the same as Case A, except PCC3 is not monitored. For this situation, it is expected that the harmonic distortion from Load 4 is included in the supply’s harmonic background distortion. Similar to the previous case, the technique stores four sets of linearly independent voltage and current phasor values to estimate the HTIs, which are shown in Figure 10.

In analyzing the results of Figure 10, it can be seen that all HTIs are the same, except for the transfer impedance from HV GRID to MV Bus, which has changed to accommodate the harmonic contribution of the non-monitored bus. The bottom right figure shows the difference in HTI from HV GRID to MV Bus between the monitored and non-monitored cases. Note that the HTIs from Loads 1 and 2 to the MV bus have not changed.

Table 4. Individual harmonic contribution generated by each system source in Case C.

Harmonic Order	HV GRID IHC (%)	LOAD 1 IHC (%)	LOAD 2 IHC (%)	LOAD 4 IHC (%)
Case C
1st	100.00	0.00	0.00	-
5th	78.09	14.90	7.00	-
7th	65.85	27.57	6.57	-
11th	105.83	−16.84	11.10	-
13th	98.52	-5.48	6.96	-
Case B (repeated here for convenience)
1st	100.00	0.00	0.00	0.00
5th	67.29	14.89	7.00	10.81
7th	55.31	27.57	6.57	10.54
11th	73.50	−16.84	11.11	32.32
13th	52.37	−5.47	6.96	46.14

shows the IHC for Case C. The results for Case B are repeated in this table to facilitate the interpretation of the new results. From this table, it is possible to verify that the IHC for Load 4 was incorporated into the supply side when PCC3 was not monitored. This is a very interesting and important result, as it indicates that it is possible to extract the IHC from the supply side by monitoring suspected loads.

4.2. Real Transmission System

To evaluate the proposed technique using real data, a 230 kV transmission system in the northern region of Brazil was used [43]. This system was chosen because it supplies power to a nearby aluminum processing industry that represents harmonic source for the system. The system’s versatility demonstrates that the method can also be applied to a high-voltage system. The single-line diagram of this system is shown in Figure 11.

Firstly, a computational simulation of this system was carried out using Matlab/Simulink. The main parameters of this system are presented in Appendix A.

In Figure 12, the top plot shows the HTI from Bus 1 to Bus B, while the bottom plot shows the HTI from Bus 2 to Bus B. The TVE is presented in Figure 13. Note that the maximum error was close to 6% for the Least Squares approach and lower than 2% for the “Conditioning Number” approach, which again was the technique with the best performance. HTI 1-B exhibits resonance around the seventh harmonic, and HTI 2-B exhibits resonance around the fifth harmonic. This means that if Bus 2 generates a fifth harmonic current, the voltage harmonic distortion at Bus B will be amplified by this resonance. The same can be said for the seventh harmonic current generated in Bus 1.

4.3. Field Data Results

Commercial PMUs commonly installed in transmission systems can accurately estimate the fundamental phasor. However, synchronized harmonic phasor estimation, or harmonic PMUs (HPMUs), investigated since the late 20th century [44,45,46], remains a challenging topic [47,48]. Currently, HPMUs are uncommon in real power systems. The main challenges in achieving an accurate and low-cost HPMU include the following: low harmonic energies in the presence of noise making phase estimation difficult; small variations in nominal frequency significantly affecting the phase estimation of high-order terms; and the dynamic behavior of harmonic phasors requires more sophisticated and computationally complex digital algorithms. These challenges have sparked considerable interest among researchers [25,49,50].

Points A and B of the real power system presented in Figure 11 were monitored using a synchronized PQ monitoring system described in [43]. Besides estimating PQ parameters, the developed system has an oscillograph feature that can be triggered by several PQ disturbances (e.g., THD, sag, transient, etc.). When an event is detected, the oscillograph stores 100 cycles of the voltage and current signals, which are synchronously sampled using GPS. The sampling rate of the monitoring system is 256 samples per cycle. Some events simultaneously affect Buses A and B, so there are synchronized oscillographies available. From the acquired waveform, synchronized harmonic phasors are estimated and, consequently, the HTI as well.

The time interval for this analysis spanned from 8 April 2021 to 18 April 2021. During this period, 36 events affecting both buses were identified. For each event, just under 100 harmonic phasors were obtained offline after discarding the signal transient part.

The average of the individual voltage harmonics monitored at Bus B during this period is shown in Figure 14 for illustration purposes. Note from this figure the predominance of harmonic orders 3, 5 and 7. The high value of the total harmonic distortion on Bus B (about 3.3%) can be explained by the proximity of aluminum processing industries in the region. Only harmonic phasors with energy above 0.5% were used to estimate the HIT, i.e., the third, fifth, seventh, and ninth harmonics, as well as the fundamental phasor.

The current phasor measured at PCC (${\dot{I}}_{PCC}$) and the source current (${\dot{I}}_S$), estimated by the ICA, are shown in Figure 15, where it is possible to note that their shapes are slightly different, as well as their average value.

The single-point ICA technique was applied to determine the Norton equivalent for each of the 36 sets of phasors, for Buses A and B. Remember that each set has about 100 harmonic phasors. From each set of phasors, the method estimates the Norton equivalent impedances, as well as the primary harmonic phasors from the utility and customer sides. The estimated primary harmonic phasors have many samples as the number of measuring samples used to estimate the Norton impedances, so, there is the need to apply some technique to scrutinize them to build the matrix ${\mathbf{I}}_{ss}$, in Equation (11). Each set of harmonic phasors estimated by the ICA gives rise to one HTI; therefore, 36 HTIs matrices are generated.

Figure 16 shows the average value of the transfer impedance amplitudes for each of the aforementioned harmonics, and for the four algorithms. The HTIs obtained by Matlab/Simulink, which were used for comparison purposes, are also included in this plot. The impedance $Z_{AA}$ represents the self-transfer impedance of Bus A, that is, how much harmonic current generated by Bus A induces harmonic voltage in the same bus. The impedance $Z_{AB}$ represents the transfer impedance from Bus A to Bus B, that is, how much harmonic current generated by Bus A induces harmonic voltage in Bus B. The impedance $Z_{BA}$ represents the transfer impedance from Bus B to Bus A, that is, how much harmonic current generated by Bus B induces harmonic voltage in Bus A. The impedance $Z_{BB}$ represents the self-transfer impedance of Bus B, that is, how much harmonic current generated by Bus B induces harmonic voltage in the same bus. Note the mean values for the four methods considered are very close to the simulated one.

The plots in Figure 16 show the average value of the HTI for the 10-day interval. The questioning that can arise is regarding the variation in the system parameters along these 10 days. It is possible to obtain some insight for this questioning by observing the TVE presented in Figure 17.

The bar errors in these plots were generated considering the 36-phasor set spread along 10 days. For clarity, only the Conditioning Number (the best method) and the Elastic Net (the comparison method) methods are presented in the figure. The relatively small standard deviation in the TVE, specially to the methodology of the Conditioning Number (lower than 2%), shows that the network does not experience significant changes in the considered set of data. This behavior can be explained because the monitored portion of the power system corresponds to a part of a short terminal system that is subjected to small variations. Another explanation is due to the fact the 36-phasor set used signals are spread along the 10 days, so there are several time slots not considered in the analysis. To illustrate the time variation behavior of the harmonic impedances, Figure 18 shows the third and fifth harmonics, along the 10-day measurement period.

Finally, with the currents estimated at each bus via ICA and the transfer impedances found, it is possible to estimate the IHC at Buses A and B, as shown in Figure 3.

Table 5. Average individual harmonic contribution for the results of field-collected data.

	${\mathit{IHC}}_{\mathit{AA}}$ (%)	${\mathit{IHC}}_{\mathit{AB}}$ (%)	${\mathit{IHC}}_{\mathit{BA}}$ (%)	${\mathit{IHC}}_{\mathit{BB}}$ (%)
1st (60 Hz)	101.01	−1.01	100.98	−0.98
3th (180 Hz)	89.22	10.78	14.64	85.36
5th (300 Hz)	90.08	9.92	14.83	85.17
7th (420 Hz)	82.91	17.09	25.24	74.76
9th (540 Hz)	100.41	−0.41	15.25	84.75

presents the average results of individual harmonic distortion obtained from field data using the Conditioning Number metric. In the columns in the table named $IH{C}_{ij}$, $i,j\in A,B$ represent the percentage of harmonic distortion at bus i that comes from bus j. Thus, $IH{C}_{AA}$ represents the voltage harmonic distortion in percentage generated by Bus A concerning the total harmonic voltage at Bus A, $IH{C}_{AB}$ represents the harmonic contribution in percentage generated by Bus B concerning the total harmonic voltage at Bus A, $IH{C}_{BA}$ represents the harmonic contribution in percentage generated by Bus A concerning the total harmonic voltage at Bus B, and $IH{C}_{BB}$ represents the harmonic contribution in percentage generated by Bus B concerning the total harmonic voltage at Bus B.

The results present in Table 5 enable the evaluation of harmonic contributions between load buses using the proposed method. For example, these results can be interpreted as follows:

• For the fundamental frequency (60 Hz), let us first analyze the results associated with Bus A. In this case, approximately 100% of the power comes from this point ($IH{C}_{AA}=100\%$), as the entire network demands the power from the main source at Bus A. This is indicated by the negative value of $IH{C}_{AB}$, meaning that Bus B serves as one of the consumers of the power generated at Bus A. Now, in considering the results from Bus B, while $IH{C}_{BB}$ is negative since it is a load bus, $IH{C}_{BA}$ is approximately 100%, indicating that the power comes from Bus A, and it is consumed at Bus B.
• For the third harmonic order (180 Hz), for example, let us first analyze the results from Bus A: 89.22% of the harmonic distortion is attributed to the load connected at Bus A, while 10.78% represents the harmonic contribution from Bus B into Bus A. On the other hand, if Bus B is under analysis, it can be noted that 85.36% is the contribution of the load at Bus B, while the other remaining 14.64% is the percentage value associated with the contribution from Bus A.
• A different situation is observed for the ninth harmonic order (540 Hz). For Bus A, approximately 100% of the distortion is attributed to it, while the other buses of the network absorb the currents generated from Bus A. In contrast, for Bus B, 84.75% represents the contribution of the load connected at Bus B, while the remaining percentage value (15.25%) indicates the contribution from Bus A into the harmonic distortion at B.

5. Conclusions

This paper presents a new method for estimating the HTI and the IHC in an electrical power system with multiple harmonic sources. The method requires synchronized current and voltage phasors estimated at the bus under analysis and at every PCC where harmonic sources are expected. Initially, single-point ICA is applied at each monitored PCC to estimate the Norton equivalent model, thereby accessing the Norton impedances and the primary harmonic current injected from both the utility and customer sides. The Norton model is estimated at several time instants, generating a temporal sequence of time-varying harmonic phasors. Then, through mining a synchronized subset of time-varying harmonics at all monitored points, it is possible to estimate the HTI and the IHC.

The main advantages of the proposed method are the following: (i) it does not rely on harmonic models of the power system to accurately estimate the HTI and, consequently, the IHC; (ii) these parameters are dynamically estimated using synchronized phasors; (iii) if a harmonic source is not monitored, its harmonic contribution is incorporated into the supplier’s contribution; and (iv) a data mining methodology is employed to obtain a well-conditioned system, thereby enhancing estimation accuracy.

In the Results section, firstly, a benchmark test system for harmonic analysis was simulated to validate the methodology. The TVE for the HTI estimation was less than 0.4%. A real case of a transmission system was also used in this study. In this real system, the synchronized harmonic phasors were obtained by processing synchronized waveforms acquired by a monitoring system. With the phasors, the HTI was estimated with an error lower than 5%. The results presented demonstrate that the proposed methodologies are suitable for estimating HTI in electrical power systems (distribution or transmission systems), and consequently the individual harmonic contribution along time.

For a fair assignment of harmonic distortion responsibility, the continuous monitoring of synchronized harmonic phasors is required. However, this was not feasible in the current study due to the unavailability of such monitoring capabilities. Therefore, for future work, we plan to develop harmonic PMUs to enable continuous access to harmonic phasors and to design a statistical approach for evaluating the IHC over time.