Dynamic Impedance Estimation: Challenges and Considerations

: The objective of this paper was to examine the dynamic impedance estimation of electrical systems from online measurements. The paper makes several considerations and highlights the challenges to obtain a precise estimation. Transducer equalization and harmonic synchrophasor estimation (HSpE) are reviewed and discussed. The use of online and adaptive equalization for transducers proves to be a viable solution for improving voltage transducer’s (VT’s) and current transducer’s (CT’s) frequency response. Additionally, the use of oversampling algorithms can mitigate the effects of noise in the HSpE. Furthermore, methods for harmonic impedance estimation are discussed. The independent component analysis ICA-based dynamic impedance estimation is proposed and results presented, which yields excellent agreement. Finally, harmonic modeling and simulation of injected harmonic currents are used to observe resonances through the ampliﬁcation and attenuations and, consequently, the opportunity to conﬁrm the system self and transfer impedances of a test system. Dynamic impedance estimation will continue to be a great challenge for the power systems engineer as the system complexity increases with the massive insertion of power electronic inverters and the associated required ﬁltering. Real-time signal processing will be an effective tool to determine the dynamic self or transfer impedance.


Introduction
Estimating the magnitude and phase of harmonics in electrical power systems (EPS) still presents significant challenges, mainly when we focus on medium voltage (MV) and high voltage (HV) systems and in dynamic scenarios. There are still limitations for these systems, given the complexity and costs of developing measurement equipment that meets the voltage levels and frequency range determined by standards [1,2]. As a consequence, the harmonic measurements performed in the secondary of the transformers do not faithfully represent the distortions present in the primary, thus limiting the analysis of harmonic propagation in MV and HV and their applications, such as: the identification and quantification of the harmonic sources responsibilities, the estimation of the electrical system parameters, the correlation of errors in the protection system with the power quality (PQ), among others.
Methods for correcting the frequency response for voltage and current transducers (VT and CT, respectively) are discussed in the literature. Most of them are based on frequency domain correction [2,3] and require laboratory tests to obtain the frequency response and correction. This procedure makes it impossible to equalize the VT and CT already installed. Additionally, as reported in [1], the transducer's frequency response 2 of 18 depends on temperature and load impedance, among other factors, which often makes the calibration process in the laboratory unsatisfactory.
In addition to the errors caused by VT and CT, errors are arising from the harmonic component estimation algorithms. The estimation of harmonics is generally performed using the discrete Fourier transform (DFT), specifically through its fast algorithm, the fast Fourier transform (FFT). It is known that the estimation of the harmonic phasor by FFT is strongly affected if there is variation in the nominal frequency of the system and the noise level present in the signal is high.
Although the systems are designed for operating at 60 Hz (50 Hz), the imbalance between the generated power and the consumed power causes small variations in the frequency. However, even if less than 1%, these small variations can produce very inaccurate results in high-order harmonics estimations [4]. To reduce the error due to frequency deviation, the most usual solution is the use of interpolators or resamplers that guarantee synchronous sampling, reducing errors due to spectral spreading [4].
Resampling (interpolation) methods work very well when the signal's noise is low to the level of harmonics. However, when this does not occur, significant errors are introduced. A signal-to-noise ratio (SNR) of 40 dB (typical values in EPS range from 40-60 dB) degrades the interpolation efficiency producing total vector error (TVE) higher than 12% for high order harmonics.
It is important to note that resampling methods are extensively used to estimate the magnitude of harmonics. However, more recently, the harmonic phase's estimation has aroused great interest from the scientific community and companies, having seen the new applications made possible thanks to GPS synchronization [5,6]. The synchronized estimation of the harmonics' magnitude and phase can be named harmonic synchrophasor estimation (HSpE), as a direct extension of the concept of synchrophasor estimation. Consequently, the need for improving conventional algorithms to reduce the error in the phase estimation has been identified.
Methods to improve the estimation of harmonic components in noise are also available in the literature. Among them, the ESPRIT method stands out [7,8]. The ESPRIT is a parametric method and, therefore, needs to have a predetermined model order. Although it is a robust method, it has a high computational cost in relation to noise and does not need synchronous sampling, not being suitable for real-time implementation, especially in low-cost electronic systems. Thus, it is necessary to propose some strategies to reduce the noise's influence in estimating harmonic components through low-cost computational algorithms [9,10].
In considering possible applications of the harmonic synchrophasor, it is worth to mention the harmonic responsibilities attribution and the harmonic network estimation. The first application aims to quantify the responsibility for harmonic distortion, at a given point in the system, by the various harmonic sources [11] and the second allows, for example, performing islanding detection in systems with distributed generation [12], protection system adjustment, and improvement of the system control and operation.
Estimating a network's harmonic impedance based on online measurements is not an easy task. The methods that try to solve this problem are divided into passive and active methods. Passive methods seek to estimate the impedance from the current and voltage signals available on the measurement point without causing any disturbance in the network, although they can take advantage of any disturbances that affect the system. On the other hand, active methods cause some intentional disturbance in the network [13,14], and from the observation of current and voltage, they can estimate the impedance more accurately. Compared to the passive ones, the disadvantages of active methods are the possible deterioration in the PQ and external devices' need to cause disturbances in the network.
The injection of small signals called Gaussian modulated signals (GMS), as suggested in [15], does not degrade the PQ and produces accurate impedance estimates. Besides, this injection can be implemented in the control loop of the distributed generators' power converters. In this case, there is no need for external hardware to inject the GMS pulses. However, since most of the non-linear loads connected to the power system cannot inject the GMS to estimate the impedance, it is also important to investigate methods for estimating impedance by passive approaches.
In this way, estimation methods based on the independent component analysis (ICA) are highlighted [16][17][18]. These works partially solve the problem of estimating the impedance at the point of common coupling (PCC) on the load and utility sides. The challenge is to extend the methodology for estimating the harmonic transfers considering several harmonic sources contributing to a given point in the system.
Regardless of the method used to estimate harmonic impedance of the system, synchronized measurement is required, either of the harmonic phasor or the disturbance waveform (in the case of GMS injection), seeking to mitigate the effects of the fundamental frequency variation of the noise and the errors from instrumentation systems.
In this context, this work will discuss the challenges related to the dynamic impedance estimation and present some solution proposals. Figure 1 illustrates the three main points to be addressed in this work, which will be shown in Section 2, Section 3, and Section 4 of this paper. The conclusions of the work are then presented in Section 5.
impedance by passive approaches.
In this way, estimation methods based on the indepe are highlighted [16][17][18]. These works partially solve the p ance at the point of common coupling (PCC) on the load is to extend the methodology for estimating the harmon harmonic sources contributing to a given point in the sys Regardless of the method used to estimate harmoni chronized measurement is required, either of the harm waveform (in the case of GMS injection), seeking to mitiga frequency variation of the noise and the errors from instr In this context, this work will discuss the challenges ance estimation and present some solution proposals. Fi points to be addressed in this work, which will be shown paper. The conclusions of the work are then presented in

