# A Harmonic Impedance Estimation Method Based on the Cauchy Mixed Model

Authors: Zhirong Tang; Huaqiang Li; Fangwei Xu; Qin Shu; Yue Jiang
Journal: Mathematical Problems in Engineering
 (2020)

---

## Abstract

In this paper, a new method without any tradition assumption to estimate the utility harmonic impedance of a point of common coupling (PCC) is proposed. But, the existing estimation methods usually are built on some assumptions, such as, the background harmonic is stable and small, the harmonic impedance of the customer side is much larger than that of utility side, and the harmonic sources of both sides are independent. However these assumptions are unpractical to modern power grid, which causes very wrong estimation. The proposed method first uses a Cauchy Mixed Model (CMM) to express the Norton equivalent circuit of the PCC because we find that the CMM can right fit the statistical distribution of the measured harmonic data for any PCC, by testing and verifying massive measured harmonic data. Also, the parameters of the CMM are determined by the expectation maximization algorithm (EM), and then the utility harmonic impedance is estimated by means of the CMM’s parameters. Compared to the existing methods, the main advantages of our method are as follows: it can obtain the accurate estimation results, but it is no longer dependent of any assumption and is only related to the measured data distribution. Finally, the effectiveness of the proposed method is verified by simulation and field cases.


---

## Body

## 1. Introduction

With the continual increase of the capacity of renewable energy and power electronic equipments in power grid, the utility harmonic problems are more and more serious, and the public requirements for power quality are more and more raised, and thus, the harmonic control becomes more significant [1, 2]. The harmonic suppression requires accurately acquiring the harmonic impedance of the utility side [3–7].At present, the mainstream methods for estimating the utility harmonic impedance are noninvasive, which mainly include the linear regression method [8–13], fluctuation quantity (FQ) method [14, 15], independent random vector (IRV) method [16], and independent component analysis (ICA) method [17–22]. The linear regression method uses linear regression to calculate Norton equivalent circuit equation coefficient and obtain the harmonic impedance. The FQ method uses the sign characteristics of the ratio of harmonic voltage fluctuation to harmonic current fluctuation at a point of common coupling (PCC) to estimate the utility harmonic impedance. The linear regression method and FQ method can effectively estimate the utility harmonic impedance under the condition that the background harmonic is stable and small, and the customer harmonic impedance is much larger than that of utility side. The IRV method assumes the harmonic current of the PCC is weak-correlation with the utility harmonic voltage, i.e., assumes implicitly that the harmonic impedance of customer side is much larger than that of utility side, which means that the harmonic current of the PCC is dominated by the customer. Thus, the covariance of the harmonic current of the PCC and the utility harmonic voltage can be approximated to zero, and then by solving the covariance equation, the utility harmonic impedance can be obtained. The main characteristic of ICA is that it can simultaneously figure out the harmonic impedances of both sides. However, ICA requires that the harmonic sources of both sides are independent.Generally, the existing methods need one or one more assumptions, which includes that the background harmonic is stable and small, the harmonic impedance of the customer side is much larger than that of utility side, and the harmonic sources of both sides are independent. However, for modern power grid, these assumptions are hard to be held up. For instance, the harmonic sources caused by renewable energy and power electronic equipment in power grid widely distribute, the harmonic sources among power electronic equipments interact, and harmonics amplify for transmission lines and cables. Thus, it is obviously unreasonable that the background harmonic is assumed to be very small and stable [17]. Besides, filters or reactive compensation capacitors are usually installed on the customer side, which may cause the harmonic impedances of both sides to be approximately equivalent. As for, the independence hypothesis of the harmonic sources of both sides is much idealized, and many researchers have pointed out that the assumption may not meet the practice [4, 5, 17, 20, 22].In this paper, a novel method without any assumption is proposed, which employs a Cauchy Mixed Model (CMM) to fit the distribution characteristic of the measured data at the PCC because we have found that the statistical distribution of the harmonic measurements of the PCC can be depicted with CMM through the analysis of the field harmonic data. So, a fitting CMM will be chosen to match the distribution relationship between the harmonic voltage, the harmonic current, and the utility harmonic impedance of the PCC. Also, then the utility harmonic impedance can be obtained by the expectation maximization algorithm (EM).This paper mainly studies from the following four aspects.(1)
The principle of CMM is introduced, and the Norton equivalent circuit of the PCC is expressed as a CMM(2)
The feasibility of CMM is verified by field measured data(3)
The CMM estimation processes are derived in detail(4)
The effectiveness and robustness of the proposed method are verified by comparing the estimation results of simulation and field data with existing algorithmsThis paper is organized as follows. Section2 presents the analysis for the validity of CMM. Section 3 derives the process of estimating the utility impedance. Sections 4 and 5 verify the validity of the proposed method by simulation and field case. The paper is concluded in Section 6.


## 2. Basic Principles

### 2.1. Harmonic Equivalent Analysis Model

