Harmonic Propagation and Interaction Evaluation between Small-Scale Wind Farms and Nonlinear Loads

Distributed generation is a flexible and effective way to utilize renewable energy. The dispersed generators are quite close to the load, and pose some power quality problems such as harmonic current emissions. This paper focuses on the harmonic propagation and interaction between a small-scale wind farm and nonlinear loads in the distribution grid. Firstly, by setting the wind turbines as P – Q ( V ) nodes, the paper discusses the expanding Newton-Raphson power flow method for the wind farm. Then the generalized gamma mixture models are proposed to study the non-characteristic harmonic propagation of the wind farm, which are based on Gaussian mixture models, improved phasor clustering and generalized Gamma models. After the integration of the small-scale wind farm, harmonic emissions of nonlinear loads will become random and fluctuating due to the non-stationary wind power. Furthermore, in this paper the harmonic coupled admittance matrix model of nonlinear loads combined with a wind farm is deduced by rigorous formulas. Then the harmonic propagation and interaction between a real wind farm and nonlinear loads are analyzed by the harmonic coupled admittance matrix and generalized gamma mixture models. Finally, the proposed models and methods are verified through the corresponding simulation models in MATLAB/SIMULINK and PSCAD/EMTDC.


Introduction
Following the depletion of fossil energy, renewable energy generation is a fast growing area nowadays.With the increasing installed wind power generation capacity, we need to find more effective ways to utilize wind energy.
Different for large-scale wind farms integrated into the regular high voltage transmission line [1,2], dispersed wind generators are more flexible.Small-scale wind farms usually consist of several dispersed wind generators, which are directly connected to a medium or low voltage distribution grid [3].This allows the effective utilization of the limited wind resource in cities or suburbs and shortens the transmission distance.The small-scale wind farm is however quite close to the load, which can pose some power quality problems such as voltage fluctuation and harmonic current emissions.
The system impedance in the distribution grid is highly resistive, and the voltage fluctuation at the Point of Common Coupling (PCC) is strongly affected by wind power smoothing level [4,5].Wind power data are non-stationary random signals because of the variable wind speed and wind energy conversion systems [6][7][8].At the same time, the fluctuating fundamental power can affect the controller of conversion systems and induce harmonics, so the harmonic current emissions of wind generators are analyzed by stochastic assessment.New models and methods for harmonic propagation and interaction issues can be verified in several simulation systems, such as MATLAB/SIMULINK, PSCAD/EMTDC and HYPERSIM [9][10][11].
Harmonic current emissions of wind generators consist of characteristic and non-characteristic harmonics.Characteristic harmonics are affected by the topology structure of converters and the switch control strategy of converter valves, whereas non-characteristic harmonics are induced by the inherent slot harmonics of wind turbines, nonsinusoidal distributions of the stator and rotor windings, and operating conditions of converters and the outside power grids [12,13].Non-characteristic harmonics of different wind turbines are almost independent low-order harmonics, even for the doubly fed induction generator and direct-drive permanent magnet generator.
In order to analyze the complex harmonic current emissions, Generalized Gamma Mixture Models (GGMM) are proposed in this paper.GGMM consists of Gaussian mixture models, improved phasor clustering and generalized Gamma models.GGMM divide the complicated probability distribution functions (pdfs) of harmonic sources for the harmonic propagation and interaction study.On the other hand, the port characteristics of nonlinear devices in the distribution gird become variable due to the installed wind generation.Power electronic devices are normally nonlinear [14,15], and the AC/DC converter is important for them.The DC side part can be a linear load, battery, and motor drives with DC/AC inverters [16,17].These nonlinear loads generate several kinds of harmonic emissions into the power grid, and iterative harmonic analysis is the common way to study harmonic power flow [18].In order to analyze the harmonic power flow, accurate harmonic source modeling is necessary.The predetermined and harmonic-voltage dependent current source models are both defective [14,18].Reference [15] proposes a harmonic admittance matrix to reflect the relationship between the AC side harmonic voltages and currents of a rectifier.The harmonic coupled admittance matrix model is strictly deduced and adjusted for the harmonic interaction in this paper.
Therefore this paper focuses on the harmonic propagation and interaction between a small-scale wind farm and nonlinear loads in the distribution grid.Harmonic emissions of nonlinear loads become random and fluctuant due to the non-stationary wind power, meanwhile wind harmonic emissions make it more complex.In this paper, first the power flow expanding model of wind turbines based on P − Q(V) characteristic is studied, and GGMM of harmonic current are proposed and described in detail.The switching harmonic currents produced by the grid-side and rotor-side converters of are also analyzed in depth.Then the harmonic admittance matrix models of nonlinear loads combined with the wind farm are established to analyze the harmonic propagation and interactions.Finally, the proposed models and methods are verified through a simulation model in MATLAB/SIMULINK and PSCAD/EMTDC.

