Study on Harmonic Impedance Estimation Based on Gaussian Mixture Regression Using Railway Power Supply Loads

: There are a huge number of harmonics in the railway power supply system. Accurately estimating the harmonic impedance of the system is the key to evaluating the harmonic emission level of the power supply system. A harmonic impedance estimation method is proposed in this paper, which takes the Gaussian mixture regression (GMR) as the main idea, and is dedicated to calculating the harmonic impedance when the load changes or the background harmonic changes in the traction power supply system. First, the harmonic voltages and currents are measured at the point of common coupling (PCC); secondly, a Gaussian mixture model (GMM) is established and optimized parameters are obtained through the EM algorithm; ﬁnally, a Gaussian mixture regression is performed to obtain the utility side harmonic impedance. In the simulation study, different harmonic impedance estimation models with uniform distribution and Gaussian distribution are established, respectively, and the harmonic impedance changes caused by different system structures in the railway power supply system are simulated. At the same time, the error is compared with the existing method to judge the accuracy and robustness of this method. In the case analysis, the average value, average error, standard deviation and other indicators are used to evaluate this method. Among them, the average error and standard deviation of this method are about one-ﬁfth to one-third of those of the binary linear regression (BLR) method and the independent random vector (IRV) method. At the same time, its index is slightly better than that of the support vector machine (SVM) method.


Introduction
Electrified railways contain a large number of impactful, single-phase and nonlinear loads. Harmonics have become common power quality problems in the traction power supply system [1,2]. The problem of harmonic pollution has become an important problem in the research of the traction power supply system. Serious harmonic pollution will affect the stable operation of electrified railways, and even cause equipment damage [3,4].
Internationally, harmonic pollution is often reduced through "reward and punishment programs", that is, relevant policies are introduced to impose economic penalties on harmonic emission sources, and at the same time, harmonic victims are financially compensated [5,6]. In heavy load, this kind of scheme has some room for "negotiation". Therefore, the division of responsibilities for harmonics is very important, which affects the implementation of the "reward and punishment programs". At the same time, the correct acquisition of harmonic impedance is the main difficulty in the division of harmonic responsibility [7,8].
Harmonic impedance estimation is a power quality analysis problem, which has been widely concerned by experts and scholars. The "non-intervention methods" are often used in the calculation of harmonic impedance of traction power supply system. Harmonic voltages and currents at the PCC are used to estimate the harmonic impedances, in this (1) The fluctuation method directly uses the ratio of harmonic voltage change rate to harmonic current change rate to replace the harmonic impedance at a specific time [17]. (2) The linear regression method takes harmonic voltage and current as variables, and the impedance can be obtained by solving the regression parameters [18]. (3) Independent random vector method, which takes advantage of the property of the probability theory that the covariance of two independent random vectors is 0. The impedance value is solved by canceling the background harmonic variation term [19]. (4) The blind source separation method uses the statistical principle to separate the harmonic mixed signal, and solves the harmonic impedance by separating the signal [20]. (5) The machine learning method, the internal connection between the input value and the output value is established. Taking the harmonic voltage and harmonic current collected at the PCC as input samples and the reference harmonic impedance as the output value, the calculation model is obtained and the harmonic impedance is solved [21].
In addition to the above "non-intervention" method, the harmonic impedance of railway power supply system can also be obtained by the "intervention" method [22][23][24], that is, the harmonic impedance is estimated by controlling the switching of each branch or injecting harmonic current into the traction power supply system [25]. Such methods will affect the stability of the system operation, and the actual use will be limited by the on-site environment.
Obviously, the existing methods all have various deficiencies that have been shown by research. For example, the BLR method mentioned in the literature [26] has high requirements on the measurement accuracy of the harmonic voltage and current, and is greatly affected by the background harmonic fluctuation. Reference [27] proposed the IRV method, which needs to determine whether the customer's harmonic source is connected with the harmonic source on the utility side. If there is a strong connection, this method is not practical. Reference [28] has proposed the SVM method, which may get trapped in local optima during operation and will affect the stability of the results.
To sum up, the harmonic impedance calculation of the railway power supply system is still restricted by problems such as the large fluctuation of background harmonics and the uncertainty of harmonic source correlation. Therefore, to overcome these problems in harmonic impedance calculation, a method of measuring the harmonic impedance of the railway power supply system via GMR is proposed [29]. GMR has been widely used in time series, automatic control, image denoising and other fields [30]. The proposed method works well in the case of uncertain harmonic source correlation, and has certain progress in timeliness. The accuracy and stability of the proposed method are verified via simulation studies and case analysis.
The work arrangement of each part of this paper is as follows. In Section 2, the basic network structure of the railway power supply system is first reviewed, and the basic model of harmonic impedance estimation is introduced according to the Norton equivalent circuit. Then, GMM and GMR are introduced, and the algorithm description and calculation steps are given. In Section 3, the harmonic impedance estimation models of uniform distribution and Gaussian distribution are built, respectively, and the harmonic impedance of the two models is calculated by the method proposed in this paper and three other existing methods. The conclusions of the three methods were compared and the accuracy and stability of the proposed method were analyzed. In Section 4, the proposed method is validated with a case study of a railway power supply system. Finally, conclusions are drawn in Section 5.

