A Method for Responsibility Division of Multi-Harmonic Sources Based on Canonical Correlation Analysis

: The variable background harmonic data and incomplete phasor information make multi-harmonic source responsibility division in three-phase symmetrical power system a signiﬁcant challenge. In this paper, a background harmonic data selection method based on canonical correlation analysis is proposed to deal with multi-harmonic source responsibility division without phasor information. Firstly, the canonical correlation coefﬁcient between harmonic voltage and harmonic current is used to characterize the ﬂuctuations of background harmonic voltage. Then, the sliding window method is adopted to select the harmonic voltage and harmonic current with small ﬂuctuations. Next, the canonical correlation results for selected data are used to calculate the harmonic responsibility index via the linear regression method. The harmonic responsibility index in the form of percentage represents the harmonic responsibility division. Finally, several experimental results demonstrate that the proposed method has a high accuracy in calculating the harmonic responsibility division, particularly when the user side contains ﬂuctuations of unknown harmonic sources.


Introduction
In the field of electronic and electrical engineering, nonlinear load is widely used, which makes the harmonic pollution of three-phase symmetrical power system increasingly serious. Nonlinear devices, such as generators, inverters, and rectifiers, inject harmonic currents into their connected power grid, resulting in voltage distortion at the point of common coupling (PCC), which in turn affects the current at this point. Moreover, the interactions among multiple harmonic loads are complicated, making it difficult to locate the sources of distortion. Therefore, to address the issue of harmonic pollution [1], accurately identifying the "problematic" harmonic sources and determining the responsibility of each harmonic source are two extremely important aspects, the latter of which being the focus of this paper.
Existing methods of harmonic responsibility division mainly include: the volatility method [2][3][4], the covariance method [5][6][7], the blind source separation method [8][9][10][11], and the linear regression method [12][13][14][15][16][17][18][19][20]. The first three methods have shown a higher accuracy in estimating the harmonic impedance, but they all require accurate harmonic waveform sampling data to obtain the harmonic amplitude and phase angle to realize the separation and calculation of the real and imaginary parts. The linear regression method is vulnerable to the fluctuation of the background harmonic. In addition, the regression formulas derived from following the circuit principles in different works may require different data, and some of them do not require phase angle data [19,20]. The volatility method and the covariance method are both proposed based on the responsibility assessment of a single harmonic source, while the blind source separation method and linear regression method can be used for a single harmonic source and multi-harmonic sources.
At present, responsibility division for multi-harmonic sources has been widely studied, and the proposed methods are generally based on accurately estimating the harmonic impedance [21]. The blind source separation method has been widely studied in recent years. The independent component method (ICA) as well as its improved variants [8][9][10][11] are the most representative. When the harmonic source meets the precondition of being statistical independence and non-Gaussianness, individual harmonic current signals can be separated one by one from their mixed signal with merely measuring the harmonic voltage and without knowing the network harmonic impedance. The least squares method and its improved variants [12][13][14][15] are a common method to obtain the harmonic impedance through calculating the coefficient of Norton equivalent circuit of each harmonic source [16].
However, in real situations where the assessment of harmonic impedance is affected by the highly autocorrelated harmonic monitoring data, the problem may become illconditioned, and the abovementioned methods could fail. The problem is highly possible to be ill-conditioned when multiple harmonic sources exist, because a large amount of data are required for such systems. Partial least-squares method is an improvement of the least squares method, which overcomes the disadvantage of variable dependence for systems model with multiple harmonic sources [17,19]. M-estimation robust regression method [20] can reduce the influence of singular values to a certain extent and realize the division without phase angel information, but this method is mainly used for responsibility division of a single harmonic source. When used for dividing responsibility of multi-harmonic sources, data segments in which the harmonic resource to be analyzed fluctuates greatly while the fluctuations of other harmonic sources being less than 5~10% have to be selected out. In practical cases, the limited available data segments would make the corresponding result difficult to reflect the overall situation of harmonic responsibility in a long time scale.
The relevant works of harmonic responsibility division in the past five years are summarized in Table 1, from which it can be seen that most methods of responsibility division of multi-harmonic sources require data of harmonic phase angle. However, due to the limitation on communication channels and storage capacity, the harmonic phase angle information is usually not stored in existing regional power quality online monitoring systems. To this end, to apply those methods, harmonic measuring instruments have to be used to conduct special tests on site to obtain short-term measurement data. However, the obtained data is not compatible to the monitoring system, and the corresponding results are hard to reflect the overall situation of long-term harmonic responsibility.
In general, the linear regression method can find highly accurate solutions when the utility harmonic voltage is constant; however, utility harmonic voltage usually varies in real systems [22]. When the utility harmonic voltage is under considerable fluctuations, the use of linear regression methods for harmonic impedance assessment becomes less accurate [23,24]. Data analysis methods can weaken the impacts of fluctuations of the background harmonic on the accuracy but may result in a reduction in the amount of available data. Taking the linear regression method as an example, we discuss three typical data analysis methods.