Power Flow Modeling of Wind Farm
When a small-scale wind farm is installed in the distribution grid, the port characteristics of nonlinear loads will become more complex.Compared to the variable loads, the small-scale wind farm induces a fluctuating negative power flow, and supplies voltage to the distribution grid.
The most common forms of wind generation contain asynchronous generators (AGs), doubly-fed induction generators (DFIGs) and direct-drive generators.In this paper, the DFIG is modeled to study the harmonic propagation in the distribution grid.The topological diagram of the small-scale wind farm is shown in Figure 1.
The steady state circuit model of DFIG is shown in Figure 2 [19].The term ω e is the nominal angular speed, ω s , ω r , ω m are the stator, rotor, and rotating angular speed, and N = ±1 represents the positive and negative circuit.When the fundamental power flow is studied,  The reactive power of DFIG is affected by the bus voltage, which is different from the normal PQ or PV node in the power flow study, so the expanding Newton-Raphson power flow method is utilized to study the fundamental power flow, in which the wind generations are set as P − Q(V) nodes [20,21].The Q − V relationship is given by Equation (1), and the output active power is affected by wind speed: where ss s m X X X = + ; Define P is as the active power, P − Q(V) nodes are used to represent the wind power nodes, and their corresponding power equations in the expanding Newton-Raphson models are given by Equation ( 2): where:

Harmonic Source Modeling of Wind Farm
In order to analyze the harmonic emissions of DFIG more accurately, a detailed model is established in MATLAB/SIMULINK.The back to back converter is typically based on a two-level voltage source.The electrical control strategies of the rotor and grid side converters are based on reference [22].Pulse Width Modulation (PWM) is used to reduce the harmonic emissions and the topological diagram of DFIG is shown in Figure 3.
The harmonic emissions of DFIG mainly consist of three parts: (i) the inherent harmonics; (ii) the switching harmonics; and (iii), the harmonics induced by the unbalance conditions or the background harmonic voltage.The inherent harmonics are from a mechanism such as the induction machine.Slot harmonics produced by variation in the reluctance due to the slots are typical inherent harmonics [23,24].Nonsinusoidal distributions of the stator and rotor windings generate the Magnetic Motive Force (MMF) space harmonics [11,24,25].
The switching harmonics are produced by the rotor and grid side converters, and different switching strategies and frequencies cause different harmonic emission levels [26][27][28].
The disturbed operation conditions also induce harmonic DFIG emissions, such as the background harmonic voltage, the nonsinusoidal rotor injection, the unbalanced stator load scenario, and the unbalanced grid-connected operation [10,29,30].
The harmonics emissions named non-characteristic harmonics and caused by nonsinusoidal distribution and disturbed operation conditions, are nonlinear and variable.Papathanassiou, Tentzerakis and Sainz studied the harmonic emissions of practical wind turbines in references [12,13,31], and drew two important conclusions: (i) the influence of a wind farm operating point on harmonic current is small, and the pdfs of h th harmonic currents do not change significantly from one power level to another; (ii) the fifth and seventh harmonics are the dominant components, and the total harmonic current distortion mainly depends on them, so the low order harmonics are the main components of non-characteristic harmonics.In this paper the fifth and seventh harmonics are researched in depth, and the proposed models are also valid for other low order harmonic emissions.This paper focuses on the evaluation of the non-characteristic and switching harmonic emissions of a small-scale wind farm.The orders of the non-characteristic harmonics are the fifth and seventh, and the order of the switching harmonic is affected by the switching frequency.