A. Transducer Equalization Review
With the increasing number of PQ monitors installe traditional instrument transformers are also being used fo ponents [1,2,19,20]. However, it is well known that traditio whether in distribution or transmission systems, presen response, both in magnitude and in phase, which introd harmonic components.
The process of correcting the frequency response of paper, will be called instrument equalization, using the

A. Transducer Equalization Review
With the increasing number of PQ monitors installed in the transmission networks, traditional instrument transformers are also being used for measuring the harmonic components [1,2,19,20]. However, it is well known that traditional measurement transformers, whether in distribution or transmission systems, present distortions in their frequency response, both in magnitude and in phase, which introduces measurement errors in the harmonic components.
The process of correcting the frequency response of measuring instruments, in this paper, will be called instrument equalization, using the terminology equivalent to the equalization of communication channels, whose purpose is exactly to correct the frequency distortions introduced by the communication channels.
Taking this into account, several methods for correcting these errors are further investigated. In [1], the authors investigate two possible approaches: the invasive and the non-invasive methods. The invasive method provides a more accurate correction; however, it requires disconnecting the instrument from the power grid and injecting test signals to obtain the frequency response, which is not always feasible. On the other hand, the noninvasive method is limited to harmonics already present in the investigated network that the instrument is connected to. Additionally, in the non-invasive method, it is necessary to temporarily connect a reference meter (typically a resistive divider) to equalize the already installed meters.
In [2], the authors present a laboratory calibration for three instrument transformers, one of which is used as a reference for measuring the frequency response. The work is based on a two-step procedure that uses isolated high voltage gas capacitors and a digital bridge, a digital current comparator. This procedure is adopted to explore the harmonics' amplification due to the differentiation circuit to obtain accurate measurements of high order and low magnitude harmonics.
The CTs are covered in more detail in [19]. The authors use the frequency coupling matrix technique, widely used in studies of harmonic penetration in power systems. The experimental results use the compensation technique in commercial CTs up to 10 kHz. The CTs are then characterized, modeled, and tested in the laboratory. Thus, it is possible to analyze their non-linearity, the phase dependence between harmonics, and propose compensation that allows the CT to maintain its accuracy class for a wide frequency range.
A problem scarcely addressed in the literature is CT's and VT's online and adaptive equalization. Many works, such as those mentioned above, work with offline techniques. In [3], the authors present an online compensation method, although it is not adaptive. The compensation is performed in real-time, but it is based on the transducer model's laboratory identification. Once the infinite impulse response (IIR) equalization filter has been identified, the reverse filter's coefficients are incorporated into the measurement system and are not updated anymore. Thus, if any characteristic of the transducer changes, the equalization will not be significant.