Method Data Requirement Data Selection Method Feature
Covariance method [5][6][7] Harmonic voltage amplitude; branch harmonic current amplitude; phase angle difference between harmonic voltage and harmonic current (without harmonic initial phase angle).
Not needed.
Applicable for harmonic responsibility division of a single harmonic source.
Data selection required; limited available data after selection.
Applicable for harmonic responsibility division of multi-harmonic sources.
Data selection required; limited available data after selection.
Applicable for harmonic responsibility division of one or multi-harmonic sources.

The proposed method
Harmonic voltage amplitude at PCC; branch harmonic current amplitude; (all available in grid online monitoring systems).
Data selection required; sufficient available data after selection.
Accurate even when the background harmonic fluctuates greatly; applicable for harmonic responsibility division of multi-harmonic sources.
The first solution is to classify the background harmonic voltage into different categories, and then perform linear regression for each category. For classification of the background harmonic voltage, Zang et al. [25] divided the harmonic data points into segments that take constant utility harmonic voltage, which can be determined from hierarchical K-means clustering. Zhang and coworkers [26] considered the effects of background harmonic fluctuations, and proposed a data selection method based on the improved Hampel weight function. The basic idea of this method is to estimate the utility background harmonic voltage by solving the phasor equation, which avoids the use of iterative methods. However, in real systems where the measurement of phase data is not mandatory, the background harmonic voltage cannot be estimated if the phase information is unknown. The accuracy of the utility background harmonic voltage relies on the accuracy of the utility harmonic impedance; its accuracy is usually low and difficult to calculate when the system harmonic impedance changes [27].
The second solution is to perform linear regressions first and then classifies the background harmonic voltage using iterative methods. As iterative methods are applied in this solution, the computational (both space and time) complexity is significantly increased [28]. The drawback of these two solutions is that both require the harmonic phase data, which are usually not available directly in power quality monitoring systems [29].
To avoid this limitation, the third solution based on data correlation is proposed to assess the harmonic responsibilities by comparing the similarity between the harmonic voltage and the harmonic current and then selecting the strongly relevant data for linear regression [30,31]. However, these methods cannot reflect the degree of correlation for multiple sets of data combinations and are used only for estimating the harmonic emission capability of harmonic sources on both sides of the PCC. For cases where multiple harmonic sources are connected to the PCC, it remains unsolved as to how to analyze the correlations between the harmonic currents and the harmonic voltage.
However, the above-mentioned multi-harmonic source division method requires a large amount of calculations and cumbersome steps, the responsibility division for multiharmonic sources needs further improvement.
The main contributions of this paper can be highlighted as follows: (1) A harmonic responsibility division method based on canonical correlation analysis is proposed. Compared with the existing methods, harmonic phase angle data is not required while still achieving a higher division accuracy even when the background harmonic fluctuates relatively greatly. (2) One-to-many correlation analysis between the harmonic voltage at PCC and the data of multi-harmonic sources is realized with canonical correlation analysis method, with which even under great fluctuations of the background harmonic, appropriate data can be selected out. (3) Merely using available data from the online power quality monitoring system and without involving additional harmonic measuring instruments, harmonic responsibility division on a longer time scale is realized focusing on long-term stable operation of power grids. I PCC are harmonic voltage and harmonic current, respectively, at the PCC;