Norton equivalent circuit of the utility side and customer side is shown in Figure1.Figure 1
The Norton equivalent circuit.In Figure1, I˙s is the utility harmonic current, and Zs is the harmonic impedance on the utility side. U˙PCC and I˙PCC are the harmonic voltage and current measured at the PCC, respectively.According to the Norton equivalent circuit, the equation can be written as(1)ZsI˙s=U˙PCC+ZsI˙PCC.We have(2)ΔI˙sn=I˙sn+1−I˙sn,ΔI˙PCCn=I˙PCCn+1−I˙PCCn,ΔU˙PCCn=U˙PCCn+1−U˙PCCn,n=1,2,…,N−1,where, n is the location of the sampling point, and N is the number of sampling points. ΔI˙sn is the fluctuation of the utility harmonic current source, ΔI˙PCCn is the harmonic current fluctuation at the PCC, and ΔU˙PCCn is the harmonic voltage fluctuation at the PCC.Then, equation (1) can be expressed as(3)ZsΔI˙s=ΔU˙PCC+ZsΔI˙PCC,where(4)ΔI˙s=ΔI˙s1,ΔI˙s2,⋯,ΔI˙sN−1,ΔU˙PCC=ΔU˙PCC1,ΔU˙PCC2,⋯,ΔU˙PCCN−1,ΔI˙PCC=ΔI˙PCC1,ΔI˙PCC2,⋯,ΔI˙PCCN−1.Thus, equation (3) can be rewritten as(5)Zs−x−Zs−yZs−yZs−xΔIs−xΔIs−y=ΔUPCC−xΔUPCC−y+Zs−x−Zs−yZs−yZs−xΔIPCC−xΔIPCC−y,where ΔIs−x and ΔIs−y are the real and imaginary parts of the utility harmonic current, respectively. ΔUPCC−x and ΔUPCC−y are the real and imaginary parts of the harmonic voltage measured at the PCC, respectively. ΔIPCC−x and ΔIPCC−y are the real and imaginary parts of the harmonic current measured at the PCC, respectively. Zs−x and Zs−y are the real and imaginary of the harmonic impedance on the utility side, respectively.Therefore,(6)ΔUs=ΔUPCC+Z^sΔIPCC,where(7)ΔUs=ΔUsxΔUsyT=Zs−x−Zs−yZs−yZs−xΔIs−xΔIs−y,ΔUPCC=ΔUPCC−xΔUPCC−yT,Z^s=Zs−x−Zs−yZs−yZs−x,ΔIPCC=ΔIPCC−xΔIPCC−yT.


### 2.2. Cauchy Mixed Model

AD-dimensional Cauchy probability density function is(8)φxμ,∑=Γ1+D/2Γ1/2πD/2∑1/21+x−μT∑−1x−μ1+D/2,where μ is the data center point, ∑ is the symmetric semidefinite matrix, D is the data dimension, and Γ⋅ is the Gamma function.If given observation dataX=x1,x2,…,xN, the CMM is(9)PXα,μ,∑=∑k=1KαkφXμk,Σk,where K, α=α1,α2,⋯,αK, μ=μ1,μ2,…,μK, and ∑=∑1,∑2,…,∑K are the number of models, weighted values, data center points, and covariance matrices, respectively.


### 2.3. Feasibility Analysis of the Model

This section introduces the fitting methods of data distribution and the test of fitting results by CMM. First, we introduce the Copula function and Sklar’s theorem [23, 24]:Definition 1.
A copula is a functionC of D variables on the unit D-cube [0,1]D with the following properties:(1)
The range ofC is the unit interval [0,1](2)
C(x) is zero for all x = (x1, …, xD) ∈ [0,1]D for which at least one coordinate equals zero(3)
C(x) = xd if all coordinates of x are 1 except the dth one(4)
C is D-increasing in the sense that for every a ≤ b in [0,1]D, the measure ΔCba assigned by C to the D-box [a, b] = [a1, b1] × … × [aD, bD] is nonnegative, i.e.,(10)ΔCba=∑ε1,…,εD∈0,1−1ε1+…+εDCε1a1+1−ε1b1,…,εDaD+1−εDbD≥0.
The key theorem due to Sklar clarifies the relations of dependence and the copula of a distribution:Theorem 1.
LetF denote a D-dimensional distribution function with margins F1, …, FD. Then, there exists a D-copula C such that for all real x = (x1, …, xD),(11)Fx=CF1x1,…,FDxD.
Therefore, theC is unique, that is,(12)Cx=FF1−1x1,…,FD−1xD,where F1−1⋅,…,FD−1⋅ are quasi-inverses of the marginal distribution functions.
To verify whether the CMM can fit the probability distribution of the harmonic voltageU˙PCC and current I˙PCC data, we combine the two-dimensional t-copula function [25] with the Monte Carlo method to generate a set of random data conforming to CMM. Then, the Cumulative Distribution Functions (CDF) of field data and generated data are calculated and expressed as curves, respectively. Finally, we use the Kolmogorov–Smirnov (KS) test [26] to detect the maximum distance between two CDF curves.
The steps of generating the CMM data are as follows:(1)
K models are selected, and the CMM is adopted to PCC data. Then, the corresponding weightαk and covariance matrix ∑k are obtained.(2)
Convert the obtained covariance matrix∑k into a correlation matrix ρk.(3)
Thet-Copula function is used to generate random two-dimensional Cauchy number of K groups associated with ρk.(4)
Generate N random numbersni from 0 to maxαk. If ni<αk, select the random number set of the kth two-dimensional Cauchy model.
According to the abovementioned method, the CMM is used to fit the distribution of two groups’ field data. The first set of field data is conducted at the site where a 150 kV bus bar feeds a 100 MW Electric Arc Furnace (EAF). The total sample time is 10 hours, and the calculated 3rd harmonic sampling values are averaged per minute, then the number of the all samples is 600. CMM random numbers are generated corresponding to the EAF data. The CDF of the data generated by CMM and the EAF data are shown in Figure2.
The results of the KS test are shown in Table1.
In Table1, H = 0 indicates acceptance of the hypothesis at the significance level of 5%, P is the probability of acceptance, and KSSTAT is the value of the test statistic. It can be seen from Figure 2 and Table 1 that the CMM can fit EAF data.
The second set of field data is taken from a city power grid with multiple DC terminals, and the measured data are collected from 600 data points (20 sampled data per minute) at a 500 kV bus of the DC terminal for the 11th harmonic voltage and current. CMM random numbers are generated corresponding to the DC terminal data. The CDF of the data generated by CMM and DC terminal data are shown in Figure3.
The results of the KS test are shown in Table2.
The meanings of the parameters in Table2 are the same as those in Table 1. According to Figure 3 and Table 2, it is obvious that the CMM can fit DC terminal data even though there is a slight deviation between the CDF curve of the CMM and that of DC terminal.Figure 2
The CDF of the data generated by CMM and the EAF data. (a) Real part of current. (b) Imaginary part of current. (c) Real part of voltage. (d) Imaginary part of voltage.
(a)
(b)
(c)
(d)Table 1
KS test result of CDF of EAF.
Parameter
Current
Voltage
Real
Imaginary
Real
Imaginary
H
0
0
0
0
P
0.7584
0.7416
0.7710
0.8479
KSSTAT
0.0299
0.0303
0.0295
0.0272Figure 3
The CDF of the data generated by CMM and DC terminal data. (a) Real part of current. (b) Imaginary part of current. (c) Real part of voltage. (d) Imaginary part of voltage.
(a)
(b)
(c)
(d)Table 2
KS test result of CDF of DC terminal.
Parameter
Current
Voltage
Real
Imaginary
Real
Imaginary
H
0
0
0
0
P
0.6559
0.8298
0.9853
0.7184
KSSTAT
0.0353
0.0300
0.0208
0.0317