Adaptive Transducers Equalization
Online equalization of transducers in an adaptive way is a viable strategy but not yet explored in the literature. A possible solution suggested in this work is illustrated in Figure 2. In [2], the authors present a laboratory calibration for three instrument transformers, one of which is used as a reference for measuring the frequency response. The work is based on a two-step procedure that uses isolated high voltage gas capacitors and a digital bridge, a digital current comparator. This procedure is adopted to explore the harmonics' amplification due to the differentiation circuit to obtain accurate measurements of high order and low magnitude harmonics.
The CTs are covered in more detail in [19]. The authors use the frequency coupling matrix technique, widely used in studies of harmonic penetration in power systems. The experimental results use the compensation technique in commercial CTs up to 10 kHz. The CTs are then characterized, modeled, and tested in the laboratory. Thus, it is possible to analyze their non-linearity, the phase dependence between harmonics, and propose compensation that allows the CT to maintain its accuracy class for a wide frequency range.
A problem scarcely addressed in the literature is CT's and VT's online and adaptive equalization. Many works, such as those mentioned above, work with offline techniques. In [3], the authors present an online compensation method, although it is not adaptive. The compensation is performed in real-time, but it is based on the transducer model's laboratory identification. Once the infinite impulse response (IIR) equalization filter has been identified, the reverse filter's coefficients are incorporated into the measurement system and are not updated anymore. Thus, if any characteristic of the transducer changes, the equalization will not be significant.

Adaptive Transducers Equalization
Online equalization of transducers in an adaptive way is a viable strategy but not yet explored in the literature. A possible solution suggested in this work is illustrated in Figure   In this methodology, a reference transducer Hp(z) is installed in substation p. This transducer has a frequency response close to the ideal for the frequency range of interest. For example, a resistive divider can be used as a reference transducer for the entire distribution network. A constant k can approximate the transfer function (TF) of this transducer.
The system transfer function between the p and q buses is considered known at this stage and is described in the z domain by Hp,q(z). However, the frequency response of the transducer Hq(z) is unknown. The objective is to introduce an adaptive equalizer, He(z), to equalize the measurement transformer connected to the bus q, as shown in Figure 2  In this methodology, a reference transducer H p (z) is installed in substation p. This transducer has a frequency response close to the ideal for the frequency range of interest. For example, a resistive divider can be used as a reference transducer for the entire distribution network. A constant k can approximate the transfer function (TF) of this transducer.
The system transfer function between the p and q buses is considered known at this stage and is described in the z domain by H p,q (z). However, the frequency response of the transducer H q (z) is unknown. The objective is to introduce an adaptive equalizer, H e (z), to equalize the measurement transformer connected to the bus q, as shown in Figure 2. This adaptive equalizer has as inputs the signal u[n], the signal d[n] (a scaled replica of the measured signal on bus p), and error e[n]. Note that as the bus p and q are geographically distant, the signal d[n] and u[n] must be synchronized by a Global Positioning System (GPS) system for the proper functioning of the adaptive algorithm. Additionally, d[n] must be sent to the transducer equalizer installed on the q bus, so the adaptive algorithms, such as the least mean squares (LMS) or recursive least squares (RLS) [21], can be used to identify the system equalizer. This procedure can be performed online throughout the day, ensuring adaptive equalization.

Illustrative Results
To illustrate the methodology's feasibility, a 15-bar system [22] was simulated in Simulink/SimPowerSys. A resistive divider was fixed at bus 2, and inductive VT was connected at bus 10. Figure 3 shows the frequency response of the inductive VT. Note that the transducer has a resonance in the 2.711 kHz component. Thus, the upper harmonics will have their magnitude and phase significantly altered by the transducer.
ci. 2021, 11, x FOR PEER REVIEW identify the system equalizer. This procedure can be performed day, ensuring adaptive equalization.

Illustrative Results
To illustrate the methodology's feasibility, a 15-bar system [E not found.] was simulated in Simulink/SimPowerSys. A resistive d 2, and inductive VT was connected at bus 10. Figure 3 shows the the inductive VT. Note that the transducer has a resonance in the Thus, the upper harmonics will have their magnitude and phase the transducer. where A is the rated voltage and is the fundamental frequency The transfer function between nodes p and q was previousl duced in the methodology presented in Figure 2. The samples of quired to be used in equalization. The adaptive RLS algorithm sho used to estimate the parameters of the equalizer filter. After the co the equalizer is applied, and the harmonic phasors are estimated sults obtained are then compared with those expected, as describe Algorithm 1 The switching time instants to update the paramete multiples of a prefixed basic minimum constant sampling interva Now, consider the voltage signal measured on the p-bus, where the reference transducer was connected, is described by: where A is the rated voltage and ω 0 is the fundamental frequency. The transfer function between nodes p and q was previously estimated and introduced in the methodology presented in Figure 2. The samples of d[n] and u[n] were acquired to be used in equalization. The adaptive RLS algorithm shown in Algorithm 1 was used to estimate the parameters of the equalizer filter. After the coefficients' convergence, the equalizer is applied, and the harmonic phasors are estimated using the FFT. The results obtained are then compared with those expected, as described in Equation (1). N size for FFT 10 For n = 1, 2, . . . . Do: No convergence obtained then go to 3.1 19 Convergence is obtained at n = n 0 : 20 Go to 3.1 Figure 4 shows the TVE (total vector error) before and after equalization. The IEEE standard C37.118.2-2011 [23] specifies that the maximum TVE should be 1%. Before equalization, the TVE (%) for the transducer (black lines) had high values for interest harmonics. For example, harmonic 34 was close to 53%, as can be seen in the graphic. After the equalization (red lines), it is possible to observe a significant reduction in the TVE for all harmonics, which confirms the transducer frequency response's improvement. In the zoom, it is possible to note that the errors after the equalization are lower than 0.001%. Compute FFT(Y) 22 Go to 3.1 Figure 4 shows the TVE (total vector error) before and after equalization standard C37.118. 2-2011 [2323] specifies that the maximum TVE should be equalization, the TVE (%) for the transducer (black lines) had high values for in monics. For example, harmonic 34 was close to 53%, as can be seen in the gra the equalization (red lines), it is possible to observe a significant reduction in t all harmonics, which confirms the transducer frequency response's improvem zoom, it is possible to note that the errors after the equalization are lower than

