Study on the Estimation of Utility Harmonic Impedance Based on Minimum Norm of Impedance Difference

The existing methods for estimating harmonic impedance usually assume that harmonic impedance is constant in a measurement period. However, for a variety of practical reasons, such as the load change of power grid, utility harmonic impedance may be changed. So a reasonable hypothesis is that utility harmonic impedance is time-varying, but the variation of the impedance is not large at adjacent sampling time. Thus, this paper proposes a novel method to estimate utility harmonic impedance by taking the criterion of minimum norm of impedance differences at adjacent sampling times with a measurement period. According to the proposed criterion, an objective function is built with the norm of the harmonic impedance differences and utility harmonic voltage differences of adjacent moments. And then, the utility harmonic impedance can be obtained by minimizing the objective function. The proposed method has high accuracy even while background harmonic emission level is high and harmonic impedance of customer side is not much larger than that of utility side. The effectiveness is verified by the simulation and field cases. In the simulation, the estimation errors of most cases are within 10% and the maximum error is not exceed 16%.


Abbreviation List
With the widespread use of power electronic devices in power system, harmonic distortions have become one of the main concerns [1]- [5]. Because multi-infeed HVDC terminals, renewable energy converters, and electrical vehicle charging stations, etc., can inject a mass of harmonic into power system, which may cause the waveform distortions of The associate editor coordinating the review of this manuscript and approving it for publication was Ahmed F. Zobaa . voltage and current in the grid, and harms the stable operation of grid [6]- [8].
The accurate estimation of utility harmonic impedance is a prerequisite of assessing harmonic emission levels, quantifying harmonic responsibilities [9], [10], designing filters, verifying harmonic limit compliance and predication of system resonance [11]. However, due to the direct immeasurability of utility harmonic impedance, the accurate estimation of utility harmonic impedance is significant for suppressing the impact of harmonics and promoting the reward and punishment mechanism of harmonic management.
In practice, only the measured data of harmonic voltage and current at the PCC can be obtained, the utility harmonic impedance and utility harmonic voltage are unknown. Therefore, the measurement equations of the PCC are underdetermined at each measurement moment, and have no analytic solution for the utility harmonic impedance. Thus, the utility harmonic impedance can only be solved by the statistical methods.
In this regard, a variety of methods for utility harmonic impedance estimation have been developed, which can be divided into two categories: the invasive methods and noninvasive methods [12], [13]. Since the invasive methods are uneconomic and may influence the normal operation of power grid, the noninvasive methods are more common-used [14].
The existing noninvasive methods include: the fluctuation methods [15]- [17], the regression methods [18]- [22], the method based on covariance characteristic of random vectors [23], ICA methods [11], [24], [25], and other methods [26]- [30]. The fluctuation methods and regression methods can get accurate results with the preconditions, that is background harmonic is stable and customer harmonic impedance is much larger than that of the utility side. The method based on covariance characteristic of random vectors assumes that the weak correlation between harmonic current at a PCC and utility harmonic voltage, which can weaken the impact of background harmonic to a certain extent. When the harmonic impedances of both sides are approximately equivalent, the method will be invalid. ICA methods, a widely-used algorithm of blind source separation, are applied for estimating harmonic impedance of both sides. Based on the hypothesis that the harmonic sources of both sides are independent, the source signals can be separated from the mixing signal robustly. Then the harmonic impedance and contribution of each source can be calculated. However, ICA methods are susceptible to the correlation of harmonic sources and highly influenced by the size of sampling data. When the correlation of the harmonic sources of both sides is strong, the source signals are not independent and thus the separated signals may be inconsistent to the original sources. Although [11], [24] consider that the fluctuations of the harmonic sources of both sides are independent, it is still missing for proving the independence with rigorous mathematical derivation.
In addition, the aforementioned methods [15]- [30] are all based on this assumption that the harmonic impedance of utility side is invariable in a measurement period. However, in practice, the harmonic impedance of utility side is time-varying provided the loads of the grid change continuously, but the variation of the impedance is very small at adjacent sampling time. Based on this fact, a novel method of estimating utility harmonic impedance is proposed in this paper. An objective function is consisted of the norm of the utility harmonic impedance difference and the utility harmonic voltage. By minimizing the objective function, utility harmonic impedance can be obtained. The proposed method can get accurate results while background harmonic emission level is high and the harmonic impedances of both sides are approximately equivalent. The effectiveness of the proposed method is verified by simulation and field cases. This paper is organized as follows. Section II presents the harmonic circuit model and the rigorous mathematical derivation of the proposed method. Section III and IV utilize the simulation and actual data to verify the effectiveness of the proposed method, and the conclusions are drawn in section V.  Figure 1 shows the equivalent harmonic source circuit at a PCC. The circuit of utility-side is equivalent to Thevenin circuit and the circuit of customer-side is equivalent to Norton circuit. Equation can be established in accordance with the equivalent circuits in Figure 1:

