**2. Sources and problems of harmonics [2]**

Harmonic sources are divided into two categories:


While the established sources of harmonics are still present on the system, the power network is also subjected to new harmonic sources:


The most damaging frequencies to power devices and machines appear to be the lower – below 5-kHz – frequency range. In years past, the magnitudes and sources of these lowerfrequency harmonics were limited and, inmost cases, power systems could tolerate them. The increase in power loss due to harmonics was also neglected because energy costs were low. These conditions no longer apply, and concern for harmonics is now becoming widespread among utilities.

For more Than 100 years, harmonics have been reported to cause operational problems to the power systems. Some of the major effects include:


These effects depend, of course, on the harmonic source, its location on the power system, and the network characteristics that promote propagation of harmonics.

#### **3. Estimation of harmonics and sub-harmonics; the static case**

#### **3.1 Time domain model [3]**

In this model, it is assumed that the waveform under consideration consists of a fundamental frequency component and harmonic components with order of integral multiples of the fundamental frequency. It is also assumed that the frequency is known and constant during the estimation period. Consider a non-sinusoidal voltage given by a Fourier-type equation:

$$w(t) = \sum\_{n=0}^{N} V\_n \sin \left( n a o\_0 t + \phi\_n \right) \tag{1}$$

where

4 Power Quality Harmonics Analysis and Real Measurements Data

A review of the literature indicates that the known sources of harmonics include: 1. Tooth ripple or ripples in the voltage waveform of rotating machines. 2. Variations in air-gap reluctance over synchronous machine pole pitch. 3. Flux distortion in the synchronous machine from sudden load changes.

4. Non-sinusoidal distribution of the flux in the air gap of synchronous machines.

4. Interconnection of wind and solar power converters with distribution systems.

storage batteries, and fuel cells that require dc/ac power converters.

Today's power system harmonic problems can be traced to a number of factors:

which create load-generated harmonics throughout the system.

6. Network nonlinearities from loads such as rectifiers, inverters, welders, arc furnaces,

While the established sources of harmonics are still present on the system, the power

1. Energy conservation measures, such as those for improved motor efficiency and load matching, which employ power semiconductor devices and switching for their operation. These devices often produce irregular voltage and current waveforms that

5. Static var compensators which have largely replaced synchronous condensors as

6. The development and potentially wide use of electric vehicles that require a significant

7. The potential use of direct energy conversion devices, such as magneto-hydrodynamics,

1. The substantial increase of nonlinear loads resulting from new technologies such as silicon-controlled rectifiers (SCRs), power transistors, and microprocessor controls

2. A change in equipment design philosophy. In the past, equipment designs tended to be under-rated or over-designed. Now, in order to be competitive, power devices and equipment are more critically designed and, in the case of iron-core devices, their operating points are more into nonlinear regions. Operation in these regions results in a

The most damaging frequencies to power devices and machines appear to be the lower – below 5-kHz – frequency range. In years past, the magnitudes and sources of these lowerfrequency harmonics were limited and, inmost cases, power systems could tolerate them. The increase in power loss due to harmonics was also neglected because energy costs were low. These conditions no longer apply, and concern for harmonics is now becoming

For more Than 100 years, harmonics have been reported to cause operational problems to

**2. Sources and problems of harmonics [2]**  Harmonic sources are divided into two categories:

1. Established and known 2. New and Future

5. Transformer magnetizing currents.

are rich in harmonics.

sharp rise in harmonics.

widespread among utilities.

the power systems. Some of the major effects include:

voltage controllers, frequency converters, etc.

2. Motor control devices such as speed controls for traction. 3. High-voltage direct-current power conversion and transmission.

amount of power rectification for battery charging.

8. Cyclo-converters used for low-speed high-torque machines. 9. Pulse-burst-modulated heating elements for large furnaces.

network is also subjected to new harmonic sources:

continuously variable-var sources.


Equation (1) can be written as

$$w(t) = \sum\_{n=0}^{N} \left( V\_n \cos \phi\_n \sin \alpha\_0 t + V\_n \sin \phi\_n \cos \alpha\_0 t \right) \tag{2}$$

Define

$$\mathbf{x}\_n = V\_n \cos \phi\_n \tag{3a}$$

$$y\_u = V\_n \sin \phi\_u \tag{3b}$$

Electric Power Systems Harmonics - Identification and Measurements 7

The above estimation procedures are simple and straight forward if the voltage and/or current waveforms under investigation are stationary and in steady state, but if there is a sudden variation in the power system operation, transient operation, such as fault, lighting and sudden loading to the system or sudden switching off a large load, the voltage signal waveform may contain, for a few cycles, a dc component, which if it is neglected, will affect the harmonics estimation content in the waveform. To overcome this problem, the voltage

> 0 0 1

The exponential term in equation (10) can be expanded using Taylor series and its first two

 <sup>0</sup> 0 0 1

*n <sup>V</sup> vt V t V n t*

0

 11 12 <sup>0</sup> <sup>0</sup> 1

If the voltage *v*(*t*) is samples at a pre-selected rate *t*, then *m* sample would be obtained at *t*1,

*Zt BtY t*

*B*(*t*) is *m* (2*N* +2) measurement matrix whose elements depend on the initial and sampling times and its order depends on the number of harmonics and the number of

If *m* > (2*N* +2), we obtain over determined set of equation and the non-recursive least error

<sup>1</sup> \* *T T Y B tBt B tZt*

*Y* is (2*N* +2) 1 parameters vector to be estimated containing *x11,* x12 and *xn, yn,*

sin cos

*n n*

12 *V*

*x*

*N*

*n v t x tx x n t*

*t*2 = *t*1 + *t*, …, *tm* = 1 + (*m* – 1)*t*, in this case equation (13) becomes

terms chosen from Taylor series expansion for the exponential term.

square algorithm can be used to solve this system of equation as

*N*

sin

*n n*

 (11)

 *y n t*

(15)

(13)

11 0 *x V* (12a)

(12b)

(14)

*n n v t Ve V n t* 

sin

 (10)

signal in equation (1) may be remodeled to take into account the dc component as [4]

*<sup>N</sup> <sup>t</sup>*

*V*0 is the amplitude of decaying dc component at *t* = 0

Is the time constant of the decaying dc component

where

where

terms can be used as

Define the new parameters

Then, equation (11) can be written as

*Z*(*t*) is the *m* 1 voltage samples

(*t*) is *m* 1 error vector to be minimized.

Then, equation (2) can be written as

$$w(t) = \sum\_{n=0}^{N} \left(\mathbf{x}\_n \sin n\alpha\_0 t + y\_n \cos n\alpha\_0 t\right) \tag{4}$$

If the voltage signal *v*(*t*) is sampled at a pre-selected rate, say *t*, then *m* samples would be obtained at *t*1, *t*2 = *t*1 + *t*, *t*3 = *t*1 + 2*t*, …, *tm* = *t*1 + (*m* – 1)*t*. Then, after expanding equation (4), it can be written as

$$
\begin{bmatrix} v(t\_1) \\ v(t\_2) \\ \vdots \\ v(t\_m) \end{bmatrix} = \begin{bmatrix} a\_{11}(t\_1) & a\_{12}(t\_1) & \dots & a\_{12N+1}(t\_1) \\ a\_{21}(t\_2) & a\_{22}(t\_2) & \dots & a\_{22N+1}(t\_2) \\ \vdots \\ a\_{m1}(t\_m) & a\_{m2}(t\_m) & \dots & a\_{m2N+1}(t\_m) \end{bmatrix} \begin{bmatrix} y\_0 \\ \mathbf{x}\_1 \\ \mathbf{y}\_1 \\ \vdots \\ y\_N \end{bmatrix} \tag{5}
$$

where the elements of the *A* matrix are the sine and cosine expansion of equation (4). In the a's vector form, equation (5) can be written as

$$
\underline{Z}\_v(t) = A(t)\underline{\theta}\_v + \underline{\underline{\pi}}\_v(t) \tag{6}
$$

where *Zv*(*t*) is *m* 1 vector of sampled voltage measurement, *A*(*t*) is *m* (2*N* + 1) matrix of measurement coefficients, *<sup>v</sup>* is (2*N* +1) vector to be estimated, *v*(*t*) is *m* 1 error vector to be minimized. The order of the matrix *A*(*t*) depends n the number of harmonics to be estimated. Furthermore, the elements of the matrix *A*(*t*) depend on the initial sampling time *t*1 the sampling rate *t* and the data window size used in the estimation process. The matrix *A* (*t*) can be calculated on off-line and stored.

At least (2*N* + 1) samples are required to solve the problem formulated in (6). Using 2*N* + 1 samples may produce a poor estimate, since we force *v*(*t*) to be zero. We assume that *m* > 2*N* +1, so that equation (4) represents over determined set of equations.

#### **3.1.1 Time domain estimation; least error squares estimation (LES)**

The solution to the over determined set of equations of (6) in the LES sense is given by

$$\begin{aligned} \underline{\theta}\_v^\cdot &= \left[ A^\top(t) A(t) \right]^{-1} A^\top(t) \underline{Z}\_v(t) \\ &= A^\cdot(t) \underline{Z}\_v(t) \end{aligned} \tag{7}$$

where <sup>1</sup> ( ) [ ( ) ( )] ( ) *T T A t A tAt A t* is the left pseudo inverse. Having obtained the \* *v* , the magnitude of any harmonic of order *n* can be calculated as

$$\mathbf{V}\_n = \left[\mathbf{x}\_n^2 + n\_n^2\right]^{\frac{1}{2}}; n = 1, \dots, N \tag{8}$$

while the phase angle of the *n*th harmonic is:

$$\phi\_n = \tan^{-1} \frac{\mathcal{Y}\_n}{\mathcal{X}\_n} \tag{9}$$

The above estimation procedures are simple and straight forward if the voltage and/or current waveforms under investigation are stationary and in steady state, but if there is a sudden variation in the power system operation, transient operation, such as fault, lighting and sudden loading to the system or sudden switching off a large load, the voltage signal waveform may contain, for a few cycles, a dc component, which if it is neglected, will affect the harmonics estimation content in the waveform. To overcome this problem, the voltage signal in equation (1) may be remodeled to take into account the dc component as [4]

$$w\left(t\right) = V\_0 e^{-t^{f\_T}} + \sum\_{n=1}^{N} V\_n \sin\left(n a \phi\_0 t + \phi\right) \tag{10}$$

where

6 Power Quality Harmonics Analysis and Real Measurements Data

*y V nn n* sin

0 0

If the voltage signal *v*(*t*) is sampled at a pre-selected rate, say *t*, then *m* samples would be obtained at *t*1, *t*2 = *t*1 + *t*, *t*3 = *t*1 + 2*t*, …, *tm* = *t*1 + (*m* – 1)*t*. Then, after expanding equation

> 

*<sup>y</sup> at at a t*

 

11 1 12 1 12 1 1 1 1 21 2 22 2 22 1 2 2 1

*v t x at at a t v t y*

1 2 2 1

where the elements of the *A* matrix are the sine and cosine expansion of equation (4). In the

*v v Z t At t <sup>v</sup>* 

where *Zv*(*t*) is *m* 1 vector of sampled voltage measurement, *A*(*t*) is *m* (2*N* + 1) matrix of

be minimized. The order of the matrix *A*(*t*) depends n the number of harmonics to be estimated. Furthermore, the elements of the matrix *A*(*t*) depend on the initial sampling time *t*1 the sampling rate *t* and the data window size used in the estimation process. The matrix

At least (2*N* + 1) samples are required to solve the problem formulated in (6). Using 2*N* + 1 samples may produce a poor estimate, since we force *v*(*t*) to be zero. We assume that *m* >

*A tAt A tZ t*

1

<sup>1</sup> tan *<sup>n</sup>*

*n y x*

The solution to the over determined set of equations of (6) in the LES sense is given by

where <sup>1</sup> ( ) [ ( ) ( )] ( ) *T T A t A tAt A t* is the left pseudo inverse. Having obtained the \*

*A tZ t*

*n*

 

<sup>1</sup> \* *T T v v v*

*at at a t*

*mm mm mN m*

*vt x n t y n t* 

*n n*

sin cos