Basic Structure of Electrified Railway Power Supply System
The power source in the electrified railway power supply system is the secondary transmission bus of the step-up transformer of the power plant. The whole system consists of dedicated high-voltage transmission lines, traction substations, feeders, contact system, return lines and electric locomotives. Its structure is shown in Figure 1. Among them, the traction substation is powered by a two-phase 110 kV bus, and a traction transformer is used to reduce the voltage to a voltage level suitable for the operation of electric locomotives. The circuit is then formed through the contact system, electric locomotives, tracks and return lines.
of uniform distribution and Gaussian distribution are built, respectively, and the harmonic impedance of the two models is calculated by the method proposed in this paper and three other existing methods. The conclusions of the three methods were compared and the accuracy and stability of the proposed method were analyzed. In Section 4, the proposed method is validated with a case study of a railway power supply system. Finally, conclusions are drawn in Section 5.

Basic Structure of Electrified Railway Power Supply System
The power source in the electrified railway power supply system is the secondary transmission bus of the step-up transformer of the power plant. The whole system consists of dedicated high-voltage transmission lines, traction substations, feeders, contact system, return lines and electric locomotives. Its structure is shown in Figure 1. Among them, the traction substation is powered by a two-phase 110 kV bus, and a traction transformer is used to reduce the voltage to a voltage level suitable for the operation of electric locomotives. The circuit is then formed through the contact system, electric locomotives, tracks and return lines.  Figure 1 shows the structure diagram of the electrified railway system. Among them, the traction substation and rolling stock on the secondary side can be equivalent to the customer side, and the 110 kV busbar and the upstream system can be equivalent to the utility side. In traction power systems, the Norton equivalent model is one of the most commonly used equivalents, as shown in Figure 2, where s Z and c Z represent the equivalent harmonic impedances on the utility side and customer side, respectively. s I represents the harmonic current value emitted by the utility side. c I represents the harmonic current value emitted by the customer side.  Figure 1 shows the structure diagram of the electrified railway system. Among them, the traction substation and rolling stock on the secondary side can be equivalent to the customer side, and the 110 kV busbar and the upstream system can be equivalent to the utility side. In traction power systems, the Norton equivalent model is one of the most commonly used equivalents, as shown in Figure 2, where · Z s and · Z c represent the equivalent harmonic impedances on the utility side and customer side, respectively. · I s represents the harmonic current value emitted by the utility side. · I c represents the harmonic current value emitted by the customer side. According to the basic circuit principle, the following equations can be obtained:

Harmonic Impedance Study
where pcc U is the harmonic voltage measured at the PCC; pcc I is the harmonic current measured at the PCC. Arranging Equations (1) and (2), the following equations can be According to the basic circuit principle, the following equations can be obtained: Energies 2022, 15,6952 The main idea of GMM is to classify the parameters of different probability density functions and accurately quantify things. Assuming that each sample of the GMM represents a Gaussian distribution, the GMM model can be regarded as a weighted average of multiple Gaussian distributions [31], that is, multiple single Gaussian models are weighted to approximate an arbitrary continuous distribution.
The probability density function of the GMM model of the random variable x is: where r represents the number of mixture Gaussian components, r = 1, 2, · · · , R; ϕ r represents the mixture weight of each Gaussian component, and 0 ≤ ϕ r ≤ 1 , ∑ R r=1 ϕ r = 1; N(x|µ r , ∑ r ) represents the distribution of the r-th Gaussian component; µ r represents the mean value of the r-th Gaussian component; the covariance of the r-th Gaussian component is denoted by ∑ r .The dimension of the random variable x determines the dimension of N(x), and the proportion of each Gaussian component in the mixture of Gaussians is determined by the mixture weight ϕ r .
The hyperparametric solution of GMM usually uses the expectation maximization (EM) algorithm. The core idea is to approach the objective function by continuously maximizing the lower bound function. First, the values of the model parameters are estimated from the known data samples. Then, estimate the value of missing data according to the estimated parameter value. Finally, the parameter values are re-estimated based on the estimated missing data values and the known sample data values. The above steps are repeated continuously, and the iterations are repeated until the final convergence ends.
The core of the EM algorithm is mainly divided into two steps: E − step: Guess the value of estimated z (i) , and predict that the sample parameters may come from a Gaussian component.
where ϕ j is the weight parameter; µ j is the mathematical expectation; ∑ j is the covariance matrix. The hyperparameters of the GMM model can be obtained by iteratively calculating E − step and M − step, which is expressed as:

Gaussian Mixture Regression
The EM algorithm is used to calculate the established GMM model, and the hyperparameter Q of each group of Gaussian distribution is obtained. At the same time, GMM is performed according to the linear combination of Gaussian distribution and Gaussian condition. The data vector consists of input data and output data. The input data consists of a two-dimensional vector, and the output data consists of a one-dimensional vector. Then, the Gaussian mixture joint probability density can be expressed by Equation (10): where µ n represents the two-dimensional mean vector; ∑ n is the covariance vector group. The conditional mean and variance are: The regression function expression is: The conditional variance function is: According to Formula (13), a regression model can be obtained, and a new output sample can be obtained by substituting the input sample.

Calculation of Harmonic Impedance by GMR
In this paper, GMR is used for harmonic impedance calculation. The specific process is as follows: (1) Collect harmonic samples at PCC and take the vector group consisting of · U pcc and · I pcc as the input sample. Meanwhile, the output sample consists of a one-dimensional harmonic vector · Z s . To sum up, the expressions of input variables and output variables are as follows: (2) Establish the combined density GMM model of · U pcc , · I pcc and · Z s . (3) Use the EM algorithm to obtain the optimal parameters of the GMM model. (4) Input the test data to calculate the conditional mean and variance of the test set shown in Equations (11) and (12).
(5) According to the calculation results of GMR output, evaluate the indicators of this method, and compare the calculation results of other methods.
In summary, the specific process of the proposed method is shown in Figure 3.

Simulation Study
In the simulation study, four methods are applied to estimate harmonic impedance and Matlab software is used to establish a simulation model. Method 1: GMR method which is the method in the text; Method 2: BLR method; Method 3: IRV method; Method 4: SVM method.

Simulation I: Harmonic Impedance Estimation with Gaussian Distribution
In the light of the equivalent model in Figure 1, a Gaussian distributed harmoni impedance estimation model is established; meanwhile, the following parameters are se in Simulation I.
(1) Harmonic impedance on the utility side is s Z , which obeys a Gaussian distribution with mean   (4) Harmonic current emitted by the utility side is I , which is k times larger than I