II. CIRCUIT MODEL AND ALGORITHM PROCESS
and I pcc (n) = I pcc (n + 1) − I pcc (n), n ranges from 1 to N -1.
Based on the circuit equation in Eq. (1), utility harmonic impedance can be calculated by following expression: In practice, Z u is time varying, in order to perform harmonic impedance estimation at each point, the impedance difference of adjacent two sample points is introduced as Eq. (3): For the steady-state operation of power grid, in a very short time interval of adjacent sample points, Z u (n) is a small value. For the same reason, the background harmonic can be considered as stable and the fluctuation V u (n) is very small in the time interval. Thus, a linear combination of the norms of Z u (n) and V u is taken as an optimization goal. And then the objective function can be established as following: where || · || presents the 2-norm, and λ is the coefficient of the constraint, which indicates the extent of the influence of the 207390 VOLUME 8, 2020 background harmonic voltage. The selection of λ is discussed in later section. Note that the objective function is used for only one harmonic order in each analysis.

B. ALGORITHM PROCESS
Substitute Eq. (2) into Eq.(4), and Eq.(4) can be rewritten as: Rewrite the Eq. (5) into a vector form, we have: Expand Eq.(6) to calculate the optimal solution: where the symbol H represents conjugate transposition. Rewrite Eq. (9) and we have: According to the objective function in Eq. (10), is a constant value, so when the positive term is equal to zero, we can obtain the minimum value of the function. Thus, an estimated result of V u can be calculated as Eq. (11), and Z u can be calculated by substitutingV u into Eq. (2).V In general, the most impact factors for the estimation results are background harmonic and impedance magnitude ratio of both sides. In Eq.(4), Z u 2 is a very small value because the impedance difference value of adjacent sampling points is very small, the coefficient λ may affects the estimation results to some extent. Thus, a large number of simulations are performed to show the influence of λ on the estimation error. Note that the simulation parameters should be set to cover as many scenarios as possible.
The harmonic sources and impedances of both sides are set according to Figure 1 and consistent with the parameters in Section III. And parameter k indicates the harmonic current magnitude ratio of both sides at a PCC (i.e., k = |I u | |I c |, I u = V u /Z u ), which reflects the background harmonic emission level, m indicates the harmonic impedance magnitude ratio of both sides (m = |Z c | |Z u |). Let k = 0.5 and k = 1.2, and m = 2 and m = 10, respectively, and meanwhile, λ ranges from 10 −8 to 10 8 . Note that k = 0.5 and k = 1.2 correspond to the following scenarios, respectively, that: 1) the background harmonic emission level is low; 2) the background harmonic emission level is high. And m = 2 and m = 10 correspond with: 1) impedance values of both sides are approximately equivalent; 2) the customer impedance is larger than that of the utility side, respectively. Then the simulations of above scenarios are performed. The RMSPE of the estimated utility harmonic impedance given in Eq. (12) is used to assess the estimation effects.
where C real and C est are the real value and estimated value, respectively. The variation trends of the estimation error corresponding to different k, m and λ are shown in Figure 2.
From Figure 2 (a) and (b), regardless of the harmonic impedance ratio of both sides, it can be seen that the value of λ has little influence on the estimation errors. The estimation errors are dominated by k. Also, for Figure 2 (c) and (d), a similar conclusion can be drawn that the estimation errors are dominated by m and the estimation errors has little to do with the selection of λ. Thus, later in this paper λ is set to 1 for solving utility harmonic impedance.

D. ALGORITHM GUIDELINE
The algorithm process of the proposed method is summarized as follows.
1) Obtain the measurement data of harmonic voltage and current of a certain harmonic order at a PCC. 2) Divide the measured data into appropriate subintervals.
Usually the sample size of the subinterval is less than 200.