## 2.1. Harmonic Equivalent Analysis Model

Norton equivalent circuit of the utility side and customer side is shown in Figure1.Figure 1
The Norton equivalent circuit.In Figure1, I˙s is the utility harmonic current, and Zs is the harmonic impedance on the utility side. U˙PCC and I˙PCC are the harmonic voltage and current measured at the PCC, respectively.According to the Norton equivalent circuit, the equation can be written as(1)ZsI˙s=U˙PCC+ZsI˙PCC.We have(2)ΔI˙sn=I˙sn+1−I˙sn,ΔI˙PCCn=I˙PCCn+1−I˙PCCn,ΔU˙PCCn=U˙PCCn+1−U˙PCCn,n=1,2,…,N−1,where, n is the location of the sampling point, and N is the number of sampling points. ΔI˙sn is the fluctuation of the utility harmonic current source, ΔI˙PCCn is the harmonic current fluctuation at the PCC, and ΔU˙PCCn is the harmonic voltage fluctuation at the PCC.Then, equation (1) can be expressed as(3)ZsΔI˙s=ΔU˙PCC+ZsΔI˙PCC,where(4)ΔI˙s=ΔI˙s1,ΔI˙s2,⋯,ΔI˙sN−1,ΔU˙PCC=ΔU˙PCC1,ΔU˙PCC2,⋯,ΔU˙PCCN−1,ΔI˙PCC=ΔI˙PCC1,ΔI˙PCC2,⋯,ΔI˙PCCN−1.Thus, equation (3) can be rewritten as(5)Zs−x−Zs−yZs−yZs−xΔIs−xΔIs−y=ΔUPCC−xΔUPCC−y+Zs−x−Zs−yZs−yZs−xΔIPCC−xΔIPCC−y,where ΔIs−x and ΔIs−y are the real and imaginary parts of the utility harmonic current, respectively. ΔUPCC−x and ΔUPCC−y are the real and imaginary parts of the harmonic voltage measured at the PCC, respectively. ΔIPCC−x and ΔIPCC−y are the real and imaginary parts of the harmonic current measured at the PCC, respectively. Zs−x and Zs−y are the real and imaginary of the harmonic impedance on the utility side, respectively.Therefore,(6)ΔUs=ΔUPCC+Z^sΔIPCC,where(7)ΔUs=ΔUsxΔUsyT=Zs−x−Zs−yZs−yZs−xΔIs−xΔIs−y,ΔUPCC=ΔUPCC−xΔUPCC−yT,Z^s=Zs−x−Zs−yZs−yZs−x,ΔIPCC=ΔIPCC−xΔIPCC−yT.


## 2.2. Cauchy Mixed Model

AD-dimensional Cauchy probability density function is(8)φxμ,∑=Γ1+D/2Γ1/2πD/2∑1/21+x−μT∑−1x−μ1+D/2,where μ is the data center point, ∑ is the symmetric semidefinite matrix, D is the data dimension, and Γ⋅ is the Gamma function.If given observation dataX=x1,x2,…,xN, the CMM is(9)PXα,μ,∑=∑k=1KαkφXμk,Σk,where K, α=α1,α2,⋯,αK, μ=μ1,μ2,…,μK, and ∑=∑1,∑2,…,∑K are the number of models, weighted values, data center points, and covariance matrices, respectively.


## 2.3. Feasibility Analysis of the Model