Low Order Non-Characteristic Harmonics
When the fifth harmonic current is analyzed, the steady state circuit of DFIG is in negative sequence.The parameters in Figure 2 are calculated as follows: When the seventh harmonic current is analyzed, the steady state circuit of DFIG is in positive sequence.The parameters in Figure 2 are calculated as follows: where:

( ) ( cos , sin )
Traditional Rayleigh distribution, Bessel function, and Rician distribution are all derived from Equation (4), in the hypothesis that the pdfs of X h and Y h can be approximated by Gaussian distributions [32].The hypothesis is important to study the harmonic propagation and summation [33].
For a large-scale wind farm, the number of wind power generators is large enough to ensure the Gaussian hypotheses of X h and Y h based on the central limit theorem, but for a small-scale wind farm, the number of wind power generators is limited, and the probability distribution of the non-characteristic low order harmonic may not be the normal distribution [13,31].
If the probability distribution function (pdf) of the harmonic emission becomes complex, the harmonic propagation and interaction analysis will be difficult.In this paper, Generalized Gamma Mixture Models are proposed to study the probability distributions of non-characteristic harmonics.And GGMM consist of the Gaussian mixture models, improved phasor clustering and generalized Gamma models.
(i) Gaussian mixture models Gaussian mixture models (GMMs) can accurately approach the non-negative Riemann integral functions such as different pdfs by dividing the complicated pdfs into several simplified Gaussian distributions [34].GMMs are fitted with Equation (5): where:  θ ε ε μ μ σ σ = are the parameters of GMM, and can be estimated by expectation maximization (EM) method [35,36].
(ii) Improved phasor clustering models After step (i), the pdfs of X h and Y h can be divided into several Gaussian components: But the pdf of harmonic phasor amplitude is determined by X h , Y h and their correlation coefficient.
, the pdf of harmonic phasor amplitude can be given by Equation ( 7): 7) is not adequate.The optimal phasor clustering method which can obtain the inherent characteristic of harmonic phasors is proposed to solve the correlation problem.
Reference [3] divides the harmonic phasors based on the shortest distance away from the central phasors.In this paper, the harmonic phasors are clustered by the whole distribution characteristic.Firstly, the central phasors are defined as ( xm μ , yn μ ) based on Cartesian axes, and their dimension is m n × .xm μ and yn μ are the mean values of X h and Y h based on step (i).
Then the vector space is divided into several regions according to the statistical properties of the Gaussian components in step (i).If the division rule is the shortest distance away from the central phasors, the variance characteristic will be ignored, so the optimal phasor clustering method divides the vector space by the intersection among the Gaussian components.
For example, the pdf of a phasor's X component is combined by two Gaussian functions, N(0,1) and N (10,5).The intersection point is (2.3, 0.025), as shown in Figure 4a.x = 2.3 is set as the first division axis.The pdf of Y component is just the same, and y = 2.3 is the second division axis.x = 5 and y = 5 are the division axis under the shortest distance rule.The harmonic phasors between x = 2.3 and x= 5 should belong to the N (10,5) part, so the shortest distance rule is not adequate.Furthermore, the divided regions are presented in Figure 4b.The optimal phasor clustering method contains the whole probability property of harmonic phasors.
After phasor clustering, the pdfs of X h and Y h in each cluster can be approximated by Gaussian distribution.So if ≠ , the pdf of harmonic phasor amplitude based on phasor clustering can be given by Equation ( 8): (iii) Generalized Gamma models After step (i) and (ii), the pdf of harmonic currents can be given by Equations (7,8).And the complicated pdf is divided into the ideal functions like f m,n (x + jy) or f k (x + jy) which can be approximated by GGD due to the pdfs of X h and Y h in each cluster agree with Gaussian distributions.
GGD with two parameters is given by Equation ( 9): (

( ) exp[ ( / ) ] ( )
where: (2 ) According to the above three steps, any kind of pdfs of harmonic phasors can be approximated by GGMM.And if the pdfs of X and Y are just fitted with Gaussian distribution, GGMM will be directly modeled by step (iii).

High Order Switching Harmonics
High order switching harmonics of DFIG with PWM controllers consist of two parts: the harmonic emission of the grid-side converter is directly injected into the grid system; the harmonic produced by the rotor-side converter intertwine through the electromagnetic coupling between the stator and rotor, and then is injected into the power grid through the stator.
The switching frequencies of the rotor-side and grid-side converters are set as f rc and f gc .And the harmonic voltage of the grid-side converter is calculated by Equation (10) [37]: where M is the modulation index amplitude; and V dc is the DC side voltage.Normally, V dc is constant by the optimal control strategy of the grid-side converter.Then the amplitude of the output harmonic voltage is determined by M.
As the fundamental frequency of the grid-side converter is f e = 50 Hz, the switching harmonic orders of the grid-side converter contain:  2 gc e f f ± is the only harmonic order we need to consider.Other switching frequencies are ignored due to their small amplitude.
On the other hand, the control strategy of DFIG is constant frequency with variable speed.The reference frequency of the rotor-side converter is| | e sf .Then the switching harmonic emissions of the rotor-side converter are injected into the grid through the air gap magnetic field between the stator and rotor.Similar to the harmonic emissions of the grid-side converter, the switching harmonic order of the rotor-side converter is According to the steady state of DFIG in Figure 2, Non-characteristic and switching harmonic emissions of DFIG are discussed in this section.In order to study the harmonic propagation, the harmonic source model of DFIG is given in Figure 5.It consists of three parts: induction motor model is based on Figure 2; harmonic source part is the non-characteristic or switching harmonic emissions; auxiliary load is the supply load in the wind farm, and it can also be compensation devices.

Harmonic Admittance Matrix Model of Nonlinear Loads
In this paper, AC/DC converters installed in the distribution grid are the nonlinear harmonic sources.The harmonic coupled admittance matrix models are utilized to analyze the harmonic emissions of AC/DC converters [15].The matrix models are adjusted to calculate the harmonic propagation and integration between the integrated wind farm and nonlinear loads.
The topological diagram of AC/DC converter with the equivalent load is shown in Figure 6.E dc is the equivalent voltage source on the DC side; R dc is the equivalent resistance; L is the DC side smoothing inductance, commutation angle and inductance are μ and L C ; α is the thyristor firing angle; I dc is the DC current mean value; V dc is the DC voltage mean value; V ac is the DC side phase voltage.) 2 ] 4 where U i is the mean value of fundamental phase voltage at bus i , ( i =1,..,5); L ic and μ i are the commutation inductance and angle of the number i converter; and α i is the thyristor firing angle.