III. SIMULATION
In this section, the simulations based on harmonic source equivalent model in Figure 1 are performed to examine the effectiveness of the proposed method.    The computer flow chart for the simulated work is as follow: Result errors of the following methods are compared, they are: 1) regression method, 2) method based on covariance characteristic of random vectors, 3) complex ICA, 4) the proposed method. Harmonic currents and impedances are set as follows.
(1) Utility harmonic impedance is set to 15 + j20 , while the amplitude of customer harmonic impedance is set to m times that of utility side. Both amplitudes and angles of the impedances have 10% random fluctuations during the measurement period.
(2) For harmonic currents, the current source of customer side is initialized to 100A with an initial phase angle of 30 degrees, and both amplitude and phase angle have 30% random fluctuations; The amplitude of current source in utility side is set to k times that of customer side, i.e., |I u | = k |I c | , V u = Z u I u , the initial phase angle of utility side is set to -30 degrees, and both amplitude and phase angle have 20% random fluctuations.
The value range of m and k has the important impact factors on harmonic impedance estimation, m reflects the magnitude of the utility side impedance and k reflects the background harmonic emission level. In the simulation, some typical scenarios in practical engineering are tested. In each scenario, 2000 samples of harmonic voltage and current at the PCC are generated. With 200 samples in each subinterval (i.e., N = 200), the above four methods are applied in the simulation to estimate utility harmonic impedance. To compare the estimation accuracy of each method, the RMSPE of amplitude and phase angle are calculated according to Eq. (12), and the results are shown in Table 1 to Table 4.
In Table 1, it manifests that with low background harmonic level (i.e., k = 0.5), the proposed method has lower errors comparing with other methods even m is equal to 1.5. According to Table 2, even the background harmonic emission level is significantly large (i.e., k = 1.4), the proposed method is still more accurate than other methods.
The results in Table 3 correspond to the situation that the harmonic impedance of customer side is not much larger than that of utility side. According to Table 3, the errors of method 1 and method 2 is markedly larger than the other two methods, which shows that the two methods are invalid when the harmonic impedance of two sides is close and background harmonic is significant. And the proposed method has more accurate results than ICA method under this situation. Similarly, when the impedance ratio of both sides gets large, according to Table 4, the proposed method is still valid and more effective than other methods.
Results in Table 1 to Table 4 show that method 1 is greatly influenced by background harmonic emission level. When background harmonic emission level gets high, the result accuracy of method 1 is susceptible. Method 2 reduces the impact of background harmonic to some extent. However, when harmonic impedances of both sides are close, the covariance between I pcc and V u might not be equal to zero, which leads to the significant increase of the estimation errors of method 2. ICA method reduces the impact of background harmonic and reconstructs the original signals of harmonic sources of both sides. When ICA algorithm is applied to harmonic impedance estimation, the mixing matrix is composed of harmonic impedances. Obviously, the estimated Z u for ICA method is constant, while the actual impedance has 10% random fluctuation in both amplitude and phase angle during the whole segment, which leads to a deviation for the estimated result. Moreover, the harmonic sources of both sides may not completely independent, so the separated signal of harmonic source may not exactly match the original signal.
For further verify the effectiveness of the proposed method, set impedance ratio m range from 1.5 to 21.5, and k range from 0.1 to 1.2, which can cover the practical engineers. Performing the simulations 100 times for each m and k, and the average values of estimation errors are calculated. Figure 4 shows the estimation results. In Figure 4, it can be seen that the estimation errors are almost less than 15% in all considered scenarios, which indicates the proposed method is valid.
In Figure 5, the estimation errors of 100 times of simulations are shown in the form of boxplot, which corresponds to two scenarios: (a) m = 8 and k = 0.5, (b) m = 4 and k = 1.2. In the figure, the red line represents the median of the errors, the red ''+'' indicates an outlier, the values of lower and upper sides of the black rectangle are respectively first and third quartile, the two black horizontal lines connected with the black rectangle are respectively non-outlying minimum and maximum values [31]. It is observed that the estimated errors of the two scenarios are concentrated on a low error level, which is almost less than 10%, and the distance between the first and third quartile are small and less than 5%, besides, there are few outliers in both boxplots (a) and (b), which proves the stability and accuracy of the proposed method.