This section introduces the fitting methods of data distribution and the test of fitting results by CMM. First, we introduce the Copula function and Sklar’s theorem [23, 24]:Definition 1.
A copula is a functionC of D variables on the unit D-cube [0,1]D with the following properties:(1)
The range ofC is the unit interval [0,1](2)
C(x) is zero for all x = (x1, …, xD) ∈ [0,1]D for which at least one coordinate equals zero(3)
C(x) = xd if all coordinates of x are 1 except the dth one(4)
C is D-increasing in the sense that for every a ≤ b in [0,1]D, the measure ΔCba assigned by C to the D-box [a, b] = [a1, b1] × … × [aD, bD] is nonnegative, i.e.,(10)ΔCba=∑ε1,…,εD∈0,1−1ε1+…+εDCε1a1+1−ε1b1,…,εDaD+1−εDbD≥0.
The key theorem due to Sklar clarifies the relations of dependence and the copula of a distribution:Theorem 1.
LetF denote a D-dimensional distribution function with margins F1, …, FD. Then, there exists a D-copula C such that for all real x = (x1, …, xD),(11)Fx=CF1x1,…,FDxD.
Therefore, theC is unique, that is,(12)Cx=FF1−1x1,…,FD−1xD,where F1−1⋅,…,FD−1⋅ are quasi-inverses of the marginal distribution functions.
To verify whether the CMM can fit the probability distribution of the harmonic voltageU˙PCC and current I˙PCC data, we combine the two-dimensional t-copula function [25] with the Monte Carlo method to generate a set of random data conforming to CMM. Then, the Cumulative Distribution Functions (CDF) of field data and generated data are calculated and expressed as curves, respectively. Finally, we use the Kolmogorov–Smirnov (KS) test [26] to detect the maximum distance between two CDF curves.
The steps of generating the CMM data are as follows:(1)
K models are selected, and the CMM is adopted to PCC data. Then, the corresponding weightαk and covariance matrix ∑k are obtained.(2)
Convert the obtained covariance matrix∑k into a correlation matrix ρk.(3)
Thet-Copula function is used to generate random two-dimensional Cauchy number of K groups associated with ρk.(4)
Generate N random numbersni from 0 to maxαk. If ni<αk, select the random number set of the kth two-dimensional Cauchy model.
According to the abovementioned method, the CMM is used to fit the distribution of two groups’ field data. The first set of field data is conducted at the site where a 150 kV bus bar feeds a 100 MW Electric Arc Furnace (EAF). The total sample time is 10 hours, and the calculated 3rd harmonic sampling values are averaged per minute, then the number of the all samples is 600. CMM random numbers are generated corresponding to the EAF data. The CDF of the data generated by CMM and the EAF data are shown in Figure2.
The results of the KS test are shown in Table1.
In Table1, H = 0 indicates acceptance of the hypothesis at the significance level of 5%, P is the probability of acceptance, and KSSTAT is the value of the test statistic. It can be seen from Figure 2 and Table 1 that the CMM can fit EAF data.
The second set of field data is taken from a city power grid with multiple DC terminals, and the measured data are collected from 600 data points (20 sampled data per minute) at a 500 kV bus of the DC terminal for the 11th harmonic voltage and current. CMM random numbers are generated corresponding to the DC terminal data. The CDF of the data generated by CMM and DC terminal data are shown in Figure3.
The results of the KS test are shown in Table2.
The meanings of the parameters in Table2 are the same as those in Table 1. According to Figure 3 and Table 2, it is obvious that the CMM can fit DC terminal data even though there is a slight deviation between the CDF curve of the CMM and that of DC terminal.Figure 2
The CDF of the data generated by CMM and the EAF data. (a) Real part of current. (b) Imaginary part of current. (c) Real part of voltage. (d) Imaginary part of voltage.
(a)
(b)
(c)
(d)Table 1
KS test result of CDF of EAF.
Parameter
Current
Voltage
Real
Imaginary
Real
Imaginary
H
0
0
0
0
P
0.7584
0.7416
0.7710
0.8479
KSSTAT
0.0299
0.0303
0.0295
0.0272Figure 3
The CDF of the data generated by CMM and DC terminal data. (a) Real part of current. (b) Imaginary part of current. (c) Real part of voltage. (d) Imaginary part of voltage.
(a)
(b)
(c)
(d)Table 2
KS test result of CDF of DC terminal.
Parameter
Current
Voltage
Real
Imaginary
Real
Imaginary
H
0
0
0
0
P
0.6559
0.8298
0.9853
0.7184
KSSTAT
0.0353
0.0300
0.0208
0.0317


## 3. Utility Harmonic Impedance Estimation