Harmonic Power Flow Modeling with Wind Farm and Nonlinear Load
Fundamental voltage phase angle and magnitude of power flow are the basis of the parameter evaluation of the admittance matrix of converters.According to Section 2, the wind generation is set as P − Q(V) node, and the converters are set as PQ nodes in the fundamental power flow. 1 ∠ is calculated from the fundamental power flow, and α i , μ i , E idc are given by Equations ( 11) and (12).Then the harmonic coupled admittance matrix is built based on Equation (13).The other devices in the network are linear, and the node admittance matrix is shown in Equation ( 17): The harmonic power flow results are given by Equation ( 18) which reflects the admittance matrix of the distribution network, nonlinear converters and wind farm: (18) where and I h_wind reflect the harmonic characteristic of nonlinear converters and wind farm.

Introduction of the Simulation Network
In order to verify the proposed models and method, a medium voltage distribution network is built based on the network benchmarks of CIGRE Task Force C6.04.02 [40].The network benchmarks have been used in some European projects [41].In this paper, the average voltage is set as 10 kV, and five converters are connected to the medium voltage network by the 10/0.4transformers, as shown in Figure 7.The simulation system is established in PSCAD/EMTDC.The linear load and line parameters of the simulation system are given in reference [41].A small-scale wind farm is installed at bus 3 to study the harmonic propagation between wind turbines and nonlinear loads.The wind farm consists of 3 × 1.5 MW DFIG wind turbines connected to the medium voltage grid by the 10/0.69kV transformers, and the parameter values of DFIG are given in Table 1.