Projection Coefficient Calculation
. I u ,Z u , .
I u represent the system side equivalent harmonic current of the harmonic source, equivalent harmonic impedance, and branch harmonic current, respectively; .
I u , represent the counterparts on the user side [18]. As proved in [32], it is feasible to use the actual measured current . I instead of the theoretical current . I to carry out the responsibility assessment of centralized multi-harmonic sources, and most current studies default to do so. According to the superposition theorem, . U PCC can be easily calculated as where n is the number of major harmonic sources connected to the PCC; . U k is the harmonic voltage at the PCC generated by the harmonic source k; Z k is the equivalent harmonic impedance excluding the branch k; . U 0 is the background harmonic voltage for other components at the PCC, including the system harmonic voltage and the user-side nonmajor harmonic source. Assuming three main harmonic sources (n = 3), the relationship between the harmonic source voltage for each user and the harmonic voltage vector at the PCC is shown in Figure 2.
. U PCC and . I k are the modulus measured at the PCC. h k (k = 1, 2, · · · , n) is the projection coefficient, calculated by the linear regression method.
Since all the terms in Equation (4) are scalars, hk can be calculated from linear regression, which avoids solving vector equations. Note that the solution of Equation (4) will have high accuracy if U 0 cos α and hk are both constant. However, U 0 cos α varies in real systems, which may cause large errors when using the linear regression method.

Principle of Data Selection
As it is shown in Figure 1, where the system side and the user side are considered as two different parts, the harmonic voltage Equation at the PCC is expressed as Since Z u << Z c , and Z c is the integrated equivalent harmonic impedance of the user side. Then, we have: Based on Equations (5) and (6), the background harmonic voltage of the system side can be expressed as The phasor . U 0 can be calculated from the vector Equation (7), and then the data for linear regression are selected via the cluster analysis [21].
The cluster analysis methods for data selection have two major drawbacks: (i) vector calculation is required, and (ii) these methods only select the data with small background harmonic voltage fluctuations from the system side. When the user side contains unknown harmonic sources that fluctuate, these methods consider them to be constant, resulting in errors that are unavoidable.
The data with small background harmonic voltage fluctuations can be directly selected using the data correlation methods; the linear equation is then solved according to the selected data, and the projection coefficient is calculated. These methods require no estimation of the system side harmonic voltage, and they also minimize the influence of the unknown harmonic sources on the user side.
Although the correlation between different sequences has been widely studied in the existing data correlation methods, in this paper we propose a new data selecting method based on the CCA. Equation (4) shows that when U 0 cos α is a constant, U PCC is linearly related to I k (k = 1, 2, . . . , n). This linear relationship between U PCC and I k (k = 1, 2, . . . , n) is affected by the term U 0 cos α if it fluctuates. In other words, as the fluctuation of U 0 cos α increases, the linearity between U PCC and the linear combination of I k (k = 1, 2, . . . , n) decreases, and vice versa. The data can thereby be selected according to the degree of linearity between U PCC and the linear combination of I k (k = 1, 2, . . . , n). The higher the degree of linearity is, the smaller the fluctuation of U 0 cos α is, and vice versa. Since the traditional correlation analysis methods, such as the Pearson correlation coefficient evaluation method, cannot be applied to cases with two groups of multiple random variables, in this paper we use the CCA method to analyze the correlations between the harmonic voltage and the harmonic currents for multiple harmonic sources.

The Method of CCA Data Selection
The analysis in the previous section indicates that the background voltage is considered constant when the harmonic voltage and the harmonic currents are highly relevant, and thus a multivariate correlation analysis is need to analyze the relationship between the harmonic voltage and the harmonic currents. The CCA is such a method used to study the correlation among multiple sets of data. The two sets of variables studied are X = (x 1 ,x 2 , . . . ,x n ) and Y = (y 1 ,y 2 , . . . ,y n ), where x i and y i are both m-dimensional vectors, and the two groups of variables are standardized. X' and Y' are defined as comprehensive variables where X' and Y' are both m-dimensional vectors; these two comprehensive variables are to replace the original variables according to the optimization method. More details about the optimization method can be found in [33].
The correlation between the input parameters X and Y can be found by calculating the Pearson correlation coefficient. The expression is shown in Equation (8) as where comprehensive variable X and Y , X = X, Y is the optimal linear combination of the variables in Y, and the principle is to make the correlation of X and Y strongest; cov(X , Y ) represents the covariance between X and Y ; D(X ) and D(Y ) are the variances of X and Y , respectively. If the correlation between X and Y satisfies the criterion in Equation (9), they are considered strongly correlated [34].
In order to illustrate the effects of the CCA on data selection, we set two harmonic sources on the user side and set the harmonic source on the system side fluctuating. Figure 3 is a schematic diagram of the data selection results, where U PCC is the harmonic voltage at the PCC, I 1 and I 2 represent the harmonic current of harmonic sources 1 and 2, respectively. For the plots on left of Figure 3, the data variance is 0.1, and the background harmonic voltage fluctuations are considered small; while the plots on the right show relatively large background harmonic voltage fluctuations with a variance of 0.5. When the background harmonic voltage fluctuations are small, there is no obvious similarity between the harmonic currents of the two harmonic sources and the harmonic voltage at the PCC, and the comprehensive variable have a strong similarity with the harmonic voltage. From the plots on the right of Figure 3, there is no obvious correlation between the comprehensive variable and the harmonic voltage at the PCC, which is in agreement with the theoretical analysis in Section 2 that as the background harmonic voltage fluctuations increase, the correlation between the comprehensive variable and the PCC harmonic voltage decreases. Based on the CCA, the data with small background harmonic voltage fluctuations can be screened for subsequent analysis and calculations.