Simulation Study
In the simulation study, four methods are applied to estimate harmonic impedance, and Matlab software is used to establish a simulation model. Method 1: GMR method, which is the method in the text; Method 2: BLR method; Method 3: IRV method; Method 4: SVM method.

Simulation I: Harmonic Impedance Estimation with Gaussian Distribution
In the light of the equivalent model in Figure 1, a Gaussian distributed harmonic impedance estimation model is established; meanwhile, the following parameters are set in Simulation I.
(1) Harmonic impedance on the utility side is · Z s , which obeys a Gaussian distribution with mean 3 + j8 and variance 0.3 + j0.7, and the unit is Ω.
(2) Harmonic impedance on the customer side is · Z c , which obeys a Gaussian distribution with mean 25 + j168 and variance 0.8 + j11, and the unit is Ω.
(3) Harmonic current emitted by the customer side is · I c , which obeys a Gaussian distribution with mean 7 + j12 and variance 0.4 + j0.7, and the unit is A. (4) Harmonic current emitted by the utility side is · I s , which is k times larger than · I c , and the value interval of k is 0.2~1.0, and the value interval is 0.2.
The number of samples collected in simulation I is 1440. Based on the sample set {x n , y n } 1440 n=1 , four methods are used in Simulation I to calculate harmonic impedance. It is worth noting that methods 2 and 3 require the use of continuous calculation, i.e., continuous calculation of harmonic impedance every 200 sampling points.
According to the calculation results of each method, 1000 calculation results are taken and each 100 points are taken as a data group, and a box plot is drawn, where "n" represents the number of data sets. The extreme values, medians, quartiles, and outliers for each group of results are visualized in the image. Since k has different values, this paper only shows the calculation result images of each method when k = 0.2, as shown in        Comparing Figures 3-6, it is not difficult to see that the results of method 1 and method 4 are less volatile. After careful comparison, it is found that the fluctuation range of the mean value of each segment in method 1 is smaller than that in method 4. Meanwhile, the fluctuation range of extreme values in method 1 is obviously smaller than that in method 4. Therefore, in this simulation, method 1 has the best stability among the four methods. To further emphasize the superiority of the proposed method, the error histogram of several methods when the value of k is different is made, in which the error size takes the absolute value, as shown in Figure 8.
Energies 2022, 15,6952 Comparing Figures 3-6, it is not difficult to see that the results of method method 4 are less volatile. After careful comparison, it is found that the fluctuation of the mean value of each segment in method 1 is smaller than that in met Meanwhile, the fluctuation range of extreme values in method 1 is obviously s than that in method 4. Therefore, in this simulation, method 1 has the best st among the four methods. To further emphasize the superiority of the proposed m the error histogram of several methods when the value of k is different is made, in the error size takes the absolute value, as shown in Figure 8. As shown in Figure 8, if the value of k increases, the errors of method 2 and m 3 also increase. The error of the real impedance part of method 2 is smaller than method 3, and the error of the imaginary impedance part of method 3 is smaller tha of method 2. The error of method 1 and method 4 fluctuates very little when k ch When k = 0.2, the error of method 2 is close to that of method 4. Since method 1 h error in the model, a small box is used in the figure to show its error. When k different values, the mean value of the real part harmonic impedance obtained th various methods is recorded in Table 1, and the mean value of the imaginar harmonic impedance is recorded in Table 2, in which three decimal places are rese  As shown in Figure 8, if the value of k increases, the errors of method 2 and method 3 also increase. The error of the real impedance part of method 2 is smaller than that of method 3, and the error of the imaginary impedance part of method 3 is smaller than that of method 2. The error of method 1 and method 4 fluctuates very little when k changes. When k = 0.2, the error of method 2 is close to that of method 4. Since method 1 has less error in the model, a small box is used in the figure to show its error. When k takes different values, the mean value of the real part harmonic impedance obtained through various methods is recorded in Table 1, and the mean value of the imaginary part harmonic impedance is recorded in Table 2, in which three decimal places are reserved.  Based on the results in Tables 1 and 2, it is clear that, method 1 is the most stable and most accurate of the several methods. The real part and imaginary part of the calculation result of method 2 both increase with the increase in k value. In the error results of method 2, when k = 1, the maximum error of the real part is 64%, and the maximum error of the imaginary part is 65%. The trend of the calculation results of method 3 is consistent with that of method 2. When k = 1, the real part error of method 3 reaches the maximum value of 77%, and the imaginary part error reaches the maximum value of 43%. Method 4 is also very stable in this model, but slightly less accurate than method 1.