In the previous section, the feasibility of the CMM method is verified. In this section, the impedance estimation method based on CMM is described.For the observation dataΔUPCC=ΔUPCC1,ΔUPCC2,…,ΔUPCCN−1ΔIPCC=ΔIPCC1,ΔIPCC2,…,ΔIPCCN−1, the joint probability function of the CMM under utility impedance Z^s=Z^s1,Z^s2,…,Z^sK is(13)Pθ=PΔUPCC,ΔIPCCθ=∏n=1N−1PΔUPCCn,ΔIPCCnθ,where θ=θ1,θ2,…,θK,θk=μk;∑k;Z^k. By taking the logarithm likelihood function for Pθ, we know(14)lθ=argmaxθlogPθ=argmaxθlog∏n=1N−1PΔUPCCn,ΔIPCCnθ,=argmaxθ∑n=1N−1logPΔUPCCn,ΔIPCCnθ,=argmaxθ∑n=1N−1log∑k=1KαkφΔUPCCn,ΔIPCCnθk.It is very difficult to directly find the value ofθ from equation (14). Since the EM algorithm has stable estimation ability for the mixture model, the EM algorithm is applied here. The EM algorithm consists of two steps, namely, E-step: calculating the expectation of the log likelihood function of the complete data and M-step: obtaining a new parameter value by maximizing the value of the intermediate amount.E-step: according to Bayesian formula, we can find the posterior probability of thenth group of samples ΔUPCCn,ΔIPCC,n generated by the kth Cauchy model.(15)γnk=αkφΔUPCC,ΔIPCCθk∑k=1KαkφΔUPCC,ΔIPCCθk.For equation (15), γnk is a known probability value and satisfies equation (16).(16)∑k=1Kγnk=1.Then, we obtain the following functionlθ based on the principle of Jensen’s inequality:(17)lθ=argmaxθ∑n=1N−1log∑k=1KγnkαkφΔUPCCn,ΔIPCCnθkγnk≥argmaxθ∑n=1N−1∑k=1KγnklogαkφΔUPCCn,ΔIPCCnθkγnk.By adjusting the lower bound function of equation (17), the maximum value of the lower bound function can be obtained, so that equation (17) can gradually approach the optimal solution or obtain the local optimal solution as follows:(18)lθ=argmaxθ∑n=1N−1∑k=1KγnklogαkφΔUPCCn,ΔIPCCnθγnk,=argmaxθ∑n=1N−1∑k=1KγnklogαkΓ1+D2−logΓ12πD/2−logγnk−logΣk1/2−1+D2log1+ΔUPCCn+Z^skΔIPCCn−μkT⋅Σk−1ΔUPCCn+Z^skΔIPCCn−μk.M-Step is to calculate a new round of the iterative parameter model. We assume that(19)Ψnk=1+ΔUPCCn+Z^skΔIPCCn−μkTΣk−1ΔUPCCn+Z^skΔIPCCn−μk.So, the partial derivative of equation (18) for μk is(20)∂lθ∂μk=−∑n=1N−1γnk⋅1+D2⋅1Ψnk⋅Σk−1+Σk−TΔUPCCn+Z^skΔIPCCn−μk.Let equation (20) be equal 0; we have(21)μk=∑n=1N−1γnk⋅ΔUPCCn+Z^skΔIPCCnΨnk∑n=1N−1γnkΨnk−1.Since equation (21) is nonlinear, it is impossible to separate μk. Therefore, in this step, equation (21) can be iterated until convergence.The partial derivative of equation (18) for ∑k is obtained as follows:(22)∂lθ∂∑k=−1+D2∑n=1N−1γnkΨnk⋅∑k−1ΔUPCCn+Z^skΔIPCCn−μk⋅ΔUPCCn+Z^skΔIPCCn−μkTΣk−1−Σk−T.Let (22) be equal to 0, and after simplification, the update equation of ∑k is(23)∑k=1+D∑n=1N−1γnk/Ψnk⋅ΔUPCCn+Z^skΔIPCCn−μkΔUPCCn+Z^skΔIPCCn−μkT∑n=1N−1γnk.Similarly, since (23) is also nonlinear, we cannot separate ∑k. Thus, in this step, (23) can be iterated until convergence.SinceZ^sk is a negative diagonal matrix, its elements cannot be directly updated as a whole, and each element needs to be updated separately. The following transformation can be performed:(24)Z^skΔIPCCn=Zs−xk−Zs−ykZs−ykZs−xkΔIPCC−xnΔIPCC−yn,=ΔIPCC−xn−ΔIPCC−xnΔIPCC−xnΔIPCC−xnZs−xkZs−yk,=WnRk,where Wn=ΔIPCC−xn−ΔIPCC−xnΔIPCC−xnΔIPCC−xn and Rk=Zs−xkZs−ykT. The cost function for Rk is(25)lRk=∑n=1N−1∑k=1Kγnklog1+ΔUPCCn+WnRk−μkT∑k−1ΔUPCCn+WnRk−μk.The partial derivative of equation (25) for Rk is(26)∂lRk∂Rk=∑n=1N−1γnkWnT∑k−1ΔUPCCn−μk+WnT∑k−1WnRkΨnk.Let equation (26) be equal to 0; then, the update equation of Rk is(27)Rk=−∑n=1N−1γnk⋅WnT∑k−1WnRkΨnk−1∑n=1N−1γnk⋅WnT∑k−1UPCCn−μkΨnk.Furthermore,(28)lα=lθ−η∑k=1Kαk−1,η=N−1.We can obtain(29)αk=1N−1∑n=1N−1γnk.For equations (16) and (29),(30)∑k=1Kαk=1.Theμk,∑k,Rk, and αk pairs of parameters are updated by iterating over equations (21), (23), (27), and (29) until convergence. Since the harmonic data corresponding to each model affects the utility harmonic impedance, we use a weighted value αk for Rk as follows:(31)Zs=1,j⋅∑k=1KαkRk.The steps of this method are in Algorithm1.Algorithm 1: CMM algorithm.
Input: Harmonic voltageU˙PCC and current I˙PCC at PCC.
Output: Utility side harmonic impedance estimateZs.
Step1: According to equations Equations (1)–(3), (5), convert U˙PCC and I˙PCC in vector form to ΔUPCC and ΔIPCC in matrix form,
Step2: Select the appropriate model parameterKK≥1 and updating E-step parameter γnk by iteration equation (15),
Step3: Update M-step parametersμk,∑k,αk by equations (21), (23), and (29),
Step4: Update parameterRk by equation (27),
Step5: EstimateZs from equation (31).


## 4. Simulation Analysis

Simulations are conducted to verify the effectiveness of the proposed method in accordance of Figure1. They are programmed by Matlab2017a. Both the customer side and the utility side adopt a parallel model of harmonic current source and harmonic impedance. The specific parameters of the simulation model are set as follows: the amplitude of the utility harmonic impedance Zs is 10 Ω, and the phase angle is 80° (radian value is 1.3963). The amplitude of the customer harmonic impedance Zc is 10mm=Zc/ZsΩ, ⋅ is the module, and the phase angle is 33° (radian value is 0.5760). I˙s and I˙c are(32)I˙sn=200q⋅1+0.5sinn+λ1rsnexp−30⋅1+0.1csn−0.05π180j,I˙cn=200⋅1+0.5sinn+λ2rcnexp30⋅1+0.1ccn−0.05π180j,λ1=0.5;λ2=0.5;n=1,2,…,N,where, λ1 and λ2 are the disturbance factors belonging to [0,1], within the value range of λ1 and λ2, we have q≈Is/Ic, sin⋅ is a sine function, rsn and rcn are random numbers that obey the standard normal distribution, and csn and ccn are random numbers that follow a uniform distribution. Also, the unit is Ampere.According to the abovementioned parameters, 300 harmonic samples are generated randomly, and the utility harmonic impedance is estimated by the CMM, Fluctuation Quantity (FQ), Binary Regression (BR), Independent Random Vector (IRV), and Independent Component Analysis (ICA) methods, respectively.


### 4.1. Influence of Background Harmonic on Utility Impedance Estimation