Harmonic Data Sliding Window Analysis
Since it is almost impossible to directly obtain the correlation of the data at a certain moment, the analysis needs to segment the data first. The CCA method proposed in this paper based on sliding window analysis. Assume that the length of each data set is m, and one sliding window includes t sets of harmonic voltage and current amplitude data.
The CCA is performed for each data segment. r j is the correlation coefficient of the jth data segment. When r j ≥ r th (r th is the correlation coefficient threshold), the background harmonic voltage fluctuations of the jth segment data are considered small, hence each harmonic current and the harmonic voltage at the PCC in the jth segment can be used to calculate the harmonic responsibility coefficient; when r j < r th , the jth data segment is not available, thus is discarded.
The CCA method based on sliding window analyzes long-term data in stages by intercepting data segments. This method analyzes the correlation between the harmonic voltage and the harmonic currents in each data segment. On this basis, the data segments with smaller background harmonic voltage fluctuations are selected, which can be used to calculate the projection coefficient by linear regression; the harmonic responsibility index can then be calculated from Equation (2). Figure 4 depicts the procedures of the harmonic responsibility division method for multi-harmonic sources, that are based on the CCA proposed in this paper. A brief description of the procedures is provided as follows. Step 1: Measure the harmonic data at the PCC for n major harmonic sources during the time periods of interest. The measurements include the harmonic voltage data at the PCC, U PCC , and the harmonic current data for each source I k (k = 1, 2, · · · n).

Harmonic Responsibility Division Procedures
Step 2: Determine the length of the sliding window t by considering the accuracy and time of the calculation. The smaller sliding window t, the higher calculation accuracy and the longer calculation time.
Step 3: Calculate, using the CCA method, the correlation coefficient, r j (1 ≤ j ≤ n), between the harmonic voltage data and the harmonic current data of the jth window. The data are selected if r j ≥ r th , and are discarded otherwise.
Step 4: Solve Equation (4) using the linear regression method to calculate the projection coefficient.
Step 5: Calculate the harmonic responsibility index for each harmonic source from Equation (2).
Step 6: Move j one point backwards, and repeat steps (3) to (5) until j = n.
Step 7: Calculate the average of the harmonic responsibility index in each period to get the total harmonic responsibility index for each user, the expression is shown in Equation (10) as where, n is the number of total data segments after screening in the analysis period.

Simulation Specifications
The case study demonstrated here is the Norton equivalent circuit with three harmonic sources (see Figure 1 with n = 3). The simulation is conducted in Power System CAD (PSCAD), and we take the fifth harmonic as an example. The parameters of the system side and three harmonic sources (user side) are listed in Table 2. The harmonic impedance is a series connection between the resistance R and the reactance L. The fluctuation of each harmonic source is assumed to follow normal distribution with zero mean value, and a variance shown in the last column of Table 2. Note that we used four different variance values to represent varying degrees of fluctuation on the system side. According to the previous analysis, we set r th to 0.7. Five different methods are cited to divide the harmonic responsibility for each user, which are the partial least squares method [25] (method 1), the mean shift method [35] (method 2), M-estimation robust regression method [20] (method 3), multiple linear regression method [19] (method 4) and the method proposed in this paper. The simulations obtain the harmonic voltage and the harmonic current of each user branch at the PCC.