*N N*

 

*<sup>v</sup>* is (2*N* +1) vector to be estimated, *v*(*t*) is *m* 1 error vector to

2 2 <sup>2</sup> *V xn n nn* ; *n N* 1, , (8)

(9)

(4)

0

 

2*N* +1, so that equation (4) represents over determined set of equations.

**3.1.1 Time domain estimation; least error squares estimation (LES)** 

magnitude of any harmonic of order *n* can be calculated as

while the phase angle of the *n*th harmonic is:

*n*

*N*

Then, equation (2) can be written as

 

*m*

*v t*

a's vector form, equation (5) can be written as

*A* (*t*) can be calculated on off-line and stored.

(4), it can be written as

measurement coefficients,

(3b)

0

(5)

*N*

(6)

(7)

*v* , the

*y*

*V*0 is the amplitude of decaying dc component at *t* = 0

Is the time constant of the decaying dc component

The exponential term in equation (10) can be expanded using Taylor series and its first two terms can be used as

$$\text{cov}\left(t\right) = V\_0 - \left(\frac{V\_0}{\tau}\right)t + \sum\_{n=1}^{N} V\_n \sin\left(n\phi\_0 t + \phi\_n\right) \tag{11}$$

Define the new parameters

$$\mathbf{x}\_{11} = V\_0 \tag{12a}$$

$$\mathcal{X}\_{12} = \frac{V\_0}{\tau} \tag{12b}$$

Then, equation (11) can be written as

$$w(t) = \mathbf{x}\_{11} - t\mathbf{x}\_{12} + \sum\_{n=1}^{N} (\mathbf{x}\_n \sin n\phi\_0 t + y\_n \cos n\phi\_0 t) \tag{13}$$

If the voltage *v*(*t*) is samples at a pre-selected rate *t*, then *m* sample would be obtained at *t*1, *t*2 = *t*1 + *t*, …, *tm* = 1 + (*m* – 1)*t*, in this case equation (13) becomes

$$
\underline{Z}(t) = B(t)\underline{Y} + \underline{\xi}(t) \tag{14}
$$

where

*Z*(*t*) is the *m* 1 voltage samples

*B*(*t*) is *m* (2*N* +2) measurement matrix whose elements depend on the initial and sampling times and its order depends on the number of harmonics and the number of terms chosen from Taylor series expansion for the exponential term.

*Y* is (2*N* +2) 1 parameters vector to be estimated containing *x11,* x12 and *xn, yn,*

(*t*) is *m* 1 error vector to be minimized.

If *m* > (2*N* +2), we obtain over determined set of equation and the non-recursive least error square algorithm can be used to solve this system of equation as

$$\underline{Y}^\* = \left[\boldsymbol{B}^T\left(\boldsymbol{t}\right)\boldsymbol{B}\left(\boldsymbol{t}\right)\right]^{-1}\boldsymbol{B}^T\left(\boldsymbol{t}\right)\underline{Z}\left(\boldsymbol{t}\right)\tag{15}$$

Electric Power Systems Harmonics - Identification and Measurements 9

The first bracket in Equation (19) presents the possible low or high frequency sinusoidal with a combination of exponential terms, while the second bracket presents the harmonics, whose frequencies, *wk*, *k* = 1, …, *M*, are greater than 50/60 c/s, that contaminated the voltage or current waveforms. If these harmonics are identified to a certain degree of accuracy, i.e. a large number of harmonics are chosen, and then the first bracket presents the error in the voltage or current waveforms. Now, assume that these harmonics are identified,

<sup>1</sup>

*i e t Ae wt Ae wt*

It is clear that this expression represents the general possible low or high frequency dynamic oscillations. This model represents the dynamic oscillations in the system in cases such as, the currents of an induction motor when controlled by variable speed drive. As a special case, if the sampling constants are equal to zero then the considered wave is just a summation of low frequency components. Without loss of generality and for simplicity, it

2 cos cos *<sup>N</sup> <sup>t</sup> it*

 

*i ii*

(20)

, *i* = 3, …, *N*.

1 1

Note that *wi wk*, but <sup>1</sup> *<sup>w</sup>*

then the error *e*(*t*) can be written as

Fig. 1. Actual recorded phase currents.

*wi <sup>i</sup>* 

Having obtained the parameters vector *Y*\*, the harmonics magnitude and phase angle can be obtained as

$$V\_n = \left[\underline{\boldsymbol{x}}\_n^2 + \underline{\boldsymbol{y}}\_n^2\right]^{\frac{1}{2}}\tag{16}$$

$$\phi\_n = \tan^{-1} \frac{\mathcal{Y}\_n}{\mathcal{X}\_n} \tag{17}$$

while the parameters of the dc component can be calculated as

$$V\_0 = \mathbf{x}\_{11} \tag{18a}$$

$$
\pi = \frac{\mathcal{X}\_{11}}{\mathcal{X}\_{12}} \tag{18b}
$$

Figure 1 gives actual recorded data for a three-phase dynamic load. The load is a variable frequency drive controlling a 3000 HP induction motor connected to an oil pipeline compressor [5]. Examining this curve reveals the following: (a) the waveform of the phase currents are not periodical; (b) there are low-frequency transients, which have frequencies not an integer number of the fundamental, we call them sub-harmonics, contaminating these waveforms, especially in the tips of the wave; and (c) the phase currents are not symmetrical. It can be concluded from these remarks that this waveform is contaminated with harmonics, as well as low frequency transients, this is due to the power electronic devices associated with the load.

#### **3.2 Modeling of sub-harmonics in time domain**

The sub-harmonics is a noise contaminated with a signal and having frequency which is not a multiple from the fundamental frequency (50/60 Hz), as given in equation (19). To measure these sub-harmonics, an accurate model is needed to present the voltage and current waves:

Assume the voltage or current waveform is contaminated with both harmonics and subharmonics. Then, the waveform can be written as

$$f(t) = \left[A\_1 e^{\sigma 1t} \cos w\_1 t + \sum\_{l=2}^{N} A\_l e^{\sigma il} \cos \left(w\_l t + \phi\_l \right) \right] + \left[\sum\_{k=1}^{M} B\_k \cos \left(w\_k t + \theta\_k \right) \right] \tag{19}$$

where

*A*1, *A*2, …, *AN* are the sub-harmonics magnitude

*B*1, *B*2, …, *Bk* are the harmonics magnitude


$$w\_{\hat{\nu}} \text{ i = 1, \dots, N} \newline \text{ are the sub-harmonic frequencies, assumed to be identified in the } \hat{\nu} \text{ frequency domain}$$

*wk*; *k* = 1, …, *M* are the harmonic frequencies assumed to be identified also in the frequency domain.

Note that *wi wk*, but <sup>1</sup> *<sup>w</sup> wi <sup>i</sup>* , *i* = 3, …, *N*.

8 Power Quality Harmonics Analysis and Real Measurements Data

Having obtained the parameters vector *Y*\*, the harmonics magnitude and phase angle can be

<sup>1</sup> tan *<sup>n</sup>*

11 12 *x x* 

Figure 1 gives actual recorded data for a three-phase dynamic load. The load is a variable frequency drive controlling a 3000 HP induction motor connected to an oil pipeline compressor [5]. Examining this curve reveals the following: (a) the waveform of the phase currents are not periodical; (b) there are low-frequency transients, which have frequencies not an integer number of the fundamental, we call them sub-harmonics, contaminating these waveforms, especially in the tips of the wave; and (c) the phase currents are not symmetrical. It can be concluded from these remarks that this waveform is contaminated with harmonics, as well as low frequency transients, this is due to the power electronic

The sub-harmonics is a noise contaminated with a signal and having frequency which is not a multiple from the fundamental frequency (50/60 Hz), as given in equation (19). To measure these sub-harmonics, an accurate model is needed to present the voltage and

Assume the voltage or current waveform is contaminated with both harmonics and sub-

( ) cos cos cos *N M t it*

*wi*; i=1, …, *N* are the sub-harmonic frequencies, assumed to be identified in the

*wk*; *k* = 1, …, *M* are the harmonic frequencies assumed to be identified also in the

 

*f t Ae wt Ae wt B wt*

<sup>1</sup>

(19)

*i ii k kk*

 

2 1

*i k*

*n y x*

*n*

while the parameters of the dc component can be calculated as

devices associated with the load.

current waves:

where

1, 2, …, 

frequency domain

frequency domain.

**3.2 Modeling of sub-harmonics in time domain** 

harmonics. Then, the waveform can be written as

*A*1, *A*2, …, *AN* are the sub-harmonics magnitude *B*1, *B*2, …, *Bk* are the harmonics magnitude

*<sup>N</sup>* are the damping constants

*<sup>i</sup>*; *i* = 1, …, *N* are the sub-harmonic phase angles

*<sup>k</sup>*; k = 1, …, *M* are the harmonic phase angles

1 1

1

2 2 <sup>2</sup> *V xy n nn* (16)

(17)

*V x* 0 11 (18a)

(18b)

obtained as

The first bracket in Equation (19) presents the possible low or high frequency sinusoidal with a combination of exponential terms, while the second bracket presents the harmonics, whose frequencies, *wk*, *k* = 1, …, *M*, are greater than 50/60 c/s, that contaminated the voltage or current waveforms. If these harmonics are identified to a certain degree of accuracy, i.e. a large number of harmonics are chosen, and then the first bracket presents the error in the voltage or current waveforms. Now, assume that these harmonics are identified, then the error *e*(*t*) can be written as

$$e(t) = A\_1 e^{\sigma \text{1} t} \cos w\_1 t + \sum\_{i=2}^{N} A\_i e^{\sigma \text{1} t} \cos \left(w\_i t + \phi\_i \right) \tag{20}$$

Fig. 1. Actual recorded phase currents.

It is clear that this expression represents the general possible low or high frequency dynamic oscillations. This model represents the dynamic oscillations in the system in cases such as, the currents of an induction motor when controlled by variable speed drive. As a special case, if the sampling constants are equal to zero then the considered wave is just a summation of low frequency components. Without loss of generality and for simplicity, it

Electric Power Systems Harmonics - Identification and Measurements 11

It is clear that this set of equations is similar to the set of equations given by equation (5).

*zt Ht t t* 

where *z*(*t*) is the vector of sampled measurements, *H(t*) is an *m* 6, in this simple case,

is an *m* 1 noise vector to be minimized. The dimensions of the previous matrices depend on the number of modes considered, as well as, the number of terms truncated from the

<sup>1</sup> \* *T T*

\* 2 11 1 \*

\*

1

3

*x*

1 \* \*2 \*2 2 4 2 35 2 \*

\* \* 5 6 2 \* \* 3 4

*x x x x*

*t H tHt H tZt*

, *<sup>x</sup> A x <sup>x</sup>* 

, *<sup>x</sup> Axx*

tan

*m*

advance, the left pseudo inverse of [*A*] can be determined for an application off-line.

uses a data window of m samples to provide an estimate of the unknowns,

In the least error squares estimates explained in the previous section, the estimated

 <sup>1</sup> \* 1 1

where [*A*]+ is the left pseudo inverse of [*A*] = [*ATA*]-1*AT*, the superscript "*m* – 1" in the equation represents the estimates calculated using data taken from *t* = *t*1 to *t* = *t*1 + (*m* – 1)*t* s, *t*1 is the initial sampling time. The elements of the matrix [*A*] are functions of the time reference, initial sampling time, and the sampling rate used. Since these are selected in

Equation 33 represents, as we said earlier, a non recursive least error squares (LES) filter that

] are calculated by taking the row products of the matrix [*A*]+ with the *m* samples. A new sample is included in the data window at each sampling interval and the oldest sample is discarded. The new [*A*]+ for the latest m samples is calculated and the estimates of [

 *A Z* 

*n nm m*

 

(28)

(*t*)

(*t*) is a 6 1 parameter vector to be estimated, and

(29)

(31)

(33)

(32)

\*(*t*), then the sub-harmonics parameters can be

(30)

. The estimates

] are

Thus an equation similar to (6) can be written as:

matrix of measurement coefficients,

**3.2.1 Least error squares estimation** 

Having obtained the parameters vector

**3.2.2 Recursive least error squares estimates** 

