CHAPTER(chapter 6)
Models for Representation and Summation of Random Harmonic Currents
R. Langella, A. Testa – Second University of Naples (Italy)
The chapter deals with the probabilistic modeling of harmonic currents. After some definitions about single current injection models, the vectorial summation of random harmonic currents characterized by distributions independent on the time is considered. The basic methods presented in literature for harmonic current summation modeling are descried: Convolution Method, the Joint Density Method and Monte Carlo Simulations. Applications of the methods to different case-studies give an idea of the accuracy of the models considered.
1.16.1Introduction
One aspect of assessing the harmonics tolerated by a power system is the estimation of the statistic figures of harmonics arising from the various time-varying / probabilistic sources.
The assessment of harmonics is not exact or uniform, since there will be unpredictable variations in either the non-linear sources and/or parameters of the system which affect the summation.
The combination of a number of harmonic time-varying sources will generally lead to less than the arithmetic sum of the maximum values due to uncertainty of magnitude and phase angle. Hence the resulting summation is extremely difficult to estimate accurately.
Before introducing the methods that can be adopted to solve the summation problem, it is necessary to introduce and discuss some basic concepts: random harmonic and interharmonic vectors and harmonic summation principle and essential nomenclature.
1.26.2Random Harmonic and Interharmonic Vectors
The analysis of interface among utilities and consumers needs representing the harmonic currents injected from disturbing loads as random vectors. The random behavior of harmonic currents is related to the parameter of influence stochastic nature such as active and reactive powers, network configuration, non-linear load operational conditions, etc.
The harmonic current of order h injected by a non linear load into the network can be represented as a vector , of amplitude Ih and phase h and of Cartesian components Xh and Yh. In any case, the statistical characterization of requires the determination of the joint statistics of a pair of real random variables (Ih,h) or (Xh,Yh). With reference to (Xh,Yh) and omitting the subscript h, the functions to be considered are:
- the joint distribution FXY(x,y), also called joint cumulative probability function (jcpf), that is the probability of the event Xx, Yy;
- the joint probability density function (jpdf), that is, by definition, the function:
;
- the marginal distributions FX(x) and FY(y), also called cumulative probability functions (cpf’s);
- the marginal probability density functions (pdf’s) fX(x) and fY(y).
The graphs in Figure 6.1 show an example of the 5th harmonic time variability magnitude and of its pdf. Figure 6.1a reports the scatter plot of the vectors edges registered and Figure 6.1b the magnitude pdf.
Figure 6.1. 5th harmonic current from measurements: a) Scatter plot, b) Probability histogram of magnitude.
1.36.3Harmonic Summation Principle and Essential Nomenclature
The basis for harmonic combination is the superposition principle. To apply the superposition principle properly, a phasorial composition should be used.
Here, the basic definitions and nomenclature are introduced together with some practical considerations about the computational problems.
1.3.16.3.1Summation of Random Harmonic and Interharmonic Vectors
The problem of modeling the injection of the sum of numerous current vectors in the network needs technical analysis and statistical elaboration of measurements for each load and also knowledge and representation of the way in which random vectors combine during the time to give a resultant.
The sum of N random harmonic vectors gives:
, , .
Great is the practical interest to obtain the pdf of Ih. The correct theoretical approach would be based on the use of the 2N-dimensional joint probability density function
, with .
Equation 6.1
The distribution of Ih is given by:
,
Equation 6.2
being the 2N-dimensional region of the hyperspace z where the constraint , is verified.
Once solved the integral form, it is trivial to obtain the pdf of Ih.
1.3.26.3.2Basic Considerations
The vectorial summation of N random harmonic and interharmonic components, in a defined scenario of space (loads, utility network characteristics, ...) and of time (the year, the annual maximum load day, ...), is in principle very simple and comprehensive if the 2N-dimensional jpdf of the 2N real random variables representing the N vectors involved is available for the scenario at hand.
In practice, also having in mind only numerical approaches that discretize the hyperspace z, assuming M discrete values for each coordinate in its definition interval, the following dramatic computational problems arise:
-the matroid to be utilized for representing the jpdf fZ given in Equation 6.1assumes dimension D=2M2N (f.i. if M=N=10, that are very low values, then D=2x1020);
-the solution of integral Equation 6.2is very time consuming;
-the determination of the jpdf is prohibitive to obtain by both experimental analyses or simulation approaches.
If the dependence amongst the different random vectors is ignored (hypothesis A) or accounted for outside the summation stage, the problem reduces to consider N matroids, each of dimension Di=2M2, to represent N bidimensional jpdf’s (M=N=10 gives D=NDi=2000).
Moreover, if the dependence amongst the pairs of real random variables utilized to represent each vector is ignored (hypothesis B) or accounted for outside the summation stage, the problem reduces to consider 2N matrices, each of dimension Di=2M, to represent 2N marginal pdf’s (M=N=10 gives D=2NDi=400).
Different procedures have been engaged for approaching the summation of stationary random vectors, in order to avoid the Monte-Carlo simulation burden or the dramatic computational problems deriving from the theoretical formulation reported in nomenclature; wide reviews are developed in [4], [5], [6], [8] and [11]. Among these procedures the first part of this section summarizes those widest diffused.
The various methods proposed have been originated from appropriate application of probability theory [18]. Most of them are based on analytical models founded on fully developed convolutions or on the central limit theorem; differences consist on the more or less restrictive hypotheses required by each of them (mainly the statistical dependence or independence amongst the random variables representing the Cartesian coordinates of the vectors and the number of components). All of the methods presented in the following are founded on two general hypotheses:
1) the random vectors are statistically independent;
2) the distributions of harmonic vectors are independent on the time.
1.46.4Basic Methods
The basic methods proposed in literature are the convolution method, the joint density method and finally Monte Carlo simulations.
1.4.16.4.1Convolution method
A first method [4] assumes that the resultant resolved components are statistically independent. The method goes on as follows:
,
,
,
,
,
,
Equation 6.3
where * denotes convolution.
In general the method requires the knowledge of each vector resolved component pdf’s. A simplification is possible when, for a large number N of harmonic vectors, the central limit theorem is applicable to the summation of the Xh,k and of the Yh,k. In this case Sh and Wh pdf’s approximate to two Gaussian distributions of known means and variances. However, the relation Equation 3requires the statistical independence between S2h and W2h.
1.4.26.4.2Joint density method
Another method [5] does not need the independence between Sh and Wh but it requires N so high to make the central limit theorem applicable. Therefore, Sh and Wh are jointly normal with their joint density given by:
,
where r is the correlation coefficient and
.
Equation 6.4
The density function of Ih is directly derived by the following relation:
.
Equation 6.5
When r equals 0, Sh and Wh are independent because they are jointly normal.
In the same hypotheses, the density function of the phase is derivable solving the following integral form:
.
Equation 6.6
1.4.36.4.3Monte Carlo method
Usually, Monte Carlo methods [17][19] are used to simulate a prescribed random behavior of the network loads. That is, random number generators are used to assign a specific probability distributions to certain parameters of the loads, thusly reflecting the random variations in the loads operating condition. In this way, deterministic models of the load can then be used to generate the random harmonic current injection. Subsequently, the statistics of the resulting harmonic voltages are numerically determined, usually from the linear propagation of these harmonic currents through the system impedances. The advantages of this approach is in the possibility of simulating a wide variety of random load characteristics until the resulting statistics agree with available field measurements. The disadvantage is that this method is computationally intensive and time consuming since it is difficult to determine how to adjust the load random models in order to produce desired results. That is, although the method is flexible, the direct relation between the load probability models and the resulting harmonics is not apparent. This means that if a given set of load random characteristics yields unsuccessful results (simulated harmonics do not match available measurements), then limited insight is gained from this simulation trial in terms of altering load models for more accurate results in the next trial. Unless equipped with accurate information on the nature of the random loads, it is extremely difficult to develop simulation models which generate adequate results.
1.56.5Magnitude and phase distributionobtainable starting from the Joint Density Method[16]
In the following some brief remarks on analytical distributions, firstly for the magnitude and then for the phase are recalled. The integral forms (Equation 6.5 and Equation 6.6) are not simply resolvable, due to the complicated structure of fSW. Nevertheless, it is worth noting that if a priori it is possible to assume r equal to zero, then S and W are independent, because they are jointly normal, and the integral forms in (Equation 6.5 and Equation 6.6) become more light to handle. In this case the integration of Equation 6.5 and Equation 6.6 can also be performed by opportune convolutions, as fully shown in [1][2][3], so obtaining a formal but not effective simplification because in the most general case difficulties arise due to numerical problems and to the sensitivity to the real axes partition choice.
The rotation aim is to determine an opportune rotation angle of the original Cartesian reference SW that gives a new Cartesian reference S’W’ in which rS’W’=0. The angle can be obtained as a function of variances and covariances of the original resolved components [3]:
.
Equation 6.7
S and W result two normal and correlated r.v.’s changed by the rotation (Equation 6.7) into S’ and W’, normal and uncorrelated, that is to say, also independent, so simplifying the integral forms (Equation 6.5 and Equation 6.6) and also solving the theoretical problem of the original distribution concerning the assumption of normal jpdf. The outcomes of the rotation are fully shown in [3], also referring to summation methods not utilizing the bivariate normal distribution.
1.5.16.5.1Magnitude
The aim is to have a comprehensive insight into the different fully analytical or empirical solutions, in order to compare the regions of applicability, the performances and, mainly, the usefulness in practical engineering problems in which sometimes only some statistic parameters are available.
1.5.1.16.5.1.1Closed form solutions holding under particular hypotheses
In the hypotheses of rSW = 0 and W = 0, the integration of Equation 6.5provides a complicated solution [6], based on a series of products of Ik, modified Bessel function of integral order k:
,
Equation 6.8
where
...
Equation 6.9
In some special cases, the resultant magnitude density given by Equation 6.8 turns into a simple form. First of all, if S = W = then the Equation 8 becomes:
,
Equation 6.10
the Equation 6.10 is called the Rician probability density function.
Moreover, when it results Si 2, that is the argument of I0(.) is large, the following approximation of Equation 6.9 arises:
,
Equation 6.11
which is a Gaussian law except for the factor (i/S)1/2, and for this reason it is called “almost Gaussian”. It is interesting to note that the previous relations can be utilized also when W0 by the substitution of S with (S2+W2)1/2.
Finally, when S = W = 0 then the Equation 6.9 becomes the Rayleigh distribution:
.
Equation 6.12
As well known, this is the case of random vectors with phases distributed in the whole interval (0,2).
Table 6.1 summarizes the methods for obtaining fI(i) from the bivariate normal distribution fSW(s,w) in closed-form solution holding under particular hypotheses.
Table 6.1 methods for obtaining fi(i) from the bivariate normal distribution fsw(s,w): closed-form solution holding under particular hypotheses
NAME / HYPOTESES / EXPRESSIONBessel / rSW = 0,
W = 0 /
Rician / rSW = 0,
W = 0,
S = W = /
“Almost
Gaussian” / rSW = 0,
W = 0,
S = W = ,
Si 2 /
Rayleigh / rSW = 0,
W = S = 0,
S = W = /
1.5.1.26.5.1.2Empirical solutions
Other methods seem to follow a not fully analytical demonstration for obtaining a closed form for the solution of Equation 6.5. Basically, they impose the equating of some moments of the actual fI with the corresponding moments of an approximate distribution. In spite of all that, they obtain very powerful and accurate results.
In the hypothesis of rSW=0, in [9]a 2 distribution is introduced so obtaining for the magnitude density:
,
Equation 6.13
with the parameters expressed by the following relations:
,
.
Instead, also in the hypothesis of rSW 0, in [15]a method for assimilating the resultant magnitude distribution with a generalized gamma distribution (GGD) is proposed:
,
Equation 6.14
where is the Gamma function and the parameters and can be expressed in terms of the bivariate normal distribution model five parameters:
It is easy to demonstrate that the expression in Equation 6.12 arises from the simplification of the more general expression in (Equation 6.13 and Equation 6.14). Parameters , , and can be derived from a moment fitting procedures as follows. The moment generating function of two jointly normal random variables X and Y is:
,
with their joint moment of order j+k given by:
.
On the other hand, the kth order moment of the GGD can be obtained by means of the following expression
.
Equating the same order moments gives the GGD model parameters.
Table 6.2summarizes the methods for obtaining fI(i) from the bivariate normal distribution fSW(s,w) giving “empirical” solutions.
Table 6.2 methods for obtaining fi(i) from the bivariate normal distribution fsw(s,w): empirical solution
NAME / HYPOTESES / EXPRESSION2 procedure / rSW = 0 /
2 procedure
+ rotation /
Generalized
Gamma
Distribution /
1.5.26.5.2Phase distribution
Attention is paid also to the resultant phase distribution which until now has received little or no importance in literature but which can considerably help in understanding the role played by clusters of harmonic injections at different buses of the network.
The density function of the phase is derivable solving the following integral form:
.
Equation 6.15
The solutions require rSW = 0. Also in this case it is possible to take advantage of an axis rotation. Nevertheless, it is necessary to consider that an angle rotation does change the phase density:
.
1.5.2.16.5.2.1Closed form solution
In the same hypotheses of the Rayleigh distribution for the resultant magnitude, the tandistribution results a Cauchy so that:
,
Equation 6.16
that is the resultant is phase uniform in (0,2).
1.5.2.26.5.2.2Almost analytical solution
In the hypothesis of rSW=0, the solution of the integral form Equation 6.15is obtained in Appendix of [10]and it is given by:
,
Equation 6.17
where the auxiliary variable , which is a function in , is defined by means of the following relation
,
Equation 6.18
and the quantities K, A and B are valuable as shown in Appendix of [16].
It is easy to demonstrate that the expression in Equation 6.16arises from the simplification of the more general expression in Equation 6.17.
1.66.6Case studies
1.6.16.6.1Numerical case-studies
All the above described procedures were tested in different case-studies, also performing Monte-Carlo simulations to obtain a reference.
Case-studies of Rayleigh and Rician magnitude distributions applicability, not reported here for the sake of brevity, were performed obtaining very good results also applying 2 and GGD procedures. Instead, the cases fully reported in the following refer:
- case 1A to completely verified Almost Gaussian applicability conditions, while case 1B and 1C move toward Almost Gaussian non-applicability conditions starting from 1A;
- case 2 to a scenario where the bivariate distribution is characterized by an elliptic symmetry; since the principle axes of the ellipse are not parallel to the Cartesian axes, the correlation coefficient is different from zero (case 2A) but it becomes equal to zero (case 2B) after a /3 angle rotation.
In Table 6.3 the parameter values utilized are reported for the original bivariate distribution, for the GGD and the 2 distribution; for the Monte-Carlo simulation a minimum of 15000 samples has been considered in order to directly solve integral forms of Equation 6.5andEquation 6.15.
Table 6.3 CASE STUDIES PARAMETERS
Bivariate Normal Distribution / GGD / 2 DistributionCase Study / Cartesian Reference / / S / S / W / W / rSW / / / /
1A / SW / - / 20.0 / 1.7 / 0.0 / 1.7 / 0.0 / 34.7 / 20.2 / 69.4 / 2.4
1B / SW / - / 10.0 / 1.7 / 0.0 / 1.7 / 0.0 / 9.3 / 10.3 / 18.5 / 2.4
1C / SW / - / 2.5 / 1.7 / 0.0 / 1.7 / 0.0 / 1.4 / 3.5 / 2.8 / 2.1
2A / SW / - / 3.0 / 1.0 / 5.2 / 1.3 / 0.6 / 4.6 / 6.2 / - / -
2B / S’W’ / /3 / 6.0 / 1.5 / 0.0 / 0.7 / 0.0 / 4.6 / 6.2 / 9.2 / 2.1
From Figure 6.2 to Figure 6.5, the scatter plot, the magnitude and phase pdfs, for the respective cases, are reported.
Figure 6.2shows that Gaussian and Almost Gaussian give results very similar in 1A, similar in 1B and appreciably wide apart in 1C. Concerning the Almost Gaussian it is worth to underline as the right tail better approximates solution, that is to say that higher percentiles are better estimated, due to the better verification of the condition Si 2, especially when this condition is not verified for all the possible values.
It is worth to note as 2 and GGD procedure results are very close to the actual (Monte-Carlo) in all of the cases considered, performing as well as the closed forms, when applicable. Moreover, 2 distribution and GGD procedures have been tested in more general scenarios as those reported in [7], [12] not necessary characterized by scatter plot symmetry, in any case giving good results.
Figure 6.2. Case studies 1A,1B and 1C, current magnitude pdf, obtained by Monte Carlo Simulation (+), 2 procedure (o), GGD model (), Almost Gaussian model () and Gaussian assumption (*).