Background Harmonic Voltage Fluctuations
The fluctuation of each harmonic source is specified according to Table 2. Here we take the variance of 0.15 as an example. Figure 5 shows the selected data and the removed data of the harmonic voltage and the harmonic currents. The background harmonic source is the system harmonic source, and its degree of fluctuation is larger than any other harmonic source at the PCC. Thus, it can be seen from Figure 5, the data with small system harmonic source fluctuations are selected by the CCA, and the data with large system harmonic fluctuations are removed. Based on the selected data, the projection coefficient is computed from the linear regression method; the harmonic responsibility index then can be calculated from Equation (2). Table 3 shows the harmonic responsibility index calculated from method 1, method 2, and our proposed method. It should be noted that both method 1 and method 2 require phasor information, while our proposed method does not. The actual value in Table 3 represents the ratio of the projection of each user's harmonic voltage on the harmonic voltage at the PCC to the PPC's harmonic voltage when the harmonic source has no fluctuation. It can be seen from Table 3 that the proposed method has higher accuracy than method 1, and its error is similar to that of method 2. Since the proposed method requires neither the estimation of the background harmonic voltage nor the phase information, it has higher practical values in engineering applications.
With the background harmonic voltage fluctuation variance set to 0.15, the results from method 3, method 4, and our proposed method are shown in Table 4. Note that neither method 3 nor method 4 requires phasor information. From Table 4, the accuracy of the proposed method is comparable to that of method 3, and is higher than that of method 4. Method 3 used the iterative solution of least squares after weight reduction to weaken the impact of abnormal data on linear regression solutions, but it needs a regression calculation for each harmonic source. Compared with method 3, the method in this paper requires no iteration, only needs one screening and regression to obtain higher accuracy results, and the algorithm is more convenient. The premise of method 4 is that the background harmonic voltage does not fluctuate, so when the background harmonic voltage fluctuates, the method of this paper is more accurate.
Data utilization rate, W, is defined as W = number of selected data total number of data × 100% Under different background voltage fluctuations, the data utilization rates of the method proposed in this paper are shown in Table 5: It can be seen from Table 5 that as the background harmonic voltage fluctuations increase, the number of the selected data for linear regression decreases. This can be explained as: with the increase of fluctuation degree, the less available data can be screened out by the CCA method. When the total amount of data is the same, the reduction of filtered data will cause the accuracy of the linear regression method to decrease.
The analysis time scale is another factor affecting the results of the harmonic responsibility index. Table 6 shows the harmonic responsibility index calculated at different analysis time scales. Since the harmonic responsibility fluctuates within a small range during the period of analysis, its actual value should be the same as the actual value in Table 3. As is indicated in Table 6, the accuracy of the harmonic responsibility index improves with the increase of the number of data. As the number of data increases, the amount of data selected for linear regression also increases, resulting in more accurate linear regression solutions, and harmonic responsibility indicators. The accuracy of the responsibility division indicator reaches a peak as the number of data is large enough; further increasing the number of data would not significantly increase the accuracy.

Systems Containing Unknown Harmonic Sources
In the system shown in Figure 1  As can be seen from Figure 6, the proposed method can still achieve high accuracy at the presence of an unknown harmonic source. In addition, as the fluctuations of the unknown harmonic source increase, the accuracy is still higher than either the partial least squares method or the mean shift method.
The proposed method selects the data with strong correlations between the harmonic voltage and the harmonic currents for linear regression. Regardless of where the harmonic voltage fluctuations are from (system side or unknown harmonic source), the proposed method is able to solve the harmonic responsibility index using the linear regression method. The partial least squares method takes all the data into account for linear regression which leads to large errors in the calculation of the responsibility index. The mean shift method estimates the background harmonic voltage vector for the system side and selects the system side background harmonic voltage data with small fluctuations for linear regression. However, if an unknown harmonic source exists, it is impossible to select the data with small fluctuations for the unknown harmonic source. In solving this linear regression problem, the voltage of the unknown harmonic source is considered as a part of the background harmonic voltage. The fluctuations of the unknown harmonic sources reduce the linearity between the harmonic voltage at the PCC and the harmonic currents for all sources, resulting in a linear regression solution with low accuracy. When the fluctuation scale of an unknown harmonic source increases, the accuracy of the harmonic responsibility index decreases.