parameters, in the three cases, take the form of

The solution to equation (28) based on LES is given as

Taylor series.

obtained as

of [

can be assumed that only two modes of equation (21) are considered, then the error *e*(*t*) can be written as (21)

$$e(t) = A\_1 e^{\sigma 1t} \cos \left( w\_1 t \right) + A\_2 e^{\sigma 2t} \cos \left( w\_2 t + \phi\_2 \right) \tag{21}$$

Using the well-known trigonometric identity

$$\cos\left(w\_2 t + \phi\_2\right) = \cos w\_2 t \cos \phi\_2 - \sin w\_2 t \sin \phi\_2$$

then equation (21) can be rewritten as:

$$e(t) = A\_1 e^{\sigma^{1t}t} \cos w\_1 t + \left(e^{\sigma^{2t}t} \cos w\_2 t\right) A\_2 \cos \phi\_2 + \left(e^{\sigma^{2t}t} \sin w\_2 t\right) A\_2 \sin \phi\_2 \tag{22}$$

It is obvious that equation (22) is a nonlinear function of *A'*s, *'*s and *'*s. By using the first two terms in the Taylor series expansion *Aie it*; *i* = 1,2. Equation (22) turns out to be

$$\begin{aligned} e(t) &= A\_1 \cos w\_1 t + \left(t \cos w\_1 t\right) \left(A\_1 \sigma\_1\right) + \left(\cos w\_2 t\right) \left(A\_2 \cos \phi\_2\right) + \left(t \cos w\_2 t\right) \left(A\_2 \sigma\_2 \cos \phi\_2\right) \\ &- \left(\sin w\_2 t\right) \left(A\_2 \sin \phi\_2\right) + \left(t \sin w\_2 t\right) \left(A\_2 \sigma\_2 \sin \phi\_2\right) \end{aligned} \tag{23}$$

where the Taylor series expansion is given by:

1 *<sup>t</sup> e t* 

Making the following substitutions in equation (23), equation (26) can be obtained,

$$\begin{cases} \mathbf{x}\_1 = A\_1; & \mathbf{x}\_2 = A\_1 \sigma\_1 \\ \mathbf{x}\_3 = A\_2 \cos \phi\_2; & \mathbf{x}\_4 = A\_2 \sigma\_2 \cos \phi\_2 \\ \mathbf{x}\_5 = A\_2 \sin \phi\_2; & \mathbf{x}\_6 = A\_2 \sigma\_2 \sin \phi\_2 \end{cases} \tag{24}$$

and

$$\begin{cases} h\_{11}\left(t\right) = \cos w\_1 t; & h\_{12}\left(t\right) = t \cos w\_1 t \\ h\_{13}\left(t\right) = \cos w\_2 t; & h\_{14}\left(t\right) = t \cos w\_2 t \\ h\_{15}\left(t\right) = -\sin w\_2 t; & h\_{16}\left(t\right) = -t \sin w\_2 t \end{cases} \tag{25}$$

$$e(\mathbf{t}) = h\_{11}(\mathbf{t})\mathbf{x}\_1 + h\_{12}(\mathbf{t})\mathbf{x}\_2 + h\_{13}(\mathbf{t})\mathbf{x}\_3 + h\_{14}(\mathbf{t})\mathbf{x}\_3 + h\_{15}(\mathbf{t})\mathbf{x}\_4 + h\_{16}(\mathbf{t})\mathbf{x}\_5 \tag{26}$$

If the function *f* (*t*) is sampled at a pre-selected rate, its samples would be obtained at equal time intervals, say *t* seconds. Considering m samples, then there will be a set of *m* equations with an arbitrary time reference *t*1 given by

$$\begin{vmatrix} \mathbf{e}\begin{pmatrix} t\_1\\ e\begin{pmatrix} t\_2\\ e\end{pmatrix} \end{pmatrix} = \begin{vmatrix} h\_{11}\begin{pmatrix} t\_1\\ t\_2 \end{pmatrix} & h\_{12}\begin{pmatrix} t\_1\\ t\_2 \end{pmatrix} & \dots & h\_{16}\begin{pmatrix} t\_1\\ t\_2 \end{pmatrix} & \dots & h\_{26}\begin{pmatrix} t\_2\\ t\_3 \end{pmatrix} \\\ \dots & \dots & \dots & \dots & \dots\\\ \dots & \dots & \dots & \dots & \dots\\\ \mathbf{h}\_{m1}\begin{pmatrix} t\_m \end{pmatrix} & h\_{m2}\begin{pmatrix} t\_m\\ t\_m \end{pmatrix} & \dots & h\_{ms}\begin{pmatrix} t\_m \end{pmatrix} \end{pmatrix} \tag{27}$$

It is clear that this set of equations is similar to the set of equations given by equation (5). Thus an equation similar to (6) can be written as:

$$
\underline{z}(t) = H(t)\underline{\theta}(t) + \underline{\xi}(t) \tag{28}
$$

where *z*(*t*) is the vector of sampled measurements, *H(t*) is an *m* 6, in this simple case, matrix of measurement coefficients, (*t*) is a 6 1 parameter vector to be estimated, and (*t*) is an *m* 1 noise vector to be minimized. The dimensions of the previous matrices depend on the number of modes considered, as well as, the number of terms truncated from the Taylor series.

#### **3.2.1 Least error squares estimation**

10 Power Quality Harmonics Analysis and Real Measurements Data

can be assumed that only two modes of equation (21) are considered, then the error *e*(*t*) can

 1 2 1 1 2 22 cos cos *t t e t Ae wt Ae wt*

cos*wt wt wt* 22 2 2 2 2

1 2 <sup>2</sup>

*e t A wt t wt A wt A t wt A*

**2**

Making the following substitutions in equation (23), equation (26) can be obtained,

;

*t t <sup>t</sup> e t Ae wt e wt A e wt A*

 

1 *<sup>t</sup> e t* 

1 1 2 11 3 2 2 4 22 2 5 2 2 6 22 2

11 1 12 1 13 2 14 2 15 2 16 2

If the function *f* (*t*) is sampled at a pre-selected rate, its samples would be obtained at equal time intervals, say *t* seconds. Considering m samples, then there will be a set of *m*

> 

1 11 1 12 1 16 1 1 2 21 2 22 2 26 2

*et h t h t h t x et h t h t h t x*

 

 

12 6

*m mm mm m m*

*et h t h t h t*

 

*h t wt h t t wt h t wt h t t wt h t wt h t t wt* 

cos ; cos cos ; cos sin ; sin

*xA xA xA xA xA xA*

 

cos ; cos sin ; sin

1 1 1 11 22 2 22 2

 

 

*et h tx h tx h tx h tx h tx h tx* 11 1 12 2 13 3 14 3 15 4 16 5 (26)

 

2

 

6 *x*

cos cos cos cos cos cos

 cos cos sin sin 

1 1 22 2 22 2 cos cos cos sin sin

   

*'*s and 

*it*; *i* = 1,2. Equation (22) turns out to be

(21)

 

(24)

(25)

(27)

**2**

(22)

*'*s. By using the first

(23)

 

It is obvious that equation (22) is a nonlinear function of *A'*s,

sin sin sin sin

*wt A t wt A*

equations with an arbitrary time reference *t*1 given by

 

22 2 22 2

Using the well-known trigonometric identity

two terms in the Taylor series expansion *Aie*

where the Taylor series expansion is given by:

then equation (21) can be rewritten as:

be written as (21)

and

The solution to equation (28) based on LES is given as

$$\underline{\boldsymbol{\theta}}^{\ast}(t) = \left[\boldsymbol{H}^{\top}(t)\boldsymbol{H}(t)\right]^{-1}\boldsymbol{H}^{\top}(t)\underline{\boldsymbol{Z}}(t) \tag{29}$$

Having obtained the parameters vector \*(*t*), then the sub-harmonics parameters can be obtained as

$$A\_1 = \stackrel{\*}{\mathfrak{x}\_{1'}^\*} \qquad \sigma\_1 = \frac{\stackrel{\*}{\mathfrak{x}\_2^\*}}{\stackrel{\*}{\mathfrak{x}\_1^\*}} \tag{30}$$

$$A\_2 = \left[\boldsymbol{\chi^\*\_{\mathfrak{z}}}^\* + \boldsymbol{\chi^\*\_{\mathfrak{z}}}^\*\right]^{\frac{1}{2}}, \qquad \sigma\_2 = \frac{\boldsymbol{\chi^\*\_{\mathfrak{z}}}}{\boldsymbol{\chi^\*\_{\mathfrak{z}}}} \tag{31}$$

$$\tan \phi\_2 = \frac{\mathbf{x}\_5^\ast}{\mathbf{x}\_3^\ast} = \frac{\mathbf{x}\_6^\ast}{\mathbf{x}\_4^\ast} \tag{32}$$

#### **3.2.2 Recursive least error squares estimates**

In the least error squares estimates explained in the previous section, the estimated parameters, in the three cases, take the form of

$$\left[\underline{\underline{\boldsymbol{\sigma}}}^{\*}\right]\_{n\times 1}^{m-1} = \left[\boldsymbol{A}\right]\_{n\times m}^{\*} \left[\underline{\underline{\boldsymbol{Z}}}\right]\_{m\times 1} \tag{33}$$

where [*A*]+ is the left pseudo inverse of [*A*] = [*ATA*]-1*AT*, the superscript "*m* – 1" in the equation represents the estimates calculated using data taken from *t* = *t*1 to *t* = *t*1 + (*m* – 1)*t* s, *t*1 is the initial sampling time. The elements of the matrix [*A*] are functions of the time reference, initial sampling time, and the sampling rate used. Since these are selected in advance, the left pseudo inverse of [*A*] can be determined for an application off-line.

Equation 33 represents, as we said earlier, a non recursive least error squares (LES) filter that uses a data window of m samples to provide an estimate of the unknowns, . The estimates of [] are calculated by taking the row products of the matrix [*A*]+ with the *m* samples. A new sample is included in the data window at each sampling interval and the oldest sample is discarded. The new [*A*]+ for the latest m samples is calculated and the estimates of [] are

Electric Power Systems Harmonics - Identification and Measurements 13

 \* *r Zt At At Zt* 

1

*m n*

as

**Step 8.** Calculate the LAV residual generated from this solution

*r r*

1

 

*m i i*

1

**Step 4.** Reject the measurements having residuals greater than the standard deviation, and

 \* <sup>1</sup> <sup>1</sup> <sup>1</sup> ˆ ˆ <sup>ˆ</sup> *<sup>n</sup> <sup>n</sup> n n*

 *<sup>A</sup> t Zt* 

Ref. 6 carried out a comparative study for power system harmonic estimation. Three algorithms are used in this study; LES, LAV, and discrete Fourier transform (DFT). The data used in this study are real data from a three-phase six pulse converter. The three techniques are thoroughly analyzed and compared in terms of standard deviation, number of samples

For the purpose of this study, the voltage signal is considered to contain up to the 13th harmonics. Higher order harmonics are neglected. The rms voltage components are given in

RMS voltage components corresponding to the harmonics

Figure 2 shows the A.C. voltage waveform at the converter terminal. The degree of the distortion depends on the order of the harmonics considered as well as the system

The variables to be estimated are the magnitudes of each voltage harmonic from the fundamental to the 13th harmonic. The estimation is performed by the three techniques while several parameters are changed and varied. These parameters are the standard

characteristics. Figure 3 shows the spectrum of the converter bus bar voltage.

0.95–2.02 0.0982. 0.0438.9 0.030212.9 0.033162.6

frequency Fundamental 5th 7th 11th 13th

**Step 5.** Recalculated the least error squares residuals generated from this new solution **Step 6.** Rank the residual and select *n* measurements corresponding to the smallest

1 2 2

**Step 3.** Calculated the standard deviation of this residual vector as

Where

1 1 *<sup>m</sup> i i r r m*

residuals

, the average residual

recalculate the LES solution

**Step 7.** Solve for the LAV estimates ˆ

**3.3 Computer simulated tests** 

and sampling frequency.

Table 1.

Table 1.

Harmonic

Voltage magnitude (p.u.)

updated by taking the row products of the updated [*A*]+ with the latest m samples. However, equation (33) can be modified to a recursive form which is computationally more efficient.

Recall that equation

$$\begin{bmatrix} \underline{Z} \end{bmatrix}\_{m \times 1} = \begin{bmatrix} A \end{bmatrix}\_{m \times n} \begin{bmatrix} \underline{\theta} \end{bmatrix}\_{n \times 1} \tag{34}$$

represents a set of equations in which [*Z*] is a vector of m current samples taken at intervals of *t* seconds. The elements of the matrix [*A*] are known. At time *t* = *t*1 + *mt* a new sample is taken. Then equation (33) can be written as

$$
\begin{bmatrix} \boldsymbol{\theta} \ \end{bmatrix}\_{n \times 1}^{m} = \begin{bmatrix} \boldsymbol{A} \ \boldsymbol{A} \\ \boldsymbol{a}\_{m} \ \end{bmatrix}\_{n \times \{m\boldsymbol{H}\}}^{\*} \begin{bmatrix} \begin{bmatrix} \boldsymbol{Z} \\ \end{bmatrix} \end{bmatrix}\_{\{m\boldsymbol{H}\} \times 1} \tag{35}
$$

where the superscript "*m*" represents the new estimate at time *t* = *t*1 + *mt*. It is possible to express the new estimates obtained from equation (34) in terms of older estimates (obtained from equation (33)) and the latest sample *Zm* as follows

$$
\begin{bmatrix} \boldsymbol{\theta}^\* \end{bmatrix}^{\boldsymbol{m}} = \begin{bmatrix} \boldsymbol{\theta}^\* \end{bmatrix}^{\boldsymbol{m}-1} + \begin{bmatrix} \boldsymbol{a}(\boldsymbol{m}) \end{bmatrix} \begin{bmatrix} \boldsymbol{Z}\_{\boldsymbol{m}} \end{bmatrix} - \begin{bmatrix} \boldsymbol{a}\_{\boldsymbol{m}\boldsymbol{i}} \end{bmatrix} \begin{bmatrix} \boldsymbol{\theta}^\* \end{bmatrix}^{\boldsymbol{m}-1} \end{bmatrix} \tag{36a}
$$

This equation represents a recursive least squares filter. The estimates of the vector [] at *t* = *t*1 + *mt* are expressed as a function of the estimates at *t* = *t*1 + (*m* – 1)*t* and the term <sup>1</sup> \* *<sup>m</sup> Z a m mi* . The elements of the vector, [(*m*)], are the time-invariant gains

of the recursive least squares filter and are given as

$$\alpha(m) = \left[A^T A\right]^{-1} \left[a\_{mi}\right]^T \left[\left[I\right] + \left[a\_{mi}\right]\left[A^T A\right]^{-1} \left[a\_{mi}\right]^T\right]^{-1} \tag{36b}$$

#### **3.2.3 Least absolute value estimates (LAV) algorithm (Soliman & Christensen algorithm) [3]**

The LAV estimation algorithm can be used to estimate the parameters vectors. For the reader's convenience, we explain here the steps behind this algorithm.

Given the observation equation in the form of that given in (28) as

$$
\underline{Z}(t) = A(t)\underline{\theta} + \epsilon \in (t).
$$

The steps in this algorithm are:

**Step 1.** Calculate the LES solution given by

$$\left[\underline{\underline{\boldsymbol{\theta}}}^{\*}\right] = \left[\boldsymbol{A}(t)\right]^{\*} \underline{\boldsymbol{Z}}(t) \cdot \left[\boldsymbol{A}(t)\right]^{\*} = \left[\boldsymbol{A}^{\top}(t)\boldsymbol{A}(t)\right]^{-1} \boldsymbol{A}^{\top}(t)$$

**Step 2.** Calculate the LES residuals vector generated from this solution as

$$\underline{\underline{r}}^\* = \underline{Z}(t) - A(t) \underline{\underline{\underline{\underline{\underline{\underline{\underline{\underline{\underline{\underline{\underline{\underline{\cdot}}}}}}}}}}}{}{}} \underline{\underline{\underline{\underline{\cdot}}}}}(t)}$$

**Step 3.** Calculated the standard deviation of this residual vector as

$$\sigma = \frac{1}{m - n + 1} \left[ \sum\_{i=1}^{m} \left( r\_i - \overline{r} \right)^2 \right]^{\frac{1}{2}}$$

Where 1 1 *<sup>m</sup> i i r r m* , the average residual

12 Power Quality Harmonics Analysis and Real Measurements Data

updated by taking the row products of the updated [*A*]+ with the latest m samples. However, equation (33) can be modified to a recursive form which is computationally more

*m mn n* 1 1 *Z A*

represents a set of equations in which [*Z*] is a vector of m current samples taken at intervals of *t* seconds. The elements of the matrix [*A*] are known. At time *t* = *t*1 + *mt* a new sample

*<sup>n</sup> mi n mH <sup>m</sup> mH A Z a Z*

where the superscript "*m*" represents the new estimate at time *t* = *t*1 + *mt*. It is possible to express the new estimates obtained from equation (34) in terms of older estimates (obtained

 1 1 \* \* \* *m m <sup>m</sup>*

This equation represents a recursive least squares filter. The estimates of the vector [

. The elements of the vector, [

**3.2.3 Least absolute value estimates (LAV) algorithm (Soliman & Christensen** 

*t*1 + *mt* are expressed as a function of the estimates at *t* = *t*1 + (*m* – 1)*t* and the term

<sup>1</sup> 1 1 *T T T T m AA a I a AA a mi mi mi*

The LAV estimation algorithm can be used to estimate the parameters vectors. For the

*Zt At t* 

*At Zt* , <sup>1</sup> *T T At A tAt A t*

   

*mZ a m mi*

1

(34)

(35)

(*m*)], are the time-invariant gains

(36a)

] at *t* =

(36b)

efficient.

Recall that equation

is taken. Then equation (33) can be written as

\* 1

from equation (33)) and the latest sample *Zm* as follows

of the recursive least squares filter and are given as

reader's convenience, we explain here the steps behind this algorithm. Given the observation equation in the form of that given in (28) as

\*

**Step 2.** Calculate the LES residuals vector generated from this solution as

<sup>1</sup> \* *<sup>m</sup>*

The steps in this algorithm are:

**Step 1.** Calculate the LES solution given by

*Z a m mi*

**algorithm) [3]** 

*m*


$$
\left[\left[\hat{\theta}\right]\right]\_{n\times 1} = \left[\hat{A}\left(t\right)\right]\_{n\times n}^{-1} \hat{Z}\_{n\times 1}\left(t\right).
$$

**Step 8.** Calculate the LAV residual generated from this solution

#### **3.3 Computer simulated tests**

Ref. 6 carried out a comparative study for power system harmonic estimation. Three algorithms are used in this study; LES, LAV, and discrete Fourier transform (DFT). The data used in this study are real data from a three-phase six pulse converter. The three techniques are thoroughly analyzed and compared in terms of standard deviation, number of samples and sampling frequency.

For the purpose of this study, the voltage signal is considered to contain up to the 13th harmonics. Higher order harmonics are neglected. The rms voltage components are given in Table 1.


Table 1.

Figure 2 shows the A.C. voltage waveform at the converter terminal. The degree of the distortion depends on the order of the harmonics considered as well as the system characteristics. Figure 3 shows the spectrum of the converter bus bar voltage.

The variables to be estimated are the magnitudes of each voltage harmonic from the fundamental to the 13th harmonic. The estimation is performed by the three techniques while several parameters are changed and varied. These parameters are the standard

Electric Power Systems Harmonics - Identification and Measurements 15

Figure 4 shows the effects of number of samples on the fundamental component magnitude using the three techniques at a sampling frequency = 1620 Hz and the measurement set is

corrupted with a noise having standard deviation of 0.1 Gaussian distribution.

Fig. 4. Effect of number of samples on the magnitude estimation of the fundamental

harmonic magnitudes. Examining these figures reveals the following remarks.

the fact that the number of samples per number of cycles is not an integer.

The LAV method gives better estimates for most of the number of samples.

At a low number of samples, the LES produces poor estimates.

It can be noticed from this figure that the DFT algorithm gives an essentially exact estimate of the fundamental voltage magnitude. The LAV algorithm requires a minimum number of samples to give a good estimate, while the LES gives reasonable estimates over a wide range of numbers of samples. However, the performance of the LAV and LES algorithms is improved when the sampling frequency is increased to 1800 Hz as shown in Figure 5. Figure 6 –9 gives the same estimates at the same conditions for the 5th, 7th, 11th and 13th

 For all harmonics components, the DFT gives bad estimates for the magnitudes. This bad estimate is attributed to the phenomenon known as "spectral leakage" and is due to

As the number of samples increases, the LES method gives a relatively good performance.

However, as the sampling frequency increased to 1800 Hz, no appreciable effects have changed, and the estimates of the harmonics magnitude are still the same for the three

harmonic (sampling frequency = 1620 Hz).

techniques.

deviation of the noise, the number of samples, and the sampling frequency. A Gaussiandistributed noise of zero mean was used.

Fig. 2. AC voltage waveform

Fig. 3. Frequency spectrums.

deviation of the noise, the number of samples, and the sampling frequency. A Gaussian-

.

distributed noise of zero mean was used.

Fig. 2. AC voltage waveform

Fig. 3. Frequency spectrums.

Figure 4 shows the effects of number of samples on the fundamental component magnitude using the three techniques at a sampling frequency = 1620 Hz and the measurement set is corrupted with a noise having standard deviation of 0.1 Gaussian distribution.

Fig. 4. Effect of number of samples on the magnitude estimation of the fundamental harmonic (sampling frequency = 1620 Hz).

It can be noticed from this figure that the DFT algorithm gives an essentially exact estimate of the fundamental voltage magnitude. The LAV algorithm requires a minimum number of samples to give a good estimate, while the LES gives reasonable estimates over a wide range of numbers of samples. However, the performance of the LAV and LES algorithms is improved when the sampling frequency is increased to 1800 Hz as shown in Figure 5.

Figure 6 –9 gives the same estimates at the same conditions for the 5th, 7th, 11th and 13th harmonic magnitudes. Examining these figures reveals the following remarks.


However, as the sampling frequency increased to 1800 Hz, no appreciable effects have changed, and the estimates of the harmonics magnitude are still the same for the three techniques.

Electric Power Systems Harmonics - Identification and Measurements 17

Fig. 7. Effect of number of samples on the magnitude estimation of the 7th harmonic

Fig. 8. Effect of number of samples on the magnitude estimation of the 11th harmonic

(sampling frequency = 1620 Hz).

(sampling frequency = 1620 Hz).

Fig. 5. Effect of number of samples on the magnitude estimation of the fundamental harmonic (sampling frequency = 1800 Hz).

Fig. 6. Effect of number of samples on the magnitude estimation of the 5th harmonic (sampling frequency = 1620 Hz).

Fig. 5. Effect of number of samples on the magnitude estimation of the fundamental

Fig. 6. Effect of number of samples on the magnitude estimation of the 5th harmonic

harmonic (sampling frequency = 1800 Hz).

(sampling frequency = 1620 Hz).

Fig. 7. Effect of number of samples on the magnitude estimation of the 7th harmonic (sampling frequency = 1620 Hz).

Fig. 8. Effect of number of samples on the magnitude estimation of the 11th harmonic (sampling frequency = 1620 Hz).

Electric Power Systems Harmonics - Identification and Measurements 19

The CPU time for the DFT and LES algorithms are essentially the same, and that of the LAV algorithm is larger. As the number of samples increases, the difference in CPU time between

Other interesting studies have been carried out on the performance of the three algorithms when 10% of the data is missed, taken uniformly at equal intervals starting from the first data point, for the noise free signal and 0.1 standard deviation added white noise Gaussian,

Figure 11 gives the estimates of the three algorithms at the two cases. Examining this figure

For the no noise estimates, the LS and DFT produce bad estimates for the fundamental

Fig. 11. Effect of number of samples on the magnitude estimation of the fundamental harmonic for 10% missing data (sampling frequency = 1620 Hz): (a) no noise; (b) 0.1

standard deviation added white Gaussian noise.

the LAV and LS/DFT algorithm increases.

and the sampling frequency used is 1620 Hz.

harmonic magnitude, even at a higher number of samples

The LAV algorithm produces good estimates, at large number of samples.

we can notice the following remarks:

Fig. 9. Effect of number of samples on the magnitude estimation of the 13th harmonic (sampling frequency = 1620 Hz).

The CPU time is computed for each of the three algorithms, at a sampling frequency of 1620 Hz. Figure 10 gives the variation of CPU.

Fig. 10. The CPU times of the LS, DFT, and LAV methods (sampling frequency = 1620 Hz).

Fig. 9. Effect of number of samples on the magnitude estimation of the 13th harmonic

The CPU time is computed for each of the three algorithms, at a sampling frequency of 1620

Fig. 10. The CPU times of the LS, DFT, and LAV methods (sampling frequency = 1620 Hz).

(sampling frequency = 1620 Hz).

Hz. Figure 10 gives the variation of CPU.

The CPU time for the DFT and LES algorithms are essentially the same, and that of the LAV algorithm is larger. As the number of samples increases, the difference in CPU time between the LAV and LS/DFT algorithm increases.

Other interesting studies have been carried out on the performance of the three algorithms when 10% of the data is missed, taken uniformly at equal intervals starting from the first data point, for the noise free signal and 0.1 standard deviation added white noise Gaussian, and the sampling frequency used is 1620 Hz.

Figure 11 gives the estimates of the three algorithms at the two cases. Examining this figure we can notice the following remarks:

For the no noise estimates, the LS and DFT produce bad estimates for the fundamental harmonic magnitude, even at a higher number of samples

The LAV algorithm produces good estimates, at large number of samples.

Fig. 11. Effect of number of samples on the magnitude estimation of the fundamental harmonic for 10% missing data (sampling frequency = 1620 Hz): (a) no noise; (b) 0.1 standard deviation added white Gaussian noise.

Electric Power Systems Harmonics - Identification and Measurements 21

Fig. 13. Effect of number of samples on the magnitude estimation of the 7th harmonic for 10% missing data (sampling frequency = 1620 Hz): (a) no noise; (b) 0.1 standard deviation

In the previous section static-state estimation algorithms are implemented for identifying and measuring power system harmonics. The techniques used in that section was the least error squares (LES), least absolute value (LAV) and the recursive least error squares algorithms. These techniques assume that harmonic magnitudes are constant during the data window size used in the estimation process. In real time, due to the switching on-off of power electronic equipments (devices) used in electric derives and power system transmission (AC/DC transmission), the situation is different, where the harmonic magnitudes are not stationary during the data window size. As such a dynamic state estimation technique is required to identifying (tracking) the harmonic magnitudes as well

In this section, we introduce the Kalman filtering algorithm as well as the dynamic least absolute value algorithm (DLAV) for identifying (tracking) the power systems harmonics

The Kalman filtering approach provides a mean for optimally estimating phasors and the

The state variable representation of a signal that includes n harmonics for a noise-free

Each frequency component requires two state variables. Thus the total number of state

*st A t i t*

*Ai*(*t*) is the amplitude of the phasor quantity representing the *i*th harmonic at time t

cos

*i i*

 (37)

1

*i*

*i* is the phase angle of the *i*th harmonic relative to a reference rotating at *i*

*n*

added white Gaussian noise.

**4. Estimation of harmonics; the dynamic case** 

as the phase angles of each harmonics component.

current or voltage signal *s*(*t*) may be represented by [7]

variable is 2n. These state variables are defined as follows

and sub-harmonics (inter-harmonics).

*n* is the harmonic order

where

ability to track-time-varying parameters.

Figure 12 –15 give the three algorithms estimates, for 10% missing data with no noise and with 0.1 standard deviation Gaussian white noise, when the sampling frequency is 1620 Hz for the harmonics magnitudes and the same discussions hold true.

#### **3.4 Remarks**

Three signal estimation algorithms were used to estimate the harmonic components of the AC voltage of a three-phase six-pulse AC-DC converter. The algorithms are the LS, LAV, and DFT. The simulation of the ideal noise-free case data revealed that all three methods give exact estimates of all the harmonics for a sufficiently high sampling rate. For the noisy case, the results are completely different. In general, the LS method worked well for a high number of samples. The DFT failed completely. The LAV gives better estimates for a large range of samples and is clearly superior for the case of missing data.

Fig. 12. Effect of number of samples on the magnitude estimation of the 5th harmonic for 10% missing data (sampling frequency = 1620 Hz): (a) no noise; (b) 0.1 standard deviation added white Gaussian noise.

Fig. 13. Effect of number of samples on the magnitude estimation of the 7th harmonic for 10% missing data (sampling frequency = 1620 Hz): (a) no noise; (b) 0.1 standard deviation added white Gaussian noise.

#### **4. Estimation of harmonics; the dynamic case**

In the previous section static-state estimation algorithms are implemented for identifying and measuring power system harmonics. The techniques used in that section was the least error squares (LES), least absolute value (LAV) and the recursive least error squares algorithms. These techniques assume that harmonic magnitudes are constant during the data window size used in the estimation process. In real time, due to the switching on-off of power electronic equipments (devices) used in electric derives and power system transmission (AC/DC transmission), the situation is different, where the harmonic magnitudes are not stationary during the data window size. As such a dynamic state estimation technique is required to identifying (tracking) the harmonic magnitudes as well as the phase angles of each harmonics component.

In this section, we introduce the Kalman filtering algorithm as well as the dynamic least absolute value algorithm (DLAV) for identifying (tracking) the power systems harmonics and sub-harmonics (inter-harmonics).

The Kalman filtering approach provides a mean for optimally estimating phasors and the ability to track-time-varying parameters.

The state variable representation of a signal that includes n harmonics for a noise-free current or voltage signal *s*(*t*) may be represented by [7]

$$s(t) = \sum\_{i=1}^{n} A\_i(t) \cos(\imath\_i \alpha t + \theta\_i) \tag{37}$$

where

20 Power Quality Harmonics Analysis and Real Measurements Data

Figure 12 –15 give the three algorithms estimates, for 10% missing data with no noise and with 0.1 standard deviation Gaussian white noise, when the sampling frequency is 1620 Hz

Three signal estimation algorithms were used to estimate the harmonic components of the AC voltage of a three-phase six-pulse AC-DC converter. The algorithms are the LS, LAV, and DFT. The simulation of the ideal noise-free case data revealed that all three methods give exact estimates of all the harmonics for a sufficiently high sampling rate. For the noisy case, the results are completely different. In general, the LS method worked well for a high number of samples. The DFT failed completely. The LAV gives better estimates for a large

Fig. 12. Effect of number of samples on the magnitude estimation of the 5th harmonic for 10% missing data (sampling frequency = 1620 Hz): (a) no noise; (b) 0.1 standard deviation

added white Gaussian noise.

for the harmonics magnitudes and the same discussions hold true.

range of samples and is clearly superior for the case of missing data.

**3.4 Remarks** 

*Ai*(*t*) is the amplitude of the phasor quantity representing the *i*th harmonic at time t

*i* is the phase angle of the *i*th harmonic relative to a reference rotating at *i*

*n* is the harmonic order

Each frequency component requires two state variables. Thus the total number of state variable is 2n. These state variables are defined as follows

Electric Power Systems Harmonics - Identification and Measurements 23

 1 1 1

1 2

1 1 1 2 2 2

1 0

1 1

<sup>0</sup> 1 1 1 1 <sup>0</sup>

1 1

cos sin sin cos *<sup>i</sup>*

*iw t iw t*

*iw t iw t* 

> *Xk kXk wk* 1

This model has constant state transition and measurement matrices. However, it assumes a stationary reference. Thus, the in-phase and quadrature phase components represent the

instantaneous values of con-sinusoidal and sinusoidal waveforms, respectively.

*x k wt wt xk wk x k wt wt xk wk* 

1 cos sin *kk k s t A t wt w t x k x k x k wt x k wt*

and *x k* <sup>2</sup> to be *A t wt k k* sin

   

1 2

 

 *n n*

(48)

*Z k HX k V k* (49)

2 1 2

cos 1

<sup>1</sup>

 

> 

2

*n*

If the signal includes *n* frequencies; the fundamental plus *n* – 1 harmonics, the state variable

*x k Z k v k x k* 

sin cos

*x k wt x k wt*

. At *tk*+2, which

(44)

(45)

(46)

*w k*

, *i n* 1, , (47)

11 2

1 sin

1 cos sin 1 sin cos

1 1 1 2 2

*x k x k <sup>M</sup> x k x k*

2 1 2 1 2 2

*x k x k <sup>M</sup> x k x k* 

*n n*

*n n*

*k k x k A t wt w t*

2 1

Thus, the state variable representation takes the following form

 

and the measurement equation then becomes

representation may be expressed as

 

where the sub-matrices *Mi* are given as

Equation (46) can be rewritten as

while equation (45) as

 

*M*

Now, consider *x k A t wt* <sup>1</sup> *k k* cos

also

is *tk* + *t*, the signal may be expressed as

$$\begin{aligned} \mathbf{x}\_1(t) &= A\_1(t) \cos \theta\_1, & \quad \mathbf{x}\_2(t) &= A\_1(t) \sin \theta\_1 \\ \mathbf{x}\_2(t) &= A\_2(t) \cos \theta\_2, & \quad \mathbf{x}\_3(t) &= A \mathbf{2}\_1(t) \sin \theta\_2 \\ &\dots & \dots & \dots \\ &\dots & \dots & \dots \\ \mathbf{x}\_{2n-1}(t) &= A\_n(t) \cos \theta\_{n'} & \quad \mathbf{x}\_{2n}(t) &= A\_2(t) \sin \theta\_n \end{aligned} \tag{38}$$

These state variables represent the in-phase and quadrate phase components of the harmonics with respect to a rotting reference, respectively. This may be referred to as model 1. Thus, the state variable equations may be expressed as:

$$
\begin{bmatrix} x\_1 \\ x\_2 \\ \vdots \\ x\_{2n-1} \\ x\_{2n} \end{bmatrix} = \begin{bmatrix} 1 & 0 & \dots & \dots & 0 \\ 0 & 1 & \dots & \dots & 0 \\ \vdots & \vdots & \vdots & \vdots & \vdots \\ 0 & 0 & \dots & 1 & 0 \\ 0 & 0 & \dots & \dots & 1 \end{bmatrix} \begin{bmatrix} x\_1 \\ x\_2 \\ \vdots \\ x\_{2n-1} \\ x\_{2n} \end{bmatrix} + \begin{bmatrix} a\_1 \\ a\_2 \\ \vdots \\ a\_{2n-1} \\ a\_{2n} \end{bmatrix} w\_k \tag{39}
$$

or in short hand

$$
\underline{X}(k+1) = \phi \underline{X}(k) + \underline{w}(k) \tag{40}
$$

where

*X* is a 2*n* 1 state vector

Is a 2*n* 2*n* state identity transition matrix, which is a diagonal matrix

*w*(*k*) is a 2*n* 1 noise vector associated with the transition of a sate from *k* to *k* + 1 instant The measurement equation for the voltage or current signal, in this case, can be rewritten as, equation (37)

$$\mathbf{s}\begin{pmatrix} k\Delta t \end{pmatrix} = \begin{bmatrix} \cos wk\Delta t & \sin wk\Delta t & \dots & \cos mwk\Delta t & \sin \begin{pmatrix} \mathbf{x}\_1 \\ \mathbf{x}\_2 \\ \vdots \\ \mathbf{x}\_{2n-1} \\ \mathbf{x}\_{2n} \end{pmatrix} + v\begin{pmatrix} k \end{pmatrix} \\ \text{(41)}$$

which can be written as

$$
\underline{Z}(k) = H(k)\underline{X}(k) + \underline{v}(k)\tag{42}
$$

where *Z*(k) is *m* 1 vector of measurements of the voltage or current waveforms, *H*(*k*) is *m* 2*n* measurement matrix, which is a time varying matrix and *v*(*k*) is *m* 1 errors measurement vector. Equation (40) and (42) are now suitable for Kalman filter application. Another model can be derived of a signal with time-varying magnitude by using a stationary reference, model 2. Consider the noise free signal to be

$$s(t\_i) = A(t\_i)\cos(wt + \theta) \tag{43}$$

Now, consider *x k A t wt* <sup>1</sup> *k k* cos and *x k* <sup>2</sup> to be *A t wt k k* sin . At *tk*+2, which is *tk* + *t*, the signal may be expressed as

$$\begin{aligned} \mathbf{s}\left(t\_{k+1}\right) &= A\left(t\_{k+1}\right)\cos\left(w t\_k + w\Delta t + \theta\right) = \mathbf{x}\_1\left(k+1\right) \\ \mathbf{x}\_1\left(k+1\right) &= \mathbf{x}\_1\left(k\right)\cos\left(w\Delta t\right) - \mathbf{x}\_2\left(k\right)\sin\left(w\Delta t\right) \end{aligned}$$

also

22 Power Quality Harmonics Analysis and Real Measurements Data

 

> 

These state variables represent the in-phase and quadrate phase components of the harmonics with respect to a rotting reference, respectively. This may be referred to as model

> *Xk Xk wk* 1

*w*(*k*) is a 2*n* 1 noise vector associated with the transition of a sate from *k* to *k* + 1 instant The measurement equation for the voltage or current signal, in this case, can be rewritten as,

where *Z*(k) is *m* 1 vector of measurements of the voltage or current waveforms, *H*(*k*) is *m* 2*n* measurement matrix, which is a time varying matrix and *v*(*k*) is *m* 1 errors measurement vector. Equation (40) and (42) are now suitable for Kalman filter application. Another model can be derived of a signal with time-varying magnitude by using a

cos *k k s t A t wt*

stationary reference, model 2. Consider the noise free signal to be

*s k t wk t wk t nwk t nwk t v k*

cos sin cos sin

*xt At xt At xt At xt A t*

 

11 1 21 1 22 2 3 1 2

cos , sin cos , 2 sin

*x t At* 2 2 *n n* sin

, (38)

1 2 

 

 

2 1

(40)

1 2

 

*x x*

> 2 1 2

(43)

*n n k*

 

*x x* 

(41)

*Zk HkXk vk* (42)

*n*

*wk*

(39)

2 2

*n n*

*x t At*

cos , *n nn*

1 1 1 0 0 0 1 0 2 2

*x x x x*

*x x n n x x*

 

0 0 1 0 2 1 2 1 0 0 1

Is a 2*n* 2*n* state identity transition matrix, which is a diagonal matrix

 

2 1

1. Thus, the state variable equations may be expressed as:

2

*X* is a 2*n* 1 state vector

or in short hand

where

equation (37)

which can be written as

*n*

 

$$\begin{aligned} \mathbf{x}\_2 \begin{pmatrix} k+1 \\ \end{pmatrix} &= A \begin{pmatrix} t\_{k+1} \\ \end{pmatrix} \sin \left( wt\_k + w\Delta t + \theta \right) \\ &= \mathbf{x}\_1 \begin{pmatrix} k \\ \end{pmatrix} \sin \left( w\Delta t \right) + \mathbf{x}\_2 \begin{pmatrix} k \\ \end{pmatrix} \cos \left( w\Delta t \right) \\ \end{aligned}$$

Thus, the state variable representation takes the following form

$$
\begin{bmatrix}
\mathbf{x}\_1(k+1) \\
\mathbf{x}\_2(k+1)
\end{bmatrix} = \begin{bmatrix}
\cos w\Delta t & -\sin w\Delta t \\
\sin w\Delta t & \cos w\Delta t
\end{bmatrix} \begin{bmatrix}
\mathbf{x}\_1(k) \\
\mathbf{x}\_2(k)
\end{bmatrix} + \begin{bmatrix}
w\_1(k) \\
w\_2(k)
\end{bmatrix} \tag{44}
$$

and the measurement equation then becomes

$$Z(k) = \begin{bmatrix} 1 & 0 \end{bmatrix} \begin{bmatrix} \mathbf{x}\_1(k) \\ \mathbf{x}\_2(k) \end{bmatrix} + v(k) \tag{45}$$

If the signal includes *n* frequencies; the fundamental plus *n* – 1 harmonics, the state variable representation may be expressed as

$$
\begin{bmatrix}
\mathbf{x}\_1(k+1) \\
\mathbf{x}\_2(k+1) \\
\vdots \\
\mathbf{x}\_{2n-1}(k+1) \\
\mathbf{x}\_{2n}(k+1) \\
\mathbf{x}\_{2n}(k+1)
\end{bmatrix} = \begin{bmatrix}
M\_1 & \dots & \dots & 0 \\
\dots & \dots & \dots & \dots \\
\dots & \dots & \dots & \dots \\
0 & \dots & \dots & M\_n
\end{bmatrix} \begin{bmatrix}
\mathbf{x}\_1(k+1) \\
\mathbf{x}\_2(k+1) \\
\vdots \\
\mathbf{x}\_{2n}(k+1) \\
\mathbf{x}\_{2n}(k+1)
\end{bmatrix} + \begin{bmatrix}
a\_1 \\
a\_2 \\
\vdots \\
a\_{2n-1} \\
a\_{2n}
\end{bmatrix} w(k) \tag{46}
$$

where the sub-matrices *Mi* are given as

$$M\_i = \begin{bmatrix} \cos iw\Delta t & -\sin iw\Delta t \\ \sin iw\Delta t & \cos iw\Delta t \end{bmatrix} \\ \text{/} i = 1, \dots, n \tag{47}$$

Equation (46) can be rewritten as

$$
\underline{X}(k+1) = \phi(k)\underline{X}(k) + \underline{w}(k)\tag{48}
$$

while equation (45) as

$$
\underline{X}(k) = H\underline{X}(k) + \underline{V}(k)\tag{49}
$$

This model has constant state transition and measurement matrices. However, it assumes a stationary reference. Thus, the in-phase and quadrature phase components represent the instantaneous values of con-sinusoidal and sinusoidal waveforms, respectively.

Electric Power Systems Harmonics - Identification and Measurements 25

Fig. 16. The first and second diagonal elements of *Pk* matrix using the 14-state model 1.

While the testing results of model 2 are given in Figures 26 –28. Figure 26 shows the first two components of Kalman gain vector. Figure 27 shows the first and second diagonal elements of *Pk*. The estimation of the magnitude of and third harmonic were exactly the

Fig. 15. Kalman gain for *x*1 and *x*2 using the 14-state model 1.

Fig. 17. Kalman gain for *x*1 and *x*2 using the 14-state model 2.

same as those shown in Figure 23.

#### **4.1 Testing the kalman filter algorithm**

The two Kalman filter models described in the preceding section were tested using a waveform with known harmonic contents. The waveform consists of the fundamental, the third, the fifth, the ninth, the eleventh, the thirteenth, and the nineteenth harmonics. The waveform is described as

$$\begin{split} s(t) &= 1.0 \cos \left( \alpha t + 10^{\circ} \right) + 0.1 \cos \left( 3 \alpha t + 20^{\circ} \right) + 0.08 \cos \left( 5 \alpha t + 30^{\circ} \right) \\ &+ 0.08 \cos \left( 7 \alpha t + 40^{\circ} \right) + 0.06 \cos \left( 11 \alpha t + 50^{\circ} \right) \\ &+ 0.05 \cos \left( 13 \alpha t + 60^{\circ} \right) + 0.03 \cos \left( 19 \alpha t + 70^{\circ} \right) \end{split}$$

The sampling frequency was selected to be 64 60 Hz.

i. Initial process vector

As the Kalman filter model started with no past measurement, the initial process vector was selected to be zero. Thus, the first half cycle (8 milliseconds) is considered to be the initialization period.


The noise variance was selected to be constant at a value of 0.05 p.u.2. This was passed on the background noise variance at field measurement.

iv. State variable covariance matrix (*Q*)

The matrix *Q* was also selected to be 0.05 p.u.

Testing results of model 1, which is a 14-state model described by equations (40) and (42) are given in the following figures. Figure 14 shows the initialization period and the recursive estimation of the magnitude of the fundamental and third harmonic. Figure 15 shows the Kalman gain for the fundamental component. Figure 16 shows the first and second diagonal element of *Pk*.

Fig. 14. Estimated magnitudes of 60 Hz and third harmonic component using the 14-state model 1.

The two Kalman filter models described in the preceding section were tested using a waveform with known harmonic contents. The waveform consists of the fundamental, the third, the fifth, the ninth, the eleventh, the thirteenth, and the nineteenth harmonics. The

1.0cos 10 0.1cos 3 20 0.08cos 5 30

As the Kalman filter model started with no past measurement, the initial process vector was selected to be zero. Thus, the first half cycle (8 milliseconds) is considered to be the

The initial covariance matrix was selected to be a diagonal matrix with the diagonal

The noise variance was selected to be constant at a value of 0.05 p.u.2. This was passed

Testing results of model 1, which is a 14-state model described by equations (40) and (42) are given in the following figures. Figure 14 shows the initialization period and the recursive estimation of the magnitude of the fundamental and third harmonic. Figure 15 shows the Kalman gain for the fundamental component. Figure 16 shows the first and second diagonal

Fig. 14. Estimated magnitudes of 60 Hz and third harmonic component using the 14-state

 

   

*t t*

*st t t t*

0.08cos 7 40 0.06 cos 11 50

on the background noise variance at field measurement.

The sampling frequency was selected to be 64 60 Hz.

0.05cos 13 60 0.03cos 19 70

*t t*

**4.1 Testing the kalman filter algorithm** 

waveform is described as

i. Initial process vector

initialization period. ii. Initial covariance matrix

values equal to 10 p.u.

iv. State variable covariance matrix (*Q*) The matrix *Q* was also selected to be 0.05 p.u.

iii. Noise variance (*R*)

element of *Pk*.

model 1.

Fig. 15. Kalman gain for *x*1 and *x*2 using the 14-state model 1.

Fig. 16. The first and second diagonal elements of *Pk* matrix using the 14-state model 1.

While the testing results of model 2 are given in Figures 26 –28. Figure 26 shows the first two components of Kalman gain vector. Figure 27 shows the first and second diagonal elements of *Pk*. The estimation of the magnitude of and third harmonic were exactly the same as those shown in Figure 23.

Fig. 17. Kalman gain for *x*1 and *x*2 using the 14-state model 2.

Electric Power Systems Harmonics - Identification and Measurements 27

The Kalman filter, however, can be applied for any number of samples over a half cycle. If the harmonic has time-varying magnitude, the Kalman filter algorithm would track the time variation after the initialization period (half a cycle). Figures 19 and 20 show the three-phase current and voltage waveforms recorded at the industrial load. Figures 21 –23 show the recursive estimation of the magnitude of the fundamental, fifth, and seventh harmonics; the eleventh and thirteenth harmonics; and the seventeenth and nineteenth harmonics, respectively, for phase A current. The same harmonic analysis was also applied to the actual recorded voltage waveforms. Figure 24 shows the recursive estimation of the magnitude of the fundamental and fifth harmonic for phase A voltage. No other voltage harmonics are

Fig. 19. Actual recorded current waveform of phases A, B, and C.

Fig. 20. Actual recorded voltage waveform of phase A, B, and C.

shown here due tot he negligible small value.

Fig. 18. The first and second diagonal elements of *Pk* matrix using the 14-state model 2.

The Kalman gain vector *Kk* and the covariance matrix *Pk* reach steady-state in about half a cycle, when model 1 is used, 1/60 seconds. Its variations include harmonics of 60 Hz. The covariance matrix in the steady-state consists of a constant plus a periodic component. These time variations are due to the time-varying vector in the measurement equation. Thus, after initialization of the model, the Kalman gain vector of the third cycle can be repeated for successive cycles.

When model 2 is selected, the components of the Kalman gain vector and the covariance matrix become constants. In both models, the Kalman gain vector is independent of the measurements and can be computed off-line. As the state transition matrix is a full matrix, it requires more computation than model 1 to update the state vector.

Kalman filter algorithm is also tested for actual recorded data. Two cases of actual recorded data are reported here. The first case represents a large industrial load served by two parallel transformers totaling 7500 KVA [5]. The load consists of four production lines of induction heating with two single-phase furnaces per line. The induction furnaces operate at 8500 Hz and are used to heat 40-ft steel rods which are cut into railroad spikes. Diodes are used in the rectifier for converting the 60 Hz power into dc and SCRs are used in the inverter for converting the dc into single-phase 8500 Hz power. The waveforms were originally sampled at 20 kHz. A program was written to use a reduced sampling rate in the analysis. A careful examination of the current and voltage waveforms indicated that the waveforms consist of (1) harmonics of 60 Hz and (2) a decaying periodic high-frequency transients. The high-frequency transients were measured independently for another purpose [6]. The rest of the waveform was then analyzed for harmonic analysis. Using a sampling frequency that is a multiple of 2 kHz, the DFT was then applied for a period of 3 cycles. The DFT results were as follows:


Fig. 18. The first and second diagonal elements of *Pk* matrix using the 14-state model 2.

successive cycles.

The Kalman gain vector *Kk* and the covariance matrix *Pk* reach steady-state in about half a cycle, when model 1 is used, 1/60 seconds. Its variations include harmonics of 60 Hz. The covariance matrix in the steady-state consists of a constant plus a periodic component. These time variations are due to the time-varying vector in the measurement equation. Thus, after initialization of the model, the Kalman gain vector of the third cycle can be repeated for

When model 2 is selected, the components of the Kalman gain vector and the covariance matrix become constants. In both models, the Kalman gain vector is independent of the measurements and can be computed off-line. As the state transition matrix is a full matrix, it

Kalman filter algorithm is also tested for actual recorded data. Two cases of actual recorded data are reported here. The first case represents a large industrial load served by two parallel transformers totaling 7500 KVA [5]. The load consists of four production lines of induction heating with two single-phase furnaces per line. The induction furnaces operate at 8500 Hz and are used to heat 40-ft steel rods which are cut into railroad spikes. Diodes are used in the rectifier for converting the 60 Hz power into dc and SCRs are used in the inverter for converting the dc into single-phase 8500 Hz power. The waveforms were originally sampled at 20 kHz. A program was written to use a reduced sampling rate in the analysis. A careful examination of the current and voltage waveforms indicated that the waveforms consist of (1) harmonics of 60 Hz and (2) a decaying periodic high-frequency transients. The high-frequency transients were measured independently for another purpose [6]. The rest of the waveform was then analyzed for harmonic analysis. Using a sampling frequency that is a multiple of 2 kHz, the DFT was then applied for a period of 3 cycles. The DFT results were as follows:

> Freq. (Hz) Mag. Angle (rad.) 60 1.0495 -0.20 300 0.1999 1.99 420 0.0489 -2.18 660 0.0299 0.48 780 0.0373 2.98 1020 0.0078 -0.78 1140 0.0175 1.88

requires more computation than model 1 to update the state vector.

Fig. 19. Actual recorded current waveform of phases A, B, and C.

The Kalman filter, however, can be applied for any number of samples over a half cycle. If the harmonic has time-varying magnitude, the Kalman filter algorithm would track the time variation after the initialization period (half a cycle). Figures 19 and 20 show the three-phase current and voltage waveforms recorded at the industrial load. Figures 21 –23 show the recursive estimation of the magnitude of the fundamental, fifth, and seventh harmonics; the eleventh and thirteenth harmonics; and the seventeenth and nineteenth harmonics, respectively, for phase A current. The same harmonic analysis was also applied to the actual recorded voltage waveforms. Figure 24 shows the recursive estimation of the magnitude of the fundamental and fifth harmonic for phase A voltage. No other voltage harmonics are shown here due tot he negligible small value.

Fig. 20. Actual recorded voltage waveform of phase A, B, and C.

Electric Power Systems Harmonics - Identification and Measurements 29

Fig. 24. Estimated magnitudes of the 60 Hz and fifth harmonic for phase A voltage.

continuously varying. The total harmonic distortion experienced similar variation.

Fig. 25. Current waveform of a continuous varying load.

can track harmonics with the time varying magnitudes.

need to be correctly defined.

The second case represents a continuous dynamic load. The load consists of two six-phase drives for two 200 HP dc motors. The current waveform of one phase is shown in Figure 25. The harmonic analysis using the Kalman filter algorithm is shown in Figure 35. It should be noted that the current waveform was continuously varying in magnitude due to the dynamic nature of the load. Thus, the magnitude of the fundamental and harmonics were

There is no doubt that the Kalman filtering algorithm is more accurate and is not sensitive to a certain sampling frequency. As the Kalman filter gain vector is time0varying, the estimator

Two models are described in this section to show the flexibility in the Kalman filtering scheme. There are many applications, where the results of FFT algorithms are as accurate as a Kalman filter model. However, there are other applications where a Kalman filter becomes superior to other algorithms. Implementing linear Kalman filter models is relatively a simple task. However, state equations, measurement equations, and covariance matrices

Fig. 21. Estimated magnitudes of the fundamental, fifth, and seventh harmonics for phase A current.

Fig. 22. Estimated magnitudes of the eleventh and thirteenth harmonics for phase A current.

Fig. 23. Estimated magnitudes of the seventeenth and nineteenth harmonics for phase A current.

Fig. 21. Estimated magnitudes of the fundamental, fifth, and seventh harmonics for phase A

Fig. 22. Estimated magnitudes of the eleventh and thirteenth harmonics for phase A current.

Fig. 23. Estimated magnitudes of the seventeenth and nineteenth harmonics for phase A

current.

current.

Fig. 24. Estimated magnitudes of the 60 Hz and fifth harmonic for phase A voltage.

The second case represents a continuous dynamic load. The load consists of two six-phase drives for two 200 HP dc motors. The current waveform of one phase is shown in Figure 25. The harmonic analysis using the Kalman filter algorithm is shown in Figure 35. It should be noted that the current waveform was continuously varying in magnitude due to the dynamic nature of the load. Thus, the magnitude of the fundamental and harmonics were continuously varying. The total harmonic distortion experienced similar variation.

Fig. 25. Current waveform of a continuous varying load.

There is no doubt that the Kalman filtering algorithm is more accurate and is not sensitive to a certain sampling frequency. As the Kalman filter gain vector is time0varying, the estimator can track harmonics with the time varying magnitudes.

Two models are described in this section to show the flexibility in the Kalman filtering scheme. There are many applications, where the results of FFT algorithms are as accurate as a Kalman filter model. However, there are other applications where a Kalman filter becomes superior to other algorithms. Implementing linear Kalman filter models is relatively a simple task. However, state equations, measurement equations, and covariance matrices need to be correctly defined.

Electric Power Systems Harmonics - Identification and Measurements 31

recorded data set. and computes the voltage and current harmonics magnitude, the voltage and current harmonics phase angles, and the fundamental power and harmonics power.

To initialize the recursive process of the proposed filter, with an initial process vector and covariance matrix *P*, a simple deterministic procedure uses the static least squares error estimate of previous measurements. Thus, the initial process vector may be computed as:

> <sup>ˆ</sup> *T T X HH Hz*

> > <sup>ˆ</sup> *<sup>T</sup> P HH*

where *H* is an *m m* matrix of measurements, and *z* is an *m* 1 vector of previous measurements, the initial process vector may be selected to be zero, and the first few

The proposed algorithm and the two models were tested using a voltage signal waveform of

 1cos 10 0.1cos 3 20 0.08cos 5 30 0.08cos 9 40

 

*f*, *f* = 60 + *f*). Figs. 24 and 29 give the results obtained for these two

 

*ttt*

The data window size is two cycles, with sampling frequency of 64 samples/cycle. That is, the total number of samples used is 128 samples, and the sampling frequency is 3840 Hz.

Using the two models, the proposed filtering algorithm estimates exactly the harmonic content of the voltage waveform both magnitudes and phase angles and the two proposed

The steady-state gain of the proposed filter is periodic with a period of 1/60 s. This time variation is due to the time varying nature of the vector states in the measurement equation.

The gain of the proposed filter reaches the steady-state value in a very short time, since the initialization of the recursive process, as explained in the preceding section, was sufficiently

The effects of frequency drift on the estimate are also considered. We assume small and large values for the frequency drift: *f* = -0.10 Hz and *f* = -1.0 Hz, respectively. In this study the elements of the matrix *H*(*k*) are calculated at 60 Hz, and the voltage signal is

frequency deviations for the fundamental and the third harmonic. Fig. 55 gives the estimated magnitude, and Fig. 29 gives the estimated phase angles. Examination of these

0.06cos 11 50 0.05cos 13 60 0.03cos 19 70 *v tt t t t*

0

0

and the corresponding covariance error matrix is:

milliseconds are considered to be the initialization period.

For this simulated example we have the following results.

Figure 54 give the proposed filter gain for *X*1 and *Y*1.

**4.3 Testing the algorithm using simulated data** 

known harmonic contents described as:

models produce the same results.

accurate.

sampled at (

 = 2

two curves reveals the following:

1

1

*Initialization of the filter* 

Kalman filter used in the previous section assumes that the digital samples for the voltage and current signal waveforms are known in advance, or at least, when it is applied on-line, good estimates for the signals parameters are assumed with a certain degree of accuracy, so that the filter converges to the optimal estimates in few samples later. Also, it assumes that an accurate model is presented for the signals; otherwise inaccurate estimates would be obtained. Ref. 8 uses the Kalman filter algorithm to obtain the optimal estimate of the power system harmonic content. The measurements used in this reference are the power system voltage and line flows at different harmonics obtained from a harmonic load flow program (HARMFLO). The effect of load variation over a one day cycle on the power system harmonics and standard are presented. The optimal estimates, in this reference, are the power system bus voltage magnitudes and phase angles at different harmonic level.

Fig. 35. Magnitude of dominant frequencies and harmonic distortion of waveform shown in Figure 34 using the Kalman filtering approach.

#### **4.2 Linear dynamic weighted least absolute estimates [11]**

This section presents the application of the linear dynamic weighted least absolute value dynamic filter for power system harmonics identification and measurements. The two models developed earlier, model 1 and model 2, are used with this filter. As we explained earlier, this filter can deal easily with the outlier, unusual events, in the voltage or current waveforms.

#### *Software implementation*

A software package has been developed to analyze digitized current and voltage waveforms. This package has been tested on simulated data sets, as well as on an actual recorded data set. and computes the voltage and current harmonics magnitude, the voltage and current harmonics phase angles, and the fundamental power and harmonics power.

#### *Initialization of the filter*

30 Power Quality Harmonics Analysis and Real Measurements Data

Kalman filter used in the previous section assumes that the digital samples for the voltage and current signal waveforms are known in advance, or at least, when it is applied on-line, good estimates for the signals parameters are assumed with a certain degree of accuracy, so that the filter converges to the optimal estimates in few samples later. Also, it assumes that an accurate model is presented for the signals; otherwise inaccurate estimates would be obtained. Ref. 8 uses the Kalman filter algorithm to obtain the optimal estimate of the power system harmonic content. The measurements used in this reference are the power system voltage and line flows at different harmonics obtained from a harmonic load flow program (HARMFLO). The effect of load variation over a one day cycle on the power system harmonics and standard are presented. The optimal estimates, in this reference, are the

power system bus voltage magnitudes and phase angles at different harmonic level.

Fig. 35. Magnitude of dominant frequencies and harmonic distortion of waveform shown in

This section presents the application of the linear dynamic weighted least absolute value dynamic filter for power system harmonics identification and measurements. The two models developed earlier, model 1 and model 2, are used with this filter. As we explained earlier, this filter can deal easily with the outlier, unusual events, in the voltage or current

A software package has been developed to analyze digitized current and voltage waveforms. This package has been tested on simulated data sets, as well as on an actual

Figure 34 using the Kalman filtering approach.

waveforms.

*Software implementation* 

**4.2 Linear dynamic weighted least absolute estimates [11]** 

To initialize the recursive process of the proposed filter, with an initial process vector and covariance matrix *P*, a simple deterministic procedure uses the static least squares error estimate of previous measurements. Thus, the initial process vector may be computed as:

$$\hat{X}\_0 = \left[H^\top H\right]^{-1} H^\top z$$

and the corresponding covariance error matrix is:

$$
\hat{P}\_0 = \left[H^T H\right]^{-1}
$$

where *H* is an *m m* matrix of measurements, and *z* is an *m* 1 vector of previous measurements, the initial process vector may be selected to be zero, and the first few milliseconds are considered to be the initialization period.

#### **4.3 Testing the algorithm using simulated data**

The proposed algorithm and the two models were tested using a voltage signal waveform of known harmonic contents described as:

$$v(t) = 1\cos\left(\alpha t + 10^{\circ}\right) + 0.1\cos\left(3\alpha t + 20^{\circ}\right) + 0.08\cos\left(5\alpha t + 30^{\circ}\right) + 0.08\cos\left(9\alpha t + 40^{\circ}\right)$$

$$+ 0.06\cos\left(11\alpha t + 50^{\circ}\right) + 0.05\cos\left(13\alpha t + 60^{\circ}\right) + 0.03\cos\left(19\alpha t + 70^{\circ}\right)$$

The data window size is two cycles, with sampling frequency of 64 samples/cycle. That is, the total number of samples used is 128 samples, and the sampling frequency is 3840 Hz. For this simulated example we have the following results.

Using the two models, the proposed filtering algorithm estimates exactly the harmonic content of the voltage waveform both magnitudes and phase angles and the two proposed models produce the same results.

The steady-state gain of the proposed filter is periodic with a period of 1/60 s. This time variation is due to the time varying nature of the vector states in the measurement equation. Figure 54 give the proposed filter gain for *X*1 and *Y*1.

The gain of the proposed filter reaches the steady-state value in a very short time, since the initialization of the recursive process, as explained in the preceding section, was sufficiently accurate.

The effects of frequency drift on the estimate are also considered. We assume small and large values for the frequency drift: *f* = -0.10 Hz and *f* = -1.0 Hz, respectively. In this study the elements of the matrix *H*(*k*) are calculated at 60 Hz, and the voltage signal is sampled at ( = 2*f*, *f* = 60 + *f*). Figs. 24 and 29 give the results obtained for these two frequency deviations for the fundamental and the third harmonic. Fig. 55 gives the estimated magnitude, and Fig. 29 gives the estimated phase angles. Examination of these two curves reveals the following:

Electric Power Systems Harmonics - Identification and Measurements 33

Fig. 29. Estimated phase angles for frequency drifts using models 1 and 2

since the filter gain *K*(*k*) does not depend on the measurements (eqn. 8).

the estimation process is begun.

**4.4 Testing on actual recorded data** 

To overcome this drawback, it has been found through extensive runs that if the elements of the matrix *H*(*k*) are calculated at the same frequency of the voltage signal waveform, good estimates are produced and the frequency drift has in this case no effect. Indeed, to perform this modification the proposed algorithm needs a frequency-measurement algorithm before

It has been found, through extensive runs that the filter gains for the fundamental voltage components, as a case study, do not change with the frequency drifts. Indeed, that is true

As the state transition matrix for model 2 is a full matrix, it requires more computation than model 1 to update the state vector. Therefore in the rest of this study, only model 1 is used.

The proposed algorithm is implemented to identify and measure the harmonics content for a practical system of operation. The system under study consists of a variable-frequency drive that controls a 3000 HP, 23 kV induction motor connected to an oil pipeline compressor. The waveforms of the three phase currents are given in Fig. 31. It has been found for this system that the waveforms of the phase voltages are nearly pure sinusoidal waveforms. A careful examination of the current waveforms revealed that the waveforms consist of: harmonics of 60 Hz, decaying period high-frequency transients, and harmonics of less than 60 Hz (sub-harmonics). The waveform was originally sampled at a 118 ms time

quality.

drift both phase angles have large changes and the estimates produced are of bad

Fig. 27. Gain of the proposed filter for *X*1 and *Y*1 using models 1 and 2.

Fig. 28. Estimated magnitudes of 60 Hz and third harmonic for frequency drifts using models 1 and 2.


Fig. 27. Gain of the proposed filter for *X*1 and *Y*1 using models 1 and 2.

Fig. 28. Estimated magnitudes of 60 Hz and third harmonic for frequency drifts using

 For a small frequency drift, *f* = -0.10 Hz, the fundamental magnitude and the third harmonic magnitude do not change appreciably; whereas for a large frequency drift, *f* = -1.0 Hz, they exhibit large relative errors, ranging from 7% for the fundamental to 25%

 On the other hand, for the small frequency drift the fundamental phase angle and the third harmonic phase angle do not change appreciably, whereas for the large frequency

models 1 and 2.

for the third harmonics.

drift both phase angles have large changes and the estimates produced are of bad quality.

Fig. 29. Estimated phase angles for frequency drifts using models 1 and 2

To overcome this drawback, it has been found through extensive runs that if the elements of the matrix *H*(*k*) are calculated at the same frequency of the voltage signal waveform, good estimates are produced and the frequency drift has in this case no effect. Indeed, to perform this modification the proposed algorithm needs a frequency-measurement algorithm before the estimation process is begun.

It has been found, through extensive runs that the filter gains for the fundamental voltage components, as a case study, do not change with the frequency drifts. Indeed, that is true since the filter gain *K*(*k*) does not depend on the measurements (eqn. 8).

As the state transition matrix for model 2 is a full matrix, it requires more computation than model 1 to update the state vector. Therefore in the rest of this study, only model 1 is used.

#### **4.4 Testing on actual recorded data**

The proposed algorithm is implemented to identify and measure the harmonics content for a practical system of operation. The system under study consists of a variable-frequency drive that controls a 3000 HP, 23 kV induction motor connected to an oil pipeline compressor. The waveforms of the three phase currents are given in Fig. 31. It has been found for this system that the waveforms of the phase voltages are nearly pure sinusoidal waveforms. A careful examination of the current waveforms revealed that the waveforms consist of: harmonics of 60 Hz, decaying period high-frequency transients, and harmonics of less than 60 Hz (sub-harmonics). The waveform was originally sampled at a 118 ms time

Electric Power Systems Harmonics - Identification and Measurements 35

Fig. 31. Estimated fundamental voltage.

Fig. 32. Estimated voltage harmonics for *V*

interval and a sampling frequency of 8.5 kHz. A computer program was written to change this sampling rate in the analysis.

Figs. 31 and 32 show the recursive estimation of the magnitude of the fundamental, second, third and fourth harmonics for the voltage of phase *A*. Examination of these curves reveals that the highest-energy harmonic is the fundamental, 60 Hz, and the magnitude of the second, third and fourth harmonics are very small. However, Fig. 33 shows the recursive estimation of the fundamental, and Fig. 34 shows the recursive estimation of the second, fourth and sixth harmonics for the current of phase A at different data window sizes. Indeed, we can note that the magnitudes of the harmonics are time-varying since their magnitudes change from one data window to another, and the highest energy harmonics are the fourth and sixth. On the other hand, Fig. 35 shows the estimate of the phase angles of the second, fourth and sixth harmonics, at different data window sizes. It can be noted from this figure that the phase angles are also time0varing because their magnitudes vary from one data window to another.

Fig. 30. Actual recorded current waveform of phases *A*, *B* and *C*.

Fig. 31. Estimated fundamental voltage.

interval and a sampling frequency of 8.5 kHz. A computer program was written to change

Figs. 31 and 32 show the recursive estimation of the magnitude of the fundamental, second, third and fourth harmonics for the voltage of phase *A*. Examination of these curves reveals that the highest-energy harmonic is the fundamental, 60 Hz, and the magnitude of the second, third and fourth harmonics are very small. However, Fig. 33 shows the recursive estimation of the fundamental, and Fig. 34 shows the recursive estimation of the second, fourth and sixth harmonics for the current of phase A at different data window sizes. Indeed, we can note that the magnitudes of the harmonics are time-varying since their magnitudes change from one data window to another, and the highest energy harmonics are the fourth and sixth. On the other hand, Fig. 35 shows the estimate of the phase angles of the second, fourth and sixth harmonics, at different data window sizes. It can be noted from this figure that the phase angles are also time0varing because their magnitudes vary from

this sampling rate in the analysis.

one data window to another.

Fig. 30. Actual recorded current waveform of phases *A*, *B* and *C*.

Fig. 32. Estimated voltage harmonics for *V*

Electric Power Systems Harmonics - Identification and Measurements 37

For this system the highest-energy harmonic component is the fundamental power, the

Fig. 35. Harmonics phase angles of *IA* against time steps at various window sizes.

power due to the fundamental voltage and current.

Fig. 36. Fundamental powers against time steps.

Fig. 33. Estimated fundamental current *IA*.

Fig. 34. Harmonics magnitude of *IA* against time steps at various window sizes.

Furthermore, Figs. 36 – 38 show the recursive estimation of the fundamental, fourth and the sixth harmonics power, respectively, for the system under study (the factor 2 in these figures is due to the fact that the maximum values for the voltage and current are used to calculate this power). Examination of these curves reveals the following results. The fundamental power and the fourth and sixth harmonics are time-varying.

Fig. 33. Estimated fundamental current *IA*.

Fig. 34. Harmonics magnitude of *IA* against time steps at various window sizes.

power and the fourth and sixth harmonics are time-varying.

Furthermore, Figs. 36 – 38 show the recursive estimation of the fundamental, fourth and the sixth harmonics power, respectively, for the system under study (the factor 2 in these figures is due to the fact that the maximum values for the voltage and current are used to calculate this power). Examination of these curves reveals the following results. The fundamental For this system the highest-energy harmonic component is the fundamental power, the power due to the fundamental voltage and current.

Fig. 35. Harmonics phase angles of *IA* against time steps at various window sizes.

Fig. 36. Fundamental powers against time steps.

Electric Power Systems Harmonics - Identification and Measurements 39

The proposed algorithm is compared with KF algorithm. Fig 39 gives the results obtained when both filters are implemented to estimate the second harmonic components of the current in phase *A*, at different data window sizes and when the considered number of harmonics is 15. Examination of the Figure reveals the following; both filters produce almost the same estimate for the second harmonic magnitude; and the magnitude of the

**4.5 Comparison with Kalman Filter (KF) algorithm** 

estimated harmonic varies from one data window to another.

Fig. 39. Estimated second harmonic magnitude using KF and WLAV.

2, and in the second Subsection the actual recorded data set is used.

magnitude instead of a constant magnitude.

In this Section the effects of outliers (unusual events on the system waveforms) are studied, and we compare the new proposed filter and the well-known Kalman filtering algorithm. In the first Subsection we compare the results obtained using the simulated data set of Section

The simulated data set of Section 4.3 has been used in this Section, where we assume (randomly) that the data set is contaminated with gross error, we change the sign for some measurements or we put these measurements equal to zero. Fig. 40 shows the recursive estimate of the fundamental voltage magnitude using the proposed filter and the wellknown Kalman filtering algorithm. Careful examination of this curve reveals the following

The proposed dynamic filter and the Kalman filter produce an optimal estimate to the fundamental voltage magnitude, depending on the data considered. In other words, the voltage waveform magnitude in the presence of outliers is considered as a time-varying

The proposed filter and the Kalman filter take approximately two cycles to reach the exact value of the fundamental voltage magnitude. However, if such outliers are corrected, the discrete least absolute value dynamic filter almost produces the exact value of the fundamental voltage during the recursive process, and the effects of the outliers are greatly

**4.5.1 Effects of outliers** 

**Simulated data** 

reduced Figure 41.

results.

Fig. 37. Fourth harmonic power in the three phases against time steps at various window sizes.

The fundamental powers, in the three phases, are unequal; i.e. the system is unbalanced. The fourth harmonic of phase *C*, and later after 1.5 cycles of phase *A*, are absorbing power from the supply, whereas those for phase *B* and the earlier phase *A* are supplying power to the network.

The sixth harmonic of phase *B* is absorbing power from the network, whereas the six harmonics of phases *A* and *C* are supplying power to the network; but the total power is still the sum of the three-phase power.

Fig. 38. Sixth harmonic powers in the three phases against time steps at various window sizes.

The fundamental power and the fourth and sixth harmonics power are changing from one data window to another.