Harmonic Evaluation of Nonlinear Loads in Different Conditions
The harmonic coupled admittance matrix model in Section 3 is used to evaluate the harmonic emissions of the AC/DC converters' nonlinear load.When there is no wind farm installed in the distribution network, the nonlinear loads are set as PQ nodes in the fundamental power flow.The fundamental voltage value of each nonlinear load is 400 V, and the active and reactive power values of the five nonlinear loads at different buses are given in Table 2.According to Equations ( 11) and ( 12), the parameters of nonlinear loads are calculated, as shown in Table 3.Based on the equivalence parameters of nonlinear loads, the harmonic coupled admittance matrix model is built by Equations (13)(14)(15)(16).According to Equations (17)(18)(19)(20), the harmonic emissions of nonlinear loads in the simulation network are calculated in MATLAB.In order to verify the harmonic coupled admittance matrix model, the harmonic emissions of the nonlinear loads are collected from the time-domain simulation in PSCAD/EMTDC.The comparison results are shown in Figure 8.
The V thd values of each nonlinear load in MATLAB and PSCAD/EMTDC are almost the same, as shown in Figure 8, so the harmonic admittance matrix model established in MATALB can accurately evaluate the harmonic distortion of nonlinear power electronic loads.
Furthermore, when the small-scale wind farm is installed in the distribution grid, the fundamental power flow becomes fluctuant, and the equivalence parameters of nonlinear loads are not constant.
The capacity of the small-scale wind farm is P wind_total = 3 × 1.5 MW = 4.5 MW, and P wind_total is set as the fundamental value of wind farm.Ten initial wind power values are from 0 to 1 p.u. in intervals of 0.1 p.u.The fundamental voltage value of each nonlinear load is 400 V.
The small-scale wind farm is set as P -Q(V) node, and the expanding Newton-Raphson power flow method is applied in the fundamental power flow calculation.And then the parameter evaluations of converters at different buses are given by Equations (11)(12)(13)(14)(15)(16).The final results of the fundamental voltage, parameter evaluation and harmonic emission at bus 12 are shown in Table 4.  From Table 4, it is found that the parameters of the nonlinear load vary with different wind power values.The differences among the thyristor firing angle α values are quite small, and the phase angle of voltage and commutation angle μ change a little.The fluctuant wind power induces the variable fundamental voltage, and then the harmonic emissions of nonlinear loads change, as shown in Table 4.And the V thd value of the nonlinear load at bus 12 changes from 4.2% to 3.7%.

High Order Switching Harmonic Analysis
In order to analysis the switching high order harmonic emissions, a detailed model of the small-scale wind farm is built in MATLAB/SIMULINK.PWM control strategies of the rotor and grid side converters are based on reference [22].The parameter values of DFIGs are given in Table 1.
Due to the sinusoidal distributions of the stator and rotor windings in the detailed model, the switching harmonics are the main components.When the output power of the small-scale wind farm is 4.5 MW, the current emission without filter is shown in Figure 9a.
From Figures 9b,c, the switching harmonic emissions of the grid-side converter are higher than the rotor-side converter.The harmonic orders of the grid-side converter are 42 Hz and 38 Hz.The harmonic orders of the rotor-side converter are 40 ± 1.63 Hz.These harmonic orders conform to the deduction in Section 2.2.Comparing the Figures 9b,c, it is found that the low pass filter reduces the high order switching harmonic effectively.In Figure 9d, the ripple of the DC voltage is apparently repressed by the filter, and then the low order harmonics are reduced.After using the filter, the high order switching harmonic currents are small enough to merge into the grid.Because wind power varies with the random wind speed, the rotor speed and slip are both variable.In Table 5, five different wind speeds are calculated in the model and the relationship between the switching harmonics and slips are studied.It is found that the slip s changes from −0.21 to 1.8 and the switching harmonic orders of the rotor-side converter are changed.Because the switching harmonic currents from the rotor-side converter are too small, only the switching harmonic currents of the grid-side converter are of interest in this paper.In practice, a wind generator may not comprise a neutral conductor or a ground return path, and the transformer windings are delta connected, so the wind farm will not inject zero-sequence harmonic  currents into the grid [42].The 42th switching harmonic currents can't be injected into the grid, so only the 38th switching harmonic currents are analyzed.From Table 5, the 38th switching harmonic current changes from 5.39 to 5.81 A, the difference is small.At the same time, the harmonic emissions of nonlinear loads are low orders, and there is little interaction between the 38th switching harmonic current and nonlinear loads, so the 38th harmonic voltage at bus 3 is directly calculated by the harmonic admittance matrix, V 38thd = 0.7%.Though the amplitude of the high order switching harmonic is not high, it may induce large fluctuant when the resonance occurs at higher frequency.