Changes in Harmonic Responsibility during the Period of Analysis
The above analysis is based on the assumption that the user's harmonic responsibility is relatively stable around a certain central value. In this section we discuss the scenarios when the harmonic responsibility fluctuates significantly. The parameters of each harmonic source are specified as in Table 2. The harmonic source of User 1 is specified to vary from 3.03∠-7.3 • to 0.71∠-7.3 • at a certain moment, the background harmonic voltage fluctuation is set to 0.1, and the window width t is 60. The harmonic responsibility index in each window width is shown as a function of time in Figure 7. The black curve in Figure 7 represents the harmonic responsibility index calculated for User 1 within each window width, while the red curve represents the average value of User 1's harmonic responsibility index. It is obvious that the harmonic responsibility has changed significantly, and hence we separate the black curve into two different periods, namely period 1 (left) and period 2 (right). The average value of the harmonic responsibility index for each period is calculated to characterize the harmonic responsibility during that period. The actual values and the average values of the harmonic responsibility indicators are shown in Table 7. According to the data utilization defined in Equation (11), the total data utilization rate is W = 78.86%. The background harmonic voltage fluctuates little in this time period, hence the data utilization rate is relatively high.
As can be seen from Table 7, the method proposed in this paper can still remain a high accuracy under the variance of responsibility. Since filtering data in the proposed method is based on correlation analysis, the selected data have small fluctuations in background harmonic voltage, thus have high accuracy when conducting linear regression. When the harmonic responsibility of the data remains the same in the period of analysis, the harmonic responsibility of the unselected data is approximately the same as that of the selected data; when the harmonic responsibility of the data varies during the period of analysis, the selected data will be the same as in the case when the harmonic responsibility data is unchanged, because the correlation between the harmonic voltage and the harmonic current is unchanged. Since this paper uses the sliding window-based responsibility division, the harmonic responsibility varies significantly with time. To improve the accuracy of the calculation, the average harmonic responsibility of different sections is solved separately.
In the method proposed in this paper, when the background harmonic voltage fluctuates greatly and the harmonic responsibility varies drastically, such data will not be selected, so in this case, a lack of responsibility will occur.

Practical Case Analysis
In this section, we use a practical case, based on the harmonic data of a 220 kV substation in Zhangzhou, Fujian Province, China, for demonstration. Figure 8 shows the diagram of a substation, and a 110 kV bus with four feeders. The power quality monitoring device measures the 95% probability value of harmonic voltage and the currents at the PCC every 3 min. (In the Chinese national standard, the 95% probability value can be determined by the following method: rearrange the measured data from big to small, then 5% data from the bigger part and select the maximum one from residual part to be 95% probability value [29]). The data used for analysis were collected within 7 days (from 20 September 2019 to 27 September 2019).  Figure 9 shows the harmonic voltage and the current data for the fifth harmonic. Since the phase angle data are unavailable in this case, the mean shift method and other methods that require the solution of the phasor equations cannot be applied. In actual systems, the background harmonic voltage is usually under fluctuations, and therefore, the accuracy of the harmonic responsibility index calculated from the linear regression method is low.
The harmonic responsibility division results of the proposed method are shown in Table 8, from which the Manan II Road should take the main harmonic responsibility, as its harmonic responsibility index is the most among the four roads. This point has been proved by the on-site investigation and the actual measurements that the main load of Manan II Road is the electric railway, which results in a large amount of the fifth harmonic emission. The loads for the rest three roads are mainly from the residential electricity or ordinary industries, so their fifth harmonic emissions are small.

Conclusions
This paper proposes a sliding window based canonical correlation analysis method that can select out appreciated data even when the background harmonics fluctuate greatly, and then divides the multi-harmonic responsibility by partial least square regression method without relying on phase angle information.
In the case study, the following observations are made: 1. Compared with applying traditional partial least squares algorithm and the mean shift method in the harmonic responsibility division, the proposed method shows higher accuracy and restricts the error within 10% in the studied cases.

2.
Compared with the M estimation robust regression method, the proposed method shows a similar accuracy that is within the allowable error range of 8%. However, the limited amount of selected data from the M estimation robust regression method in the practical cases will be difficult to reflect the overall situation of harmonic responsibility.

3.
Merely using available data from the online power quality monitoring system and focusing on long-term stable operation of power grids, the proposed method realizes harmonic responsibility division on a longer time scale, weakens the impact of shortterm measurement and analysis errors and is applicable to practical projects.
As for future work, methodologies for parameter adjustment can be developed to further improve the accuracy of responsibility division for different types of customers, especially users with severe background harmonic fluctuations.
Author Contributions: Conceptualization, Y.Z. and Y.W.; methodology, Y.Z. and Z.S.; software, Y.W. and J.G.; validation, Y.W. and J.G.; writing-original draft preparation, Y.W.; writing-review and editing, J.G., Y.Z. and Y.W.; supervision, Y.Z. and Z.S. All authors have read and agreed to the published version of the manuscript.
Funding: This research is supported in part by the National Natural Science Foundation of China (51777035).
Data Availability Statement: Some or all data, models, or codes generated or used during the study are available from the corresponding author by request.