Where to Go
With the advent of synchronized measurements in power systems, appli ing the harmonic phasor tend to multiply. However, the CTs and VTs scattere

Where to Go
With the advent of synchronized measurements in power systems, applications using the harmonic phasor tend to multiply. However, the CTs and VTs scattered throughout the system do not have an adequate frequency response for these applications. The replacement of these transducers with others with a better response is not economically viable. Thus, online and adaptive equalization seems to be a promising alternative.
The approach shown in this section served only to open the range of issues to be investigated in the coming years, among which are: What is the influence of the estimation error of the transfer function of the H p,q (z) system on the design of the equalizer? This point is essential, as it is known that the transfer function obtained from design data can present significant variations from the real system. The ideal would be to estimate these transfer functions by online procedure, but this falls into the process of recursive error (tail-chasing); 2.
On which bus, or buses, of the system should the reference meters be installed? This question seems similar to another one already addressed in the literature, which concerns the points of the network for PQ meters allocation; 3.
In telecommunications systems, it is common to use the blind equalization strategy [23]. The question is whether it is possible to use the same strategy for equalizing transducers.

Review Harmonic Synchrophasor Estimation
Estimating the harmonic phasor synchronized by satellite started in the 1990s, with a CHART (continuous harmonic analysis in real-time) [24]. However, it was only after the consolidation of the fundamental phasor estimation techniques, with the phasor measurement unit (PMU), that research for estimating the synchronized harmonic phasors had a new impulse. The harmonic synchrophasor estimation (HSpE) involves techniques and algorithms to simultaneously estimate and synchronize the harmonics' magnitude and phase.
In the method described in [28], the authors propose using FFT with the Hanning window to estimate the harmonic phasors. In the frequency domain, a correction is applied based on the system's instantaneous frequency to improve the estimated phasors' accuracy. Performance is tested up to the 9th order harmonic component only, in a limited test scenario without quantifying frequency deviations and without the presence of noise.
In [26], a contribution is presented, using KF for the individual estimation of each frequency component. The use of the TT has also been proposed for the task of estimating harmonic phasors [27][28][29], named Taylor weighted least squares (TWLS). Signal contamination by noise in these studies is considered, with cases with a SNR of 60 dB. However, the algorithms have a high computational cost, making their application in low-cost devices difficult.
In [30], the WT is used as an alternative to HSpE for the power system at the distribution level. The method first proposes applying an algorithm known as sample value adjustment (SVA) [31]. After this step, the FFT of the signal is applied, thus identifying the main frequency components. Next, wavelet filters are designed centered on the frequencies with significant energy. The SVA algorithm has the function of adjusting the signal samples, providing synchronous sampling in the FFT window to reduce spectral leakage.
Considering a phasor estimation in scenarios with significant noise contamination, more robust techniques such as ESPRIT [32] and TWLS present better results. However, they are parametric techniques and present high computational cost, being useful only for offline analysis. Besides, choosing the model order is always a challenge for parametric techniques.
In [33,34], it is proposed the use of interpolated DFT in the frequency domain (ipDFT). However, its performance lacks when it is desired to estimate low amplitude components in noise (individual harmonic content lower than 5%). In [35], ipDFT is compared with the improved SVA (MSVA). The phasor estimation results for both methods show that the interpolation in the frequency domain presents worse results concerning the time domain's interpolation strategy.
Time domain's interpolation strategy, also named resampling algorithm, seems to be one of the best candidates for HSpE because of the balance between the estimation performance and the computational effort. Figure 5 illustrates the basic idea behind this approach. Note that the algorithm is based on DFT and, consequently, suffers the spectrum leakage problem that can be reduced by the resampling algorithm. That is the key point for a successful harmonic phasor estimation: to choose a good resampling algorithm. Different algorithms have been reported in the literature for this purpose, among then the SVA and their variations [31,35], Lagrange, and B-Spline algorithms [4]. For example, the MSVA algorithm generates the interpolated samples using the following difference equation [5], The time-varying coefficients are a function of  that is a function of the actual frequency f (which varies with time) and the time instant n. This parameter is given by, Different algorithms have been reported in the literature for this purpose, among then the SVA and their variations [31,35], Lagrange, and B-Spline algorithms [4]. For example, the MSVA algorithm generates the interpolated samples using the following difference equation [5], where x[n] is the interpolated sample, x 0 [n] is the crude sample, and D i , i = −2, −1, . . . , 2 are time-varying coefficients given by, The time-varying coefficients are a function of α that is a function of the actual frequency f (which varies with time) and the time instant n. This parameter is given by, Despite the good performance of time-domain interpolation shown in the literature, the noise interference in the resampling algorithm's quality has been almost neglected. However, there is a degradation of the estimation under adverse conditions from the analyzed methodologies, especially when the SNR is less than 60 dB. In transmission networks, SNRs greater than 50 dB are common, but 45 dB values can be observed in distribution net-works. Noise interference seems to be one of the limiting factors in estimating high-order harmonic phasors.
A strategy to mitigate the effect of noise in the resampling process is using the oversampling technique [36]. Figure 6 illustrates this technique. In this case, the signal in continuous time is sampled with L times greater than the necessary. After sampling, an average moving filter (MAF) reduces the noise level, and the signal is then decimated to the required sampling frequency. It is noteworthy that some analog-to-digital converters (ADC) have this methodology as an internal function, which relieves the processor.  Figure 7 shows the performance of the FFT with a 1 cycle frequency of 15,360 kHz for HSpE, considering a signal contam and harmonics, from the 3rd to the 49th components with an am the fundamental component and variations at the system freque TVE graph is shown for each harmonic for oversampling values L = 4, L = 16, and L = 64. Each of the circles shows different offestimate the harmonics, a 3rd order Lagrange interpolator [37] w FFT. Note that the TVE error is around 6% for L = 1, reducing shows that noise reduction is an important pre-processing in the tion process.  To estimate the harmonics, a 3rd order Lagrange interpolator [37] was used, followed by the FFT. Note that the TVE error is around 6% for L = 1, reducing to 1.5% when L = 64. It shows that noise reduction is an important pre-processing in the harmonic phasor estimation process.  To estimate the harmonics, a 3rd order Lagrange interpolator [37] was used, followed by the FFT. Note that the TVE error is around 6% for L = 1, reducing to 1.5% when L = 64. It shows that noise reduction is an important pre-processing in the harmonic phasor estimation process.    MSVA, Lagrange ("Lagran" in the plot), and B-Spline ("Bspline" in the plot). The SNR was considered 60 dB, and 6000 cases were generated for each off-nominal frequency varying from 55 to 66 Hz (1 Hz of step-changing). Additionally, the phase was sorted from 0 to 2π, using a uniform distribution, and the magnitude of all harmonic, up to the 50th harmonic, was kept constant and equal to 5% off the fundamental component. Then the average and the maximum TVE were obtained. As shown in Figure 8, the B-spline interpolation algorithm presented the best results, followed by the MSVA. The maximum TVE was lower than 2% for the B-Spline algorithm and reached more than 3% for the two other algorithms. Only the B-Spline algorithm complied with the PMU standard for the fundamental frequency. Although there is no standard for the maximum TVE error for the harmonic phasor estimation, the results show a challenge to move the error to below than 1% for all harmonic components.