Simulation II: Uniformly Distributed Harmonic Impedance Estimation
In the light of the equivalent model in Figure 2, a uniformly distributed harmonic impedance estimation model is established. At this time, there is no correlation between the harmonic sources at both ends of the system. The following parameters are set in Simulation II.
(1) Harmonic impedance on the utility side is · Z s , which is uniformly distributed with mean 5 + j20. Add a standard deviation of 0.01 to the real part and a standard deviation of 0.034 to the imaginary part, and the unit is Ω.
(2) Harmonic impedance on the customer side is · Z c , which is uniformly distributed with mean 40 + j296. Add a standard deviation of 0.662 to the real part and a standard deviation of 1.664 to the imaginary part, and the unit is Ω.
(3) Harmonic current emitted by the utility side is · I s , which is set to a fixed value 0.75 + j0.75, and the unit is A.
(4) Harmonic current emitted by the customer side is · I c , which is uniformly distributed with mean 3.73 + j3.74. Add a standard deviation of 0.242 to the real part and a standard deviation of 0.238 to the imaginary part and the unit is A.
Consistent with Simulation I, the harmonic impedance in Simulation II can be calculated in four ways. Similarly, the calculation uses continuous calculation every 200 sampling points. In the calculation results of each method, take 1000 calculation results and use every 100 points as a data group to draw a box plot as shown in Figures 9-12. Energies 2022, 15,6952 of the imaginary part is 65%. The trend of the calculation results of method 3 is co with that of method 2. When k = 1, the real part error of method 3 reaches the ma value of 77%, and the imaginary part error reaches the maximum value of 43%. M is also very stable in this model, but slightly less accurate than method 1.