Non-Characteristic Low Order Harmonic Analysis
Due to the random wind speed and variable environment, the real output power data of the small-scale wind farm are fluctuant and non-stationary.In Figure 10, the active power data of the small-scale wind farm in a month is collected with a ten minute sampling interval.
The scatter plots of the fifth and seventh harmonic current X-Y projections are shown in Figure 11.The reference value is the percentage of rated current I N , and I N = 245 A. Accord to the scatter plots, the pdfs of the X and Y components of the seventh harmonic currents can be fitted with Gaussian distribution and the parameters are   Different from the seventh harmonic currents, the pdfs of X and Y components of the fifth harmonic current can't be approximated by Gaussian distribution, as shown in Figure 11a.And the correlation coefficient of the fifth harmonic current is 5 5  0.746 . The pdf of X components of the fifth harmonic current is verified by the Normal probability plot which can graphically assess whether the data could be fitted with Normal distribution.The pdf is not fitted with a Normal distribution because the plot in Figure 12b is not linear.

(i) GMM of X and Y components
The parameters of GMM are estimated by EM method, and the observed data are the fifth harmonic current emissions of wind generators, as shown in Figure 11a.The GMM order is set as M = 2, and the divided results of the fifth harmonic current are given in Table 6.Table 6.The divided results of the fifth harmonic current.
Table 7.The clustering results of the optimal phasor clustering method.It can be found that the first and second part are the main components, and their correlation coefficients are 1 0.07 r = and 2 0.44 r = − .

Coefficient
(iii) Generalized Gamma distribution Based on GMM and phasor clustering, the observed data of the fifth harmonic currents are clustered into four Gaussian components.Each component can be approximated by GGD, so GGMM of the fifth harmonic current emissions are given by Equation ( 21): The parameters in Equation ( 21) are from Table 7, and the first two parts can accurately approximate the pdfs due to their bigger weighted coefficients.Then the cdf of the fifth harmonic current is approximated by Rayleigh, Weibull, Normal distributions and GGMM, as shown in Figure 13.
The 95% value of the fifth harmonic voltages is defined as an index, CP 95 , which is used to compare the different distributions.In Figure 13, CP 95 of GGMM and the actual data distribution are both 1.5 p.u., but CP 95 of Normal distribution, Weibull and Rayleigh distributions are 1.1 and 1.2 p.u., so the cdf of the fifth harmonic current can be more accurately approximated by the GGMM.