Where to Go
In order to obtain reliable harmonic phasor estimations, it is necessary to advance the research on the following points:  Develop resampling algorithms to ensure synchronous sampling that is robust to additive noise or to evaluate noise reduction techniques in pre-processing;  Investigate the effect of the delay of the frequency estimation on the interpolation performance over time. Since the frequency value used in the resampling process is always delayed about the current frequency value;  Define standards for HSpE, such as the TVE error allowed, the dynamic test conditions, the report's rate, etc.

Dynamic Impedance Estimation
The online estimation of the network impedance is a fundamental step in the study of harmonic propagation, attribution of harmonic responsibility, identification of losses in the network, among others [17,[38][39][40][41][42][43][44][45][46]. In recent years, research has focused on estimating impedance at the PCC, both on the load side and the source side [38,39]. However, more information would be obtained if, instead of knowing the PCC's impedances, it was possible to find the system transfer functions. For example, Figure 9 shows a small section of a real transmission system, with harmonic sources in buses A and C. Knowing the transfer from A to D and from C to D, it would be possible to understand the harmonic propagation in the network, as well as to establish the responsibilities of each source. As shown in Figure 8, the B-spline interpolation algorithm presented the best results, followed by the MSVA. The maximum TVE was lower than 2% for the B-Spline algorithm and reached more than 3% for the two other algorithms. Only the B-Spline algorithm complied with the PMU standard for the fundamental frequency. Although there is no standard for the maximum TVE error for the harmonic phasor estimation, the results show a challenge to move the error to below than 1% for all harmonic components.

Where to Go
In order to obtain reliable harmonic phasor estimations, it is necessary to advance the research on the following points:

•
Develop resampling algorithms to ensure synchronous sampling that is robust to additive noise or to evaluate noise reduction techniques in pre-processing; • Investigate the effect of the delay of the frequency estimation on the interpolation performance over time. Since the frequency value used in the resampling process is always delayed about the current frequency value; • Define standards for HSpE, such as the TVE error allowed, the dynamic test conditions, the report's rate, etc.

Dynamic Impedance Estimation
The online estimation of the network impedance is a fundamental step in the study of harmonic propagation, attribution of harmonic responsibility, identification of losses in the network, among others [17,[38][39][40][41][42][43][44][45][46]. In recent years, research has focused on estimating impedance at the PCC, both on the load side and the source side [38,39]. However, more information would be obtained if, instead of knowing the PCC's impedances, it was possible to find the system transfer functions. For example, Figure 9 shows a small section of a real transmission system, with harmonic sources in buses A and C. Knowing the transfer from A to D and from C to D, it would be possible to understand the harmonic propagation in the network, as well as to establish the responsibilities of each source.
While in the Norton model in Figure 10, there is no need for satellite synchronization of measurements, since all measurements are taken from the same geographical point, the transfer functions' estimates will require synchronized measurements.
Appl. Sci. 2021, 11, x FOR PEER REVIEW Figure 9. A small section of a real transmission system.
While in the Norton model in Figure 10, there is no need for satellite s of measurements, since all measurements are taken from the same geograp transfer functions' estimates will require synchronized measurements.

Impedance Estimation at PCC
The impedances in the Norton model of Figure 10 may be estima [11,17,18,41]. The main reason for focusing on this technique is that it is a p nique to be extended to the more general estimation of the transfer functio The ICA's basic concept for impedance estimation in electrical syste measurements made in the electrical system, identifying the individual har and the mixing matrix, which approximates the network's harmonic transfe tion (5) represents the generic form of ICA. In this equation, is the meas composed by currents and voltages at various points in the network, is th harmonic sources present in the system, but generally unknown, and is trix.
To estimate impedances seen from the PCC, the measurement vector the current and voltage at PCC, taken a long time. The ICA technique then signals' statistical nature to determine both the individual sources and the Then the ICA identifies the inverse matrix of , called decoupling matrix, [11]  s Hx   While in the Norton model in Figure 10, there is no need for satellite sy of measurements, since all measurements are taken from the same geograph transfer functions' estimates will require synchronized measurements.

Impedance Estimation at PCC
The impedances in the Norton model of Figure 10 may be estimate [11,17,18,41]. The main reason for focusing on this technique is that it is a pr nique to be extended to the more general estimation of the transfer function The ICA's basic concept for impedance estimation in electrical system measurements made in the electrical system, identifying the individual harm and the mixing matrix, which approximates the network's harmonic transfer tion (5) represents the generic form of ICA. In this equation, is the measu composed by currents and voltages at various points in the network, is the harmonic sources present in the system, but generally unknown, and is th trix.
To estimate impedances seen from the PCC, the measurement vector is the current and voltage at PCC, taken a long time. The ICA technique then e signals' statistical nature to determine both the individual sources and the m Then the ICA identifies the inverse matrix of , called decoupling matrix, a [11]  s Hx  where, is the vector of the estimated sources.
As known from the literature, the ICA has scale indetermination, i.e., th

Impedance Estimation at PCC
The impedances in the Norton model of Figure 10 may be estimated using ICA [11,17,18,41]. The main reason for focusing on this technique is that it is a promising technique to be extended to the more general estimation of the transfer function.
The ICA's basic concept for impedance estimation in electrical systems consists of measurements made in the electrical system, identifying the individual harmonic sources, and the mixing matrix, which approximates the network's harmonic transfer matrix. Equation (5) represents the generic form of ICA. In this equation, x is the measurement vector composed by currents and voltages at various points in the network, s is the vector of the harmonic sources present in the system, but generally unknown, and A is the mixing matrix.
To estimate impedances seen from the PCC, the measurement vector is composed of the current and voltage at PCC, taken a long time. The ICA technique then explores these signals' statistical nature to determine both the individual sources and the mixing matrix. Then the ICA identifies the inverse matrix of A, called decoupling matrix, as described in [11] where,ŝ = Î sÎL T is the vector of the estimated sources.
As known from the literature, the ICA has scale indetermination, i.e., the sources are estimated with a scale factor, Appl. Sci. 2021, 11, 558 12 of 18 Taking (7) in (6) results, To overcome the scaling indetermination given in (8), consider the equation of the Norton equivalent circuit shown in Figure 10, Taking the inverse of (9) and comparing with (8), one arrives at the final Norton impedance estimation, without scale indetermination:

Illustrative Results
For illustrative purposes, the ICA technique was used to estimate the source and the load impedance of the IEEE 15 bus system (see Figure 11), with the addition of two capacitor banks on buses 1 and 3.
To overcome the scaling indetermination given in (8), consider the equation Norton equivalent circuit shown in Figure 10, Taking the inverse of (9) and comparing with (8), one arrives at the final Nor pedance estimation, without scale indetermination:

Illustrative Results
For illustrative purposes, the ICA technique was used to estimate the source load impedance of the IEEE 15 bus system (see Figure 11), with the addition of two itor banks on buses 1 and 3. As the ICA is a passive method, the harmonic impedance of order h will accurately estimated if this harmonic is present in the system. To assess the perfo of the estimation, the existence of all harmonic components up to 50th was assu window containing 2 min of each harmonic samples was acquired at the PCC a cessed by the ICA algorithm. Figure 12 shows a typical variation for the 5th ha magnitude. As the ICA is a passive method, the harmonic impedance of order h will only be accurately estimated if this harmonic is present in the system. To assess the performance of the estimation, the existence of all harmonic components up to 50th was assumed. A window containing 2 min of each harmonic samples was acquired at the PCC and processed by the ICA algorithm. Figure 12 shows a typical variation for the 5th harmonic magnitude.
Appl. Sci. 2021, 11, x FOR PEER REVIEW   The ICA technique has some limitations that must be considered for its app The first limitation is that the harmonic sources must be statistically indepen weakly correlate with each other. Furthermore, ICA is a passive estimation techniq consequently, it is not possible to estimate the impedance in the frequency with v energy. Finally, the impedance is assumed constant within the window used for mation, and if it changes during the estimation, the results are inconsistent.  Figure 13a shows the impedance's magnitude on the source side and Figure 13b on the load side. Note that the ICA method can estimate both impedances with good precision. The phase, although not shown, was also estimated with accuracy.
ppl. Sci. 2021, 11, x FOR PEER REVIEW  Figure 13a shows the impedance's magnitude on the source side and Figure 1 the load side. Note that the ICA method can estimate both impedances with good sion. The phase, although not shown, was also estimated with accuracy. The ICA technique has some limitations that must be considered for its applic The first limitation is that the harmonic sources must be statistically independe weakly correlate with each other. Furthermore, ICA is a passive estimation techniqu consequently, it is not possible to estimate the impedance in the frequency with ver energy. Finally, the impedance is assumed constant within the window used for th mation, and if it changes during the estimation, the results are inconsistent.

Active Impedance Estimation Method
As previously mentioned, passive methods will only estimate the impedance The ICA technique has some limitations that must be considered for its application. The first limitation is that the harmonic sources must be statistically independent or weakly correlate with each other. Furthermore, ICA is a passive estimation technique, and consequently, it is not possible to estimate the impedance in the frequency with very low energy. Finally, the impedance is assumed constant within the window used for the estimation, and if it changes during the estimation, the results are inconsistent.

Active Impedance Estimation Method
As previously mentioned, passive methods will only estimate the impedance value for harmonic h if this harmonic is present in the power signal. To overcome this limitation, some disturbance must be caused in the network in a controlled manner to excite the network at interest frequencies. These disturbances can, however, significantly affect the PQ.
A solution proposed in [41,42] is a sequential injection of small short-duration signals and low energy. In this way, it is possible to scan the network's frequency without significantly affecting the power quality. Another advantage of the injection method is that it is possible to find the harmonic transfers in the system, as long as GPS synchronizes the measurements. Note that depending on the application, one will need not the synchronized phasor but the synchronized waveform.
In the method described in [41], signals called Gaussian Modulated Signal (GMS) are injected at a given point in the network, and the response to this disturbance is captured at other points in the network, determining the harmonic transfer. The GMS is mathematically described by, where i d (t) represents the current injection of the GMS form, f 1 the fundamental frequency of the network, h is the harmonic order whose transfer wants to be estimated, k is a constant that defines the resolution in the frequency, and G controls the pulse energy. The challenge is the injection of the signal in the network. However, when there are electronic converters, it is relatively simple to insert the pulse through their control.

Illustrative Results
Using this technique, the current transfer from A to D and C to D (Figure 9) was found. As shown in Figure 14 a,b, the present transfers values higher than the 1 for several frequencies. It means that harmonic sources injecting current at bus A or C can be amplified at bus D. It is also possible to observe the behavior of distributed parameters transmission line model by the resonance points.
. Sci. 2021, 11, x FOR PEER REVIEW 14 of A solution proposed in [41,42] is a sequential injection of small short-duration signa and low energy. In this way, it is possible to scan the network's frequency without sign icantly affecting the power quality. Another advantage of the injection method is that it possible to find the harmonic transfers in the system, as long as GPS synchronizes t measurements. Note that depending on the application, one will need not the synchr nized phasor but the synchronized waveform.
In the method described in [41], signals called Gaussian Modulated Signal (GMS) a injected at a given point in the network, and the response to this disturbance is captur at other points in the network, determining the harmonic transfer. The GMS is mathem ically described by, represents the current injection of the GMS form, the fundamental fr quency of the network, ℎ is the harmonic order whose transfer wants to be estimated, is a constant that defines the resolution in the frequency, and controls the pulse energ The challenge is the injection of the signal in the network. However, when there are ele tronic converters, it is relatively simple to insert the pulse through their control.

Illustrative Results
Using this technique, the current transfer from A to D and C to D (Figure 9) w found. As shown in Figure 14 a,b, the present transfers values higher than the 1 for sever frequencies. It means that harmonic sources injecting current at bus A or C can be amp fied at bus D. It is also possible to observe the behavior of distributed parameters tran mission line model by the resonance points.  Figure 15 shows the total impedance seen on bus A. Note that a current injection this node will cause a high voltage value due to the resonance frequencies.  Figure 15 shows the total impedance seen on bus A. Note that a current injection at this node will cause a high voltage value due to the resonance frequencies.

Harmonic Amplification
The use of simulation is recurrent in analyzing harmonic distortions in ele power systems. MATLAB/Simulink software has been increasingly used to analyz trical networks' performance and power quality. For the investigation of resonanc the contribution of multiple harmonics sources, the system shown in Figure 8 wa eled using the transmission system's typical values. So, 1A of three-phase curre injected for each harmonic frequency, up to 35th order (2.1 kHz), in the emitting b and C. The harmonic currents are injected with a fixed amplitude value and a phase angle.
To observe the distortions arising from these harmonic contributions, current urement is made on buses B and D. Figure 16 shows the current at buses B and D. Note that some harmonic currents are amplified at bus D. This is consistent w current transfer presented in Figure 16, especially for harmonic 7th and 13th.
Using online impedance or transfer information, it is possibly better to unde harmonic propagation behavior at the power system and identify harmonic respon

Harmonic Amplification
The use of simulation is recurrent in analyzing harmonic distortions in electrical power systems. MATLAB/Simulink software has been increasingly used to analyze electrical networks' performance and power quality. For the investigation of resonances and the contribution of multiple harmonics sources, the system shown in Figure 8 was modeled using the transmission system's typical values. So, 1A of three-phase current was injected for each harmonic frequency, up to 35th order (2.1 kHz), in the emitting buses A and C. The harmonic currents are injected with a fixed amplitude value and a typical phase angle.
To observe the distortions arising from these harmonic contributions, current measurement is made on buses B and D. Figure 16 shows the current at buses B and D.

Harmonic Amplification
The use of simulation is recurrent in analyzing harmonic distortions in electrica power systems. MATLAB/Simulink software has been increasingly used to analyze elec trical networks' performance and power quality. For the investigation of resonances and the contribution of multiple harmonics sources, the system shown in Figure 8 was mod eled using the transmission system's typical values. So, 1A of three-phase current wa injected for each harmonic frequency, up to 35th order (2.1 kHz), in the emitting buses A and C. The harmonic currents are injected with a fixed amplitude value and a typica phase angle.
To observe the distortions arising from these harmonic contributions, current meas urement is made on buses B and D. Figure 16 shows the current at buses B and D. Note that some harmonic currents are amplified at bus D. This is consistent with th current transfer presented in Figure 16, especially for harmonic 7th and 13th.
Using online impedance or transfer information, it is possibly better to understand harmonic propagation behavior at the power system and identify harmonic responsibility among other applications.

Conclusions
This paper addressed the main challenges in harmonic synchrophasor estimation (HSpE) and presented some paths for research evolution. The first aspect to be considered Note that some harmonic currents are amplified at bus D. This is consistent with the current transfer presented in Figure 16, especially for harmonic 7th and 13th.
Using online impedance or transfer information, it is possibly better to understand harmonic propagation behavior at the power system and identify harmonic responsibility among other applications.

Conclusions
This paper addressed the main challenges in harmonic synchrophasor estimation (HSpE) and presented some paths for research evolution. The first aspect to be considered when estimating the harmonic phasor concerns the voltage and current transducers' imperfections. The transducers present in the transmission substations or the distribution networks were designed to respond accurately to the fundamental component, presenting an inadequate frequency response for the harmonic range considered today in the PQ analysis (up to 50th harmonic). This paper presented the adaptive equalization method through algorithms already established in digital signal processing theory adapted to power networks as a proposed solution. For this, the installation of a reference transducer in a previously chosen point of the network and using a satellite synchronization system and availability of a data communication channel to carry out the equalization in a supervised and online manner is suggested. The initial results in this field are auspicious. The second aspect considered was related to the estimation algorithms of the harmonic phasor itself. The harmonic synchronized phasor measurements are still incipient, and further efforts are needed to improve the magnitude and phase estimation algorithms. One of the problems pointed out in this work concerns the error introduced by the resampling system in the magnitude and phase of the harmonic presence in the presence of noise. These errors can generate a TVE greater than 12% for a 40 dB SNR, especially if the harmonics present in the network have low energy, such as less than 5%. It was observed that these errors were introduced in the resampling circuit that is necessary to reduce the leakage error due to the asynchronous sampling rate. In this regard, the proposed solution was using an oversampling technique to reduce the noise level, although the search for new HSpE algorithms is also an exciting path. Finally, two possibilities of estimating the network impedance were presented, one passive using ICA and the other active through the injection of small pulses in the system with controlled frequency content. In particular, the importance of online estimation of network impedance and harmonic transfers was highlighted. The knowledge of these quantities in an online and dynamic way will make it possible to explain phenomena of harmonic propagation in the network and identify polluting agents over time and how each agent contributes to increase or reduce the total harmonic distortion of the network.
Dynamic impedance estimation will continue to be a great challenge for the power systems engineer as the system complexity increases with the massive insertion of power electronic inverters and the associated required filtering. Real-time signal processing will be an effective tool to determine the dynamic self or transfer impedance.