To analyze the influence of background harmonic on the utility impedance estimation, we setq from 0.5 to 1.2 in steps of 0.1, and λ1 = λ2 = 0.5. When m = 3, the impedance of both sides is approximately equal. Here, we define the correlation coefficient ρ of the current of both sides as follows:(33)ρ=∑h=1NI˙sh−EI˙sI˙ch−EI˙c∗∑h=1NI˙sh−EI˙sI˙sh−EI˙s∗∑h=1NI˙ch−EI˙cI˙ch−EI˙c∗,where E⋅ is the expected value of the data and ∗ is the conjugate plural of a complex data. In the simulation, the ρ is about 0.43, and the relative errors of the harmonic impedance amplitude and the phase relative errors are shown in Figure 4, respectively.Figure 4
Relative errors of amplitude and phase ofZs (m = 3). (a) Amplitude relative errors. (b) Phase relative errors.
(a)
(b)It can be seen from Figure4, when the impedance of both sides is close, the relative errors of amplitude and phase of the five methods gradually increase with the increase of background harmonics. However, compared with the other four methods, the CMM method is the most reliable estimate of the utility harmonic impedance. When q ≥ 0.7, the FQ, BR, IRV, and ICA methods fail to estimate the utility impedance.Whenm = 10, the consumer impedance is much larger than the utility impedance, the value of impedance on each side is approximately equal. The amplitude relative errors and phase relative errors of the five methods are shown in Figure 5.Figure 5
Relative errors of amplitude and phase ofZs (m = 10). (a) Amplitude relative errors. (b) Phase relative errors.
(a)
(b)According to Figure5 and related analysis in the previous section, the five methods' estimation errors increase with the increase of amplitude ratio q. In addition, with the increase of current ratio q, the estimation accuracy of the CMM method is higher than that of other four methods. Therefore, the CMM method can effectively estimate the utility harmonic impedance when the background harmonic increases.


### 4.2. Influence of the Continuous Change of Background Harmonic and Harmonic Impedance on Utility Impedance Estimation

For further analyzing the influence of the continuous change of background harmonic and harmonic impedance on utility impedance estimation, setq ranges from 0.5 to 1.2 and the impedance ratio of both sides ranges from 1.5 to 15.0, and for every q (its step size is 0.05.) and m (its step size is 0.25). And set the disturbance factors λ1 = λ2 = 0.4, the correlation coefficient ρ of current of both sides is about 0.46. The relative errors of amplitude and phase of utility impedance are shown in Figure 6.Figure 6
Relative errors of amplitude and phase ofZs under the continuous change of background harmonic and harmonic impedance. (a) Amplitude relative errors. (b) Phase relative errors.
(a)
(b)According to Figure6, when the correlation coefficient ρ of the current amplitudes on each side is about 0.63, the relative error of the amplitude and phase angle of the CMM method is less than 12% and 10%, respectively, even under the condition that the impedances of both sides are approximately equivalent and the background harmonic is large, which verifies the robustness and effectiveness of the CMM method.


### 4.3. Influence of Correlation of Harmonic Currents of Both Sides on Utility Harmonic Impedance Estimation

To further verify the effectiveness of the CMM method when the currents of both sides are not independent or are strongly correlated, we conducted another simulation. Since changing the interference of both sides of the current can change the independence of the current, we first set the disturbance factorsλ1 = λ2 = λ and then set λ from 0.1 to 1.0 with a step size of 0.1. With the change in value of λ, the correlation coefficient ρ of the current of both sides is as shown in Table 3 (the current ratio q is fixed at 0.5, and the impedance ratio m is fixed at 5).Table 3
The value of the correlation coefficient (m = 5, q = 0.5).
Parameter
λ
ρ
λ
ρ
Value
0.1
0.9146
0.6
0.3697
0.2
0.8078
0.7
0.3093
0.3
0.6551
0.8
0.2322
0.4
0.5345
0.9
0.1228
0.5
0.4400
1.0
0.0936As can be seen from Table3, the correlation coefficient ρ decreases with the increase of disturbance parameter λ. Also, the amplitude relative errors and phase relative errors of the five methods are shown in Figure 7.Figure 7
Amplitude relative errors and phase relative errors of five algorithms with different λ values. (a) Amplitude relative errors. (b) Phase relative errors.
(a)
(b)From Figure7, when the correlation is strong between the currents of the consumer side and the utility side, the estimation error ratio of the ICA method on the utility harmonic impedance is the highest, followed by the FQ, BR, and IRV methods. In addition, it can be seen that the performance of the CMM method is significantly more stable, which does not depend on the correlation of currents on each side. Therefore, it proves that the CMM method provides reliable utility harmonic impedance regardless of the strong or weak correlation between the consumer side current and the utility side current.


## 4.1. Influence of Background Harmonic on Utility Impedance Estimation

To analyze the influence of background harmonic on the utility impedance estimation, we setq from 0.5 to 1.2 in steps of 0.1, and λ1 = λ2 = 0.5. When m = 3, the impedance of both sides is approximately equal. Here, we define the correlation coefficient ρ of the current of both sides as follows:(33)ρ=∑h=1NI˙sh−EI˙sI˙ch−EI˙c∗∑h=1NI˙sh−EI˙sI˙sh−EI˙s∗∑h=1NI˙ch−EI˙cI˙ch−EI˙c∗,where E⋅ is the expected value of the data and ∗ is the conjugate plural of a complex data. In the simulation, the ρ is about 0.43, and the relative errors of the harmonic impedance amplitude and the phase relative errors are shown in Figure 4, respectively.Figure 4
Relative errors of amplitude and phase ofZs (m = 3). (a) Amplitude relative errors. (b) Phase relative errors.
(a)
(b)It can be seen from Figure4, when the impedance of both sides is close, the relative errors of amplitude and phase of the five methods gradually increase with the increase of background harmonics. However, compared with the other four methods, the CMM method is the most reliable estimate of the utility harmonic impedance. When q ≥ 0.7, the FQ, BR, IRV, and ICA methods fail to estimate the utility impedance.Whenm = 10, the consumer impedance is much larger than the utility impedance, the value of impedance on each side is approximately equal. The amplitude relative errors and phase relative errors of the five methods are shown in Figure 5.Figure 5
Relative errors of amplitude and phase ofZs (m = 10). (a) Amplitude relative errors. (b) Phase relative errors.
(a)
(b)According to Figure5 and related analysis in the previous section, the five methods' estimation errors increase with the increase of amplitude ratio q. In addition, with the increase of current ratio q, the estimation accuracy of the CMM method is higher than that of other four methods. Therefore, the CMM method can effectively estimate the utility harmonic impedance when the background harmonic increases.