Harmonic Interaction between the Wind Farm and Nonlinear Loads
When a small-scale wind farm is installed in the distribution grid, the fundamental power flow changes with the fluctuating wind power, and then the equivalence parameters of nonlinear loads become variable, as shown in Table 4.At the same time, stochastic harmonic emissions are calculated in the harmonic power flow analysis based on Equations (11)(12)(13)(14).Because the influence of the wind farm operating point on harmonic current is small, the power level and harmonic emissions of the wind farm are analyzed independently.
First, the fluctuating wind power values are divided into ten intervals which are from 0 to 1 p.u. with intervals of 0.1 p.u.In Table 4, the changes of the parameter values and V thd of the nonlinear load in each interval are quite small, so the parameters of nonlinear loads are set as fixed values in each interval.Then the probability value of random wind power in each interval is set as the weighted coefficient.Based on the real wind power data in Figure 10, the weighted coefficient in each interval is calculated and shown in Table 8  In Table 8, the bigger weighted coefficients are in the smaller wind power interval.And the parameter values of nonlinear loads in each interval are from Table 4.For example, α in (0, 0.1 p.u.) is the average value while P wind = 0 and P wind = 0.1 p.u.. Due to the independent relationship between wind harmonic emission and power level, the same pdf of wind harmonic emission is utilized in each interval.According to the results in Section 4.3, the pdfs of the fifth and seventh harmonic current of wind power are both accurately approximated by GGMM.
Then the harmonic interaction between the wind farm and nonlinear loads are analyzed by Equations (17)(18)(19)(20).I s in Equation ( 20) reflects the harmonic characteristic of nonlinear load, and I h_wind represents the wind harmonic current emissions.And Equation ( 20) is changed into Equation ( 22): where represents the wind harmonic propagation.
In each interval the parameters of the harmonic admittance matrix are fixed, and the probability function of _ h wind V  is the linear combination of I h_wind which consists of the fifth and seventh harmonic emissions.
Based on the above complete formula derivation, statistical characteristic of each bus voltage is analyzed and reflected by the generalized gamma mixture models.For example, the statistical results of the nonlinear load at bus 3 in the first interval (0, 0.1 p.u.), in which s = 0.2, are given in the Table 9.According to Table 9, the probability distribution of the fifth harmonic voltage at bus 3 in (0, 0.1 p.u.) is combined by two GGD parts of which the coefficients are ε 1 = 0.8 and ε 2 = 0.2.This is mainly caused by the GGMM of the fifth harmonic current emissions of the wind farm.Different from the fifth harmonic voltage, the characteristics of the seventh harmonic voltage with two coefficients are almost the same.This means that the mixture pdfs of the fifth harmonic current emissions of wind farm mainly affect the fifth harmonic voltage.Furthermore, the standard deviation of the eleventh harmonic voltage is very close to zero which means that the eleventh harmonic voltage is almost constant.
So the fifth and seventh harmonic voltages at bus 3 are analyzed by probability distribution models.The higher order harmonic voltages at bus 3 are nearly constant and directly solved.There are ten intervals for the harmonic propagation and interaction evaluation, and all of them are calculated.The probability characteristics of the fifth and seventh harmonic voltage in (0, 0.1 p.u) and (0.9 p.u., 1.0 p.u.) are compared, in Figure 14.The pdf of the fifth harmonic voltage at bus 3 is approximated by GGMM, and the pdfs in (0, 0.1 p.u) and (0.9 p.u., 1 p.u.) are shown in Figure 14a (real line) and Figure 14a (dashed line).Then the pdf of the seventh harmonic voltage at bus 3 is shown in Figure 14b.The harmonic voltage distribution in the higher interval is closer to zero because the wind power in the higher interval is larger.Larger wind power improves the fundamental voltage, and then the voltage harmonic distortion level is reduced.
The index CP 95 is used to evaluate the cdfs of the harmonic voltages at bus 3, and the cdfs of harmonic voltages are shown in Figure 15.CP 95 of the fifth harmonic voltage in (0, 0.1 p.u) and (0.9 p.u., 1.0 p.u.) are 1.32% and 1.35%.The difference of CP 95 between the two intervals is quite small.Based on the Figure 14a, the main part of the pdf of the fifth harmonic voltage is closer to zero when the interval value is bigger.But the other part between 1.1% and 1.5% changes a little.The different CP 95 values reflect the final harmonic distortion level, and each part of the generalized gamma mixture models is important.
Different from the fifth harmonic voltage, V X and V Y of the seventh harmonic voltage at bus 3 are approximated by one Gaussian distribution.And the pdf of the seventh harmonic voltage is fitted with GGD which is the simplified type of GGMM.CP 95 of the seventh harmonic voltage in (0, 0.1 p.u) and (0.9 p.u., 1.0 p.u.) are 0.88% and 0.94%, as shown in Figure 15b.Similar to the harmonic analysis in (0, 0.1 p.u) and (0.9 p.u., 1.0 p.u.), pdfs and cdfs of harmonic voltage at different buses are analyzed.Then the final pdf and cdf of the harmonic voltage are the linear combination of ten intervals with different weighted coefficients.The comparison of the cdfs of the harmonic voltages at bus 3 between GGMM and simulation results in PSCAD/EMTDC is shown in Figure 16.In Figure 16a, CP 95 of the fifth harmonic voltage based on the simulation results is 1.29%, and CP 95 based on the proposed GGMM combined with the harmonic admittance matrix model is 1.33%.The difference between the proposed model and simulation result is quite small.Then CP 95 of the seventh harmonic voltage based on the proposed models and simulation results are 0.93% and 0.94%.If the pdf of the fifth harmonic current is directly approximated by Gaussian distribution, the CP 95 values of the fifth and seventh harmonic voltage are 1.18% and 0.88%.The error is quite big, and it can't guarantee the high power quality level.Compared with the Gaussian distribution, the GGMM combined with the harmonic admittance matrix model simulate and evaluate the propagation and interaction between the small-scale wind farm and nonlinear loads more accurately.