Simulation II: Uniformly Distributed Harmonic Impedance Estimation
In the light of the equivalent model in Figure 2, a uniformly distributed ha impedance estimation model is established. At this time, there is no correlation b the harmonic sources at both ends of the system. The following parameters ar Simulation II.
(1) Harmonic impedance on the utility side is s Z , which is uniformly distribut mean 5+ 20 j . Add a standard deviation of 0.01 to the real part and a s deviation of 0.034 to the imaginary part, and the unit is Ω.
(2) Harmonic impedance on the customer side is c Z , which is uniformly dis with mean   Comparing Figures 8-11, method 1 and method 4 are relatively stable. A there are some outliers in the imaginary part of method 1, the fluctuation ampl still smaller than method 4, and the volatility of method 2 and method 3 is lar further highlight the superiority of the proposed method in this paper, th histogram of several methods is drawn, in which the absolute value of the error i as shown in Figure 12. Figure 13 shows that methods 1 and 4 are more accurate than methods 2 an worth noting that the errors of the four methods in the second simulation are ver all less than 1%. The average impedance values calculated by several meth Simulation II are shown in Table 3, and the data in the table are kept to three places. Comparing Figures 8-11, method 1 and method 4 are relatively stable. A there are some outliers in the imaginary part of method 1, the fluctuation ampl still smaller than method 4, and the volatility of method 2 and method 3 is lar further highlight the superiority of the proposed method in this paper, th histogram of several methods is drawn, in which the absolute value of the error i as shown in Figure 12. Figure 13 shows that methods 1 and 4 are more accurate than methods 2 and worth noting that the errors of the four methods in the second simulation are ver all less than 1%. The average impedance values calculated by several meth Simulation II are shown in Table 3, and the data in the table are kept to three places. Comparing Figures 8-11, method 1 and method 4 are relatively stable. Although there are some outliers in the imaginary part of method 1, the fluctuation amplitude is still smaller than method 4, and the volatility of method 2 and method 3 is larger. To further highlight the superiority of the proposed method in this paper, the error histogram of several methods is drawn, in which the absolute value of the error is taken, as shown in Figure 12. Figure 13 shows that methods 1 and 4 are more accurate than methods 2 and 3. It is worth noting that the errors of the four methods in the second simulation are very small, all less than 1%. The average impedance values calculated by several methods in Simulation II are shown in Table 3, and the data in the table are kept to three decimal places.
histogram of several methods is drawn, in which the absolute value of the error i as shown in Figure 12. Figure 13 shows that methods 1 and 4 are more accurate than methods 2 and worth noting that the errors of the four methods in the second simulation are very all less than 1%. The average impedance values calculated by several meth Simulation II are shown in Table 3, and the data in the table are kept to three d places.   Longitudinal comparison of data in Tables 1-3 in Simulation I and Simulation II shows that method 1 (the method in this paper) has significant superiority in estimating the harmonic impedance of Gaussian distribution. In simulation II, a good result is obtained when the sample fitting degree is not as good as Gaussian distribution, which demonstrates the generalization performance of this method. The results of method 2 and method 3 are closer in simulation II. The change trend of the two methods is the same as that of Simulation II, and with the increase in k, the error also increases. The reason for this phenomenon is that there is a correlation between the harmonic sources at both ends of the system. The results of method 4 demonstrate the advantages of the SVM method in simulations I and II. When there is a correlation between the harmonic sources on the utility side and the customer side, it will not be significantly affected, and is only slightly lower than method 1 in terms of accuracy and stability.

Case Analysis
Based on the 110kV traction power supply substation, the power quality analyzer is capable of accurately collecting harmonic voltages and currents at the PCC. A schematic diagram of harmonic measurement is shown in Figure 14. In the figure, the fundamental wave voltage amplitude test error of the PQ monitor is: <±0.5% of Full Scale (F.S.), and the fundamental wave current amplitude test error is: <±1.0% F.S. The unit measurement time interval is 3 s, and the total measurement time is 10 h. fast Fourier transform is used to obtain each harmonic, and this paper takes the 3rd harmonic as an example. The traction power supply load will inevitably have low load points and no-load points. The fundamental current is analyzed to remove such points that affect the calculation. Figure 15 shows 600 load samples, which were used as experimental targets, and the third harmonic voltage and current with load are represented. Figure 16 shows the phase angles of the third harmonic voltage and current with load, expressed in the radian system.