IV. FIELD CASE VERIFICATIONS
In this section, the field measurement data of an arc furnace and DC terminal are utilized to verify the effectiveness of the proposed method. In case 1, the 3 rd harmonic exceeds VOLUME 8, 2020  the limits, and in case 2, 11 th and 13 th harmonics exceed the limits. Hence the 3 rd and 11 th harmonics are taken as the analyzed data.
A. CASE 1 In Figure 6, the waveforms of the 3 rd harmonic voltage and current at the PCC, where a 150kV busbar feeds a 100 MVA DC arc furnace, are shown. The system fault level at the measured site is about 7500 MVA. The sampling frequency of measured data is 6400Hz, and FFT algorithm is applied to acquire the 3 rd harmonic data. The total measurement period is 10 hours and the sampling points are 600 [25], the corresponding adjacent sampling time is 1 minute.
With 60 samples in each subinterval, the method of sliding is applied on the measured data (i.e., 1-60, 2-61. . . ). Estimation results are shown in Figure 7. According to Figure 7, the conclusion can be drawn that the amplitudes and angles of the estimated Z u by these methods are close, the reason is that the 3 rd harmonic filter is not installed for the arc furnace, which causes the 3 rd harmonic impedance of the customer side is much larger than that of the utility side. Compared with other methods, the estimated harmonic impedance amplitude and phase angle of the proposed method are more stable.
During the downtime of the arc furnace, the customer devices don't work, the voltage at the PCC is approximately equal to the background harmonic voltage. Accordingly, 95% probability value of the utility harmonic emission level can be calculated as reference value [10]. The estimated utility  harmonic emission levels are shown in Table 5. Apparently, the emission level of proposed method is closest to the reference value, which indicates the effectiveness of the proposed method.

B. CASE 2
This case involves a 500kV multi-infeed HVDC system of a city grid that contains 4 DC terminals, where 12-pulse converters are adopted in DC transmission system, which mainly injects 11 th and 13 th harmonic into the city grid with a large number of cables. Somewhere the 11 th and 13 th harmonics exceed the limits. The transmission power of the converter stations are 3843MW, 2074MW, 1836MW, 1177MW, respectively. The filters are installed on the converter station side, thus, the customer harmonic impedance are not much larger than that of the utility. Note that in the city grid with four converter stations, when a converter station is chosen as the customer side harmonic source, then the utility side contains other three converter stations and other nonlinearity loads.
The 11 th harmonic voltage and current measured at a 500kV bus of a converter station are illustrated in Figure 8. 20 sampling points per minute are obtained by the power quality monitoring equipment and the corresponding adjacent sampling time is 3s, the total data sample length is 3600.
The sampling data are divided into 30 subintervals, and then there are 120 samples in each subinterval, where the mean value of the harmonic impedance is calculated. Thus, the estimation results of the 11 th utility harmonic impedance are shown in Figure 9. According to Figure 9, method 1 and method 2 have great fluctuations in the time period. For method 3, the estimated phase angle of utility harmonic impedance also has some fluctuations. Both amplitude and  phase angle of the proposed method change slowly and comparatively centralized. The polar diagram in Figure 10 presents the distributions of the amplitudes and phase angles of these methods. It can be seen that the harmonic impedance distribution estimated by the proposed method is more concentrated in the sampling period.
Besides, a rough reference value of utility impedance can be obtained with the short circuit capacity of the system. The short circuit capacity of the system is 35428 MVA, thus reference value of 11 th harmonic impedance is 77.62 , the average harmonic impedances of the methods are 33.65 , 121.47 , 83.14 , 75.82 , respectively. The results of the ICA method and the proposed method are both close to the reference impedance. While the variations of amplitude and angle of Z u estimated by ICA method are relatively large. Conclusively, the proposed method has more accurate estimation results comparing with other methods.

V. CONCLUSION
The existing methods of estimating harmonic impedance of utility side are usually based on the precondition that the impedance is constant in a measurement period. In practice, the harmonic impedance of utility side is variable, although the impedance difference of adjacent sample points is very small. Therefore, this paper proposes a novel method to estimate utility harmonic impedance. First, an objective function is established with norms of the utility harmonic impedance and background harmonic voltage, and then, by minimizing the objective function, the utility harmonic impedance can be estimated. The proposed method is still valid when the background harmonic emission level is high and the impedances of both sides are close. Simulation and field data verify the effectiveness of the proposed method.