## 4.2. Influence of the Continuous Change of Background Harmonic and Harmonic Impedance on Utility Impedance Estimation

For further analyzing the influence of the continuous change of background harmonic and harmonic impedance on utility impedance estimation, setq ranges from 0.5 to 1.2 and the impedance ratio of both sides ranges from 1.5 to 15.0, and for every q (its step size is 0.05.) and m (its step size is 0.25). And set the disturbance factors λ1 = λ2 = 0.4, the correlation coefficient ρ of current of both sides is about 0.46. The relative errors of amplitude and phase of utility impedance are shown in Figure 6.Figure 6
Relative errors of amplitude and phase ofZs under the continuous change of background harmonic and harmonic impedance. (a) Amplitude relative errors. (b) Phase relative errors.
(a)
(b)According to Figure6, when the correlation coefficient ρ of the current amplitudes on each side is about 0.63, the relative error of the amplitude and phase angle of the CMM method is less than 12% and 10%, respectively, even under the condition that the impedances of both sides are approximately equivalent and the background harmonic is large, which verifies the robustness and effectiveness of the CMM method.


## 4.3. Influence of Correlation of Harmonic Currents of Both Sides on Utility Harmonic Impedance Estimation

To further verify the effectiveness of the CMM method when the currents of both sides are not independent or are strongly correlated, we conducted another simulation. Since changing the interference of both sides of the current can change the independence of the current, we first set the disturbance factorsλ1 = λ2 = λ and then set λ from 0.1 to 1.0 with a step size of 0.1. With the change in value of λ, the correlation coefficient ρ of the current of both sides is as shown in Table 3 (the current ratio q is fixed at 0.5, and the impedance ratio m is fixed at 5).Table 3
The value of the correlation coefficient (m = 5, q = 0.5).
Parameter
λ
ρ
λ
ρ
Value
0.1
0.9146
0.6
0.3697
0.2
0.8078
0.7
0.3093
0.3
0.6551
0.8
0.2322
0.4
0.5345
0.9
0.1228
0.5
0.4400
1.0
0.0936As can be seen from Table3, the correlation coefficient ρ decreases with the increase of disturbance parameter λ. Also, the amplitude relative errors and phase relative errors of the five methods are shown in Figure 7.Figure 7
Amplitude relative errors and phase relative errors of five algorithms with different λ values. (a) Amplitude relative errors. (b) Phase relative errors.
(a)
(b)From Figure7, when the correlation is strong between the currents of the consumer side and the utility side, the estimation error ratio of the ICA method on the utility harmonic impedance is the highest, followed by the FQ, BR, and IRV methods. In addition, it can be seen that the performance of the CMM method is significantly more stable, which does not depend on the correlation of currents on each side. Therefore, it proves that the CMM method provides reliable utility harmonic impedance regardless of the strong or weak correlation between the consumer side current and the utility side current.


## 5. Practical Cases

The feasibility and effectiveness of the proposed method are verified by simulation in Section4. In this section, the practicability of the proposed method is verified by two practical cases under different backgrounds.


### 5.1. Case 1

To verify the practicability of the proposed method, the field data of the EAF are used for verification, and the harmonic voltage and current amplitudes of the EAF are shown in Figure8. The estimation values of the five methods with 600 sample points are shown in Table 4.Figure 8
Amplitudes of Voltage and Current at the PCC of EAF. (a) Voltage amplitude. (b) Current amplitude.
(a)
(b)Table 4
Impedance values estimated by five algorithms for the EAF data.
Impedance estimation
Algorithm
CMM
FQ
BR
IRV
ICA
Amplitude (Ω)
10.1131
11.0819
10.1796
10.1199
10.1351
Phase (rad)
1.2220
1.2537
1.1847
1.2121
1.1697In Table4, the estimation results of the CMM, BR, IRV, and ICA methods are equivalent, while the amplitude estimation result of the FQ method is quite different from those of other four methods. To further verify the reliability of the proposed method in estimating the utility harmonic impedance, the abovementioned 600 data are divided into 11 subintervals. The first 10 subintervals, each of which has 54 sample points, and the last subinterval have 60 sample points. The five methods are used to estimate the utility impedance in per hour. The effectiveness of impedance estimation of the five methods is shown in Figure 9.Figure 9
Utility impedance estimation effect of five methods in per hour. (a) Utility impedance amplitude. (b) Utility impedance phase.
(a)
(b)As can be seen from Figure9, compared with the other four methods, the CMM method produces the most stable estimation results of the amplitude and phase radian of the utility harmonic impedance in each time period. The FQ method is unstable. Compared with the CMM method, the BR method still shows great fluctuations on its performances. The performance of the IRV method is roughly equivalent to the one of the CMM method, but in the 9th time period, the fluctuation of the IRV method is larger than that of the CMM method. This result proves the accuracy and robustness of the CMM method.


### 5.2. Case 2