Case Analysis
Based on the 110kV traction power supply substation, the power quality analyzer is capable of accurately collecting harmonic voltages and currents at the PCC. A schematic diagram of harmonic measurement is shown in Figure 14. In the figure, the fundamental wave voltage amplitude test error of the PQ monitor is: 0.5%  of Full Scale (F.S.), and the fundamental wave current amplitude test error is: 1.0%  F.S. The unit measurement time interval is 3 s, and the total measurement time is 10 h. fast Fourier transform is used to obtain each harmonic, and this paper takes the 3rd harmonic as an example. The traction power supply load will inevitably have low load points and no-load points. The fundamental current is analyzed to remove such points that affect the calculation. Figure 15 shows 600 load samples, which were used as experimental targets, and the third harmonic voltage and current with load are Energies 2022, 15,6952 represented. Figure 16 shows the phase angles of the third harmonic voltage and with load, expressed in the radian system. In the case analysis, four methods consistent with the simulation study are calculate the 3rd harmonic impedance. The differences of each method in case should be more directly reflected, and the results should be display three-dimensional images. The differences of each method in the case study sh presented more directly, and the results should be presented in the f three-dimensional images. Mean value, mean error, standard deviation, ma deviation of mean and calculation time are used to evaluate the results. The cal formula of the evaluation index is as follows: Energies 2022, 15,6952 represented. Figure 16 shows the phase angles of the third harmonic voltage and c with load, expressed in the radian system. In the case analysis, four methods consistent with the simulation study are u calculate the 3rd harmonic impedance. The differences of each method in case a should be more directly reflected, and the results should be display three-dimensional images. The differences of each method in the case study sho presented more directly, and the results should be presented in the fo three-dimensional images. Mean value, mean error, standard deviation, max deviation of mean and calculation time are used to evaluate the results. The calcu formula of the evaluation index is as follows: In the case analysis, four methods consistent with the simulation study are used to calculate the 3rd harmonic impedance. The differences of each method in case analysis should be more directly reflected, and the results should be displayed by three-dimensional images. The differences of each method in the case study should be presented more directly, and the results should be presented in the form of three-dimensional images. Mean value, mean error, standard deviation, maximum deviation of mean and calculation time are used to evaluate the results. The calculation formula of the evaluation index is as follows: Among them, Error ave represents the average error; Error std represents the standard deviation; Error max represents the maximum deviation of the mean; Z s,i represents the calculated harmonic impedance; Z s represents the average harmonic impedance of this method.
In method 1, the harmonic impedance of the utility side is estimated by Gaussian mixture regression, and the result is obtained by interpolation. Figure 17 shows the calculation results of method 1, while Table 4 shows the index parameters of method 1.
Energies 2022, 15, 6952 14 In method 1, the harmonic impedance of the utility side is estimated by Gau mixture regression, and the result is obtained by interpolation. Figure 17 shows calculation results of method 1, while Table 4 shows the index parameters of method  Through methods 2, 3 and 4, the real and imaginary parts of harmonic imped are estimated, as shown in Figures 18-20, respectively. The parameters of var indicators are shown in Tables 5-7.    Through methods 2, 3 and 4, the real and imaginary parts of harmonic impedance are estimated, as shown in Figures 18-20, respectively. The parameters of various indicators are shown in Tables 5-7. Energies 2022, 15, 6952 14 In method 1, the harmonic impedance of the utility side is estimated by Gau mixture regression, and the result is obtained by interpolation. Figure 17 show calculation results of method 1, while Table 4 shows the index parameters of method          By comparing Figures 17-20 and Tables 4-7, it can be seen that each evalu index of method 1 is significantly better than those of the other three methods, excep calculation time. Due to the measured data from the actual traction power supply sy the linearity of the data is low, and the result of method 2 has a very large fluctua which is greatly affected by the background harmonics. Method 3 is the worst amon methods used in this paper. The maximum deviation of the mean of the real part is n 4 times higher than its mean, and the maximum deviation of the mean of the imag part is nearly 1 times higher than its mean. Because the harmonic source on the u side has a certain relationship with the harmonic source on the customer side method is abnormal. In method 4, some peak values appear in the operation, falling the local optimum, and the results show great instability.
In order to further reflect the universality of the method in this paper, th harmonic of the measured data is analyzed again. The 7th harmonic impedance obt by the four methods is shown in Figures 21-24. At the same time, the parameters o 7th harmonic impedance of the four methods are shown in Tables 8-11.      Tables 4-7, it can be seen that each evalu index of method 1 is significantly better than those of the other three methods, excep calculation time. Due to the measured data from the actual traction power supply sys the linearity of the data is low, and the result of method 2 has a very large fluctua which is greatly affected by the background harmonics. Method 3 is the worst amon methods used in this paper. The maximum deviation of the mean of the real part is n 4 times higher than its mean, and the maximum deviation of the mean of the imag part is nearly 1 times higher than its mean. Because the harmonic source on the u side has a certain relationship with the harmonic source on the customer side, method is abnormal. In method 4, some peak values appear in the operation, falling the local optimum, and the results show great instability.
In order to further reflect the universality of the method in this paper, the harmonic of the measured data is analyzed again. The 7th harmonic impedance obta by the four methods is shown in Figures 21-24. At the same time, the parameters o 7th harmonic impedance of the four methods are shown in Tables 8-11.  By comparing Figures 17-20 and Tables 4-7, it can be seen that each evaluation index of method 1 is significantly better than those of the other three methods, except the calculation time. Due to the measured data from the actual traction power supply system, the linearity of the data is low, and the result of method 2 has a very large fluctuation, which is greatly affected by the background harmonics. Method 3 is the worst among the methods used in this paper. The maximum deviation of the mean of the real part is nearly 4 times higher than its mean, and the maximum deviation of the mean of the imaginary part is nearly 1 times higher than its mean. Because the harmonic source on the utility side has a certain relationship with the harmonic source on the customer side, this method is abnormal. In method 4, some peak values appear in the operation, falling into the local optimum, and the results show great instability.
In order to further reflect the universality of the method in this paper, the 7th harmonic of the measured data is analyzed again. The 7th harmonic impedance obtained by the four methods is shown in Figures 21-24. At the same time, the parameters of the 7th harmonic impedance of the four methods are shown in Tables 8-11.                     Previously, the literature [20] proposed a method for harmonic impedance estimation via O-GPR. This method uses the Bayesian framework to optimize Gaussian  Previously, the literature [20] proposed a method for harmonic impedance estimation via O-GPR. This method uses the Bayesian framework to optimize Gaussian process regression, and the effect is very significant. It is not difficult to find that GMR results are similar to O-GPR results after visualizing the results. Likewise, the results of the two methods in the calculation of harmonic impedance are symmetrical on the positive and negative half axis. In terms of mean error, standard deviation and mean maximum deviation, GMR is close to O-GPR. However, in terms of running time, GMR is significantly better than O-GPR.