Conclusions
The integration of a small-scale wind farm in the distribution grid brings about some challenges to the system safety and reliability.Non-stationary wind power induces a fluctuating fundamental power flow, and the non-characteristic harmonic currents of the wind farm interact with the harmonic emissions of nonlinear loads.In this paper, the fundamental power flow is calculated by the expanding Newton-Raphson method due to the P -Q(V) characteristic of wind turbines.Generalized gamma mixture models are proposed to approximate the pdf of harmonic emissions, which are more effective than general distributions such as the Gaussian, Rayleigh and Weibull ones.Then the voltage harmonic distortion level of nonlinear loads combined with the wind farm is accurately evaluated by the harmonic admittance matrix model.The harmonic interaction between the wind farm and nonlinear loads are calculated by the proposed models.The final results in MATLAB and PSCAD/EMTDC reveal that the generalized gamma mixture and harmonic admittance matrix models accurately evaluate the harmonic propagation and interaction between the small-scale wind farm and nonlinear loads in the distribution grid.

Figure 1 .
Figure 1.The topological diagram of the small-scale wind farm.

Figure 2 .
Figure 2. The steady state circuit model of DFIG.
= ; cosϕ is the power factor, s V  represents the stator voltage; s I  represents the stator current; r V ′  represents the rotator voltage; r I  represents the rotator current; r s and X s represent the stator resistance and reactance respectively; r r ′ and X r represent the rotator resistance and reactance respectively; X m is the excitation reactance; slip s = is the slip ratio.

σ
are the standard deviation values.The pdfs of I h and h φ are given by Equation (4 is the order of GMM; m ε is the weighted coefficient; ( )

Figure 4 .
Figure 4. Verification of the optimal phasor clustering.(a) Pdf and division axis of X component; (b) Comparison of two clustering method.
,39].When the switching frequency is quite high (thousands of Hz), there aren't any background harmonics, the switching harmonic current emissions are affected by r V ′ , and slip.But the switching harmonic current emissions of the grid-side converter are only related to the output harmonic voltage.

Figure 5 .
Figure 5.The diagram of the harmonic source model of DFIG.

Figure 6 .
Figure 6.The circuit diagram of AC/DC converter with the equivalent load.

Figure 7 .
Figure 7.The simulation network with the small-scale wind farm and nonlinear loads.

Figure 8 .
Figure 8.The comparison of V thd values of five nonlinear loads between harmonic coupled admittance matrix model and simulation results.

Figure 9 .
Figure 9. Switching harmonic simulation.(a) Currents of the small-scale wind farm without filter; (b) I thd without filter; (c) I thd with filter; (d) Comparison of V dc with and without filter.

Figure 10 .
Figure 10.Wind power data in 31 days with 10 min intervals.

Figure 11 .
Figure 11.Scatter plots of wind farm harmonic current X-Y projections.(a) The fifth harmonic currents; (b) The seventh harmonic currents.

Figure 12 .
Figure 12.Pdfs of the fifth harmonic current and Normal verification.(a) Pdfs of X and Y components; (b) The Normal probability verification plot.

Figure 13 .
Figure 13.Cdf of the fifth harmonic current approximated by different distributions.

Figure 14 .
Figure 14.Pdfs of the harmonic voltages at bus 3 in different intervals.(a) The fifth harmonic voltage magnitude; (b) The seventh harmonic voltage magnitude.

Figure 15 .Figure 16 .
Figure 15.Cdfs of the harmonic voltages at bus 3 in different intervals.(a) The fifth harmonic voltage magnitude; (b) The seventh harmonic voltage magnitude.

Table 1 .
The parameters of DFIG.

Table 2 .
The power values of five nonlinear loads.

Table 3 .
The parameter values of nonlinear loads without wind farm.

Table 4 .
The parameters and harmonic distortion level of the nonlinear load at bus 12.

Table 5 .
The comparison of two kinds of switching harmonics with different slips. V

Table 8 .
The weighted coefficient in each interval.

Table 9 .
Statistical characteristic of harmonic voltage at bus 3.