The amplitudes of harmonic voltage and current of the DC terminal data are shown in Figure10. In the absence of load, the utility side is directly connected into a closed loop with wires. Then, the measured utility short-circuit capacity is 34580 MVA.Figure 10
Amplitudes of voltage and current at the PCC of DC terminal. (a) Voltage amplitude. (b) Current amplitude.
(a)
(b)The estimation values of the five methods with 600 sample points are shown in Table5.Table 5
Impedance values estimated by five algorithms for the DC terminal data.
Impedance estimation
Algorithm
CMM
FQ
BR
IRV
ICA
Amplitude (Ω)
79.1246
77.6837
78.6427
78.6716
77.9453
Phase (rad)
1.0636
1.0423
1.0607
1.0626
1.0403And then, the 11th harmonic data are divided into 10 subintervals, and each subinterval has 60 sample points. Figure11 shows the estimation values with per subinterval.Figure 11
Impedance estimation effect of five methods on utility side in each data segment. (a) Amplitude. (b) Phase.
(a)
(b)According to the short-circuit capacity, the approximate value of the utility harmonic impedance is 79.5 Ω. It can be seen from Table5 that the impedance amplitude estimated by the CMM method is closest to the rough estimation of short-circuit capacity. As can be seen from Figure 11(a), the estimation results of the CMM method in each data segment are more stable than those of the other four methods. The BR, IRV, and FQ methods have the same amplitude estimation fluctuation, but the amplitude estimated by the FQ method is the smallest. The fluctuation of the ICA method is relatively large. Figure 11(b) shows that the phase estimation values of the CMM, BR, and IRV methods are fluctuating around 1.0630 rad, and the proposed method is the most robust. The phase estimation of the FQ method is the smallest. The ICA method is unstable for phase estimation.


## 5.1. Case 1

To verify the practicability of the proposed method, the field data of the EAF are used for verification, and the harmonic voltage and current amplitudes of the EAF are shown in Figure8. The estimation values of the five methods with 600 sample points are shown in Table 4.Figure 8
Amplitudes of Voltage and Current at the PCC of EAF. (a) Voltage amplitude. (b) Current amplitude.
(a)
(b)Table 4
Impedance values estimated by five algorithms for the EAF data.
Impedance estimation
Algorithm
CMM
FQ
BR
IRV
ICA
Amplitude (Ω)
10.1131
11.0819
10.1796
10.1199
10.1351
Phase (rad)
1.2220
1.2537
1.1847
1.2121
1.1697In Table4, the estimation results of the CMM, BR, IRV, and ICA methods are equivalent, while the amplitude estimation result of the FQ method is quite different from those of other four methods. To further verify the reliability of the proposed method in estimating the utility harmonic impedance, the abovementioned 600 data are divided into 11 subintervals. The first 10 subintervals, each of which has 54 sample points, and the last subinterval have 60 sample points. The five methods are used to estimate the utility impedance in per hour. The effectiveness of impedance estimation of the five methods is shown in Figure 9.Figure 9
Utility impedance estimation effect of five methods in per hour. (a) Utility impedance amplitude. (b) Utility impedance phase.
(a)
(b)As can be seen from Figure9, compared with the other four methods, the CMM method produces the most stable estimation results of the amplitude and phase radian of the utility harmonic impedance in each time period. The FQ method is unstable. Compared with the CMM method, the BR method still shows great fluctuations on its performances. The performance of the IRV method is roughly equivalent to the one of the CMM method, but in the 9th time period, the fluctuation of the IRV method is larger than that of the CMM method. This result proves the accuracy and robustness of the CMM method.


## 5.2. Case 2

The amplitudes of harmonic voltage and current of the DC terminal data are shown in Figure10. In the absence of load, the utility side is directly connected into a closed loop with wires. Then, the measured utility short-circuit capacity is 34580 MVA.Figure 10
Amplitudes of voltage and current at the PCC of DC terminal. (a) Voltage amplitude. (b) Current amplitude.
(a)
(b)The estimation values of the five methods with 600 sample points are shown in Table5.Table 5
Impedance values estimated by five algorithms for the DC terminal data.
Impedance estimation
Algorithm
CMM
FQ
BR
IRV
ICA
Amplitude (Ω)
79.1246
77.6837
78.6427
78.6716
77.9453
Phase (rad)
1.0636
1.0423
1.0607
1.0626
1.0403And then, the 11th harmonic data are divided into 10 subintervals, and each subinterval has 60 sample points. Figure11 shows the estimation values with per subinterval.Figure 11
Impedance estimation effect of five methods on utility side in each data segment. (a) Amplitude. (b) Phase.
(a)
(b)According to the short-circuit capacity, the approximate value of the utility harmonic impedance is 79.5 Ω. It can be seen from Table5 that the impedance amplitude estimated by the CMM method is closest to the rough estimation of short-circuit capacity. As can be seen from Figure 11(a), the estimation results of the CMM method in each data segment are more stable than those of the other four methods. The BR, IRV, and FQ methods have the same amplitude estimation fluctuation, but the amplitude estimated by the FQ method is the smallest. The fluctuation of the ICA method is relatively large. Figure 11(b) shows that the phase estimation values of the CMM, BR, and IRV methods are fluctuating around 1.0630 rad, and the proposed method is the most robust. The phase estimation of the FQ method is the smallest. The ICA method is unstable for phase estimation.


## 6. Conclusions

In this paper, a novel method based on CMM for estimating utility harmonic impedance is proposed. In theory, CMM can approach any probability model, and in practice, it is easy to verify that field harmonic data of PCC can effectively be matched with CMM. Therefore, we convert the Norton equivalent circuit model at the PCC into a CMM, and through adjusting the parameters of the CMM with the EM algorithm, the optimal estimation of utility harmonic impedance is calculated. The distinguishing feature of the method can obtain accurate and stable estimation results without any assumption conditions. Simulations and field cases verify the effectiveness of the proposed method.



---

*Source: 1580475-2020-03-25.xml*