Conclusions
The background harmonics of the railway power supply system are characterized by strong fluctuation and correlation between harmonics sources. The existing methods of harmonic impedance estimation may be biased. A GMR method for estimating the harmonic impedance of the utility side of the railway power supply system is presented in this paper. First, the GMM model is constructed and the hyperparameters are solved by the EM algorithm to obtain the GMR model and the weights of each Gaussian component. Then, substitute the harmonic voltage and harmonic current data to solve the harmonic impedance. Finally, the validity is verified by two different harmonic impedance estimation models and an actual case of a railway power supply system. The following conclusions are drawn: (1) This method is not sensitive to the correlation of harmonic sources and is suitable for occasions where the correlation of harmonic sources is strong. When the correlation between the two harmonic sources changes, the method can obtain better evaluation results and has good generalization performance. (2) Compared with previous similar methods, such as the GPR method, its timeliness is also better optimized. (3) In the Gaussian distributed harmonic and uniformly distributed harmonic impedance models, the stability performance of the GMR method is better than that of the BLR, IRV, and SVM methods compared in this paper. (4) In the case analysis, the average error and standard deviation of this method are about one-fifth to one-third of those of the binary linear regression (BLR) method and the independent random vector (IRV) method. At the same time, its index is slightly better than that of the support vector machine (SVM) method.
Author Contributions: W.T. and Y.X. mainly carried out, completed this work, simulated and analyzed; Y.X. provided and analyzed the field data; Y.X. supplemented and perfected this work in the later stage; W.T. wrote the text; Y.X. and W.T. reviewed the text; all authors revised and approved the final manuscript. All authors have read and agreed to the published version of the manuscript.
Funding: This paper was supported in part by the Key Scientific Research fund of Xihua University (Z1520909), the Southwest jiaotong University Traction Power State Key Laboratory Open Issue (TPL1907).

Conflicts of Interest:
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.