- Email: [email protected]

Simulated annealing optimization algorithm for power systems quality analysis S.A. Solimana,*, A.H. Mantawayb, M.E. El-Hawaryc a

Department of Electrical Engineering, University of Qatar, P.O. Box 2713, Doha, Qatar b KFUPM, Department of Electrical Engineering, Dhahran, Saudi Arabia c Department of Electrical and Computer Engineering, Dalhousie University , P.O. Box 1000, Halifax, NS, Canada

Abstract This article presents a new application of Simulated Annealing (SA) optimization algorithm for harmonics and frequency evaluation, for power system quality analysis and frequency relaying. An objective function based on the sum of the squares of the error is to be minimized to identify the amplitude and phase angle of each harmoic component as well as the fundamental frequency of the voltage signal. The problem formulated is a highly nonlinear optimization problem and does not need any approximation such as Taylor’s series expansion. The proposed algorithm uses the digitized samples of the voltage signal at the place where the power quality and frequency relaying are to be implemented. The proposed algorithm has an adaptive cooling schedule and a variable discretization to enhance the speed and convergence of the original SA algorithm. The proposed algorithm is tested on simulated data. Effects of different parameters on the estimated parameters are studied. It is shown that the proposed algorithm is able to identify harmonics spectrum in the signal. q 2004 Elsevier Ltd. All rights reserved. Keywords: Simulated annealing; Power quality; Frequency and harmonics measurements

1. Introduction Due to the widespread use of power electronic equipment in power systems, as well as the nonlinear loads, especially arc furnace makes the needs to adapt power quality criteria at the bus feeds such type of loads. Many digital algorithms have been developed for measuring harmonics and flicker in the network. Ref. [1] applies a multi-rate digital signal processing technique for harmonic analysis. An antialiasing filter is used with FFT to estimate the system harmonics. Ref. [2] proposes a scheme based on Parseval’s relation and energy concept, which defines a ‘group harmonics’ identification algorithm for the estimation of the energy distribution in the harmonics of time-varying waveforms. Using this approach many of the drawbacks of FFT are overcame. A Newton solution for the harmonic phasor analysis of ac/dc converter is proposed in Ref. [3]. The nonlinear equation for the voltage contaminated with harmonics is solved for the harmonics pharameters, using Newton iterations. The solution includes the interaction of the converter with the dc system. An experimental way to measure the contribution of City Street gas discharge lamps * Corresponding author. E-mail address: [email protected] (S.A. Soliman). 0142-0615/$ - see front matter q 2004 Elsevier Ltd. All rights reserved. doi:10.1016/S0142-0615(03)00068-1

is developed in Ref. [4]. The equipment used in this experiment is described. It is shown that an equivalent current harmonic spectrum exists which is independent of the type of gas and the lamp rated power. Ref. [5] applied a 12state Kalman filter to the voltage or current samples to obtain the instantaneous values for a maximum of six harmonics as well as the existing harmonic distortion in real time. Ref. [6] discusses problems associated with direct application of FFT to computer harmonic levels on nonsteady-state distorted waveforms, and various ways to describe recorded data in statistical terms. Fast Fourier Transform of a sequence which slides over a time-limited rectangular window is carried out in a nonrecursive manner to compute the harmonics for digital relaying is applied in Ref. [7], where only certain isolated harmonics, rather than the full spectrum are needed. The proposed algorithm is compared with the DFT and it is proven that this algorithm is less expensive than the DFT. The concept of a switching function for harmonics measurement is applied in Ref. [8] for studying ac and dc harmonics generation by 24-pulse phase controlled converters under unbalanced operation. The characteristic and noncharacteristic ac and dc harmonics using symmetrical voltage components and rectifier switching-functions theory is presented.

32

S.A. Soliman et al. / Electrical Power and Energy Systems 26 (2004) 31–36

Fuzzy linear regression for measurement of harmonics is presented in Ref. [9], where fuzzy parameters for each harmonic component are assumed, having certain middle and certain spread. The problem is converted to one of linear optimization problem, based on two inequality constraints to insure that the optimal solution is within the assumed membership. A digital power analyzer for use in power system for harmonics analyzing is described in Ref. [10]. The measurement algorithm is based on a recursive digital filter that processes input signal sample values. Ref. [11] presents a technique for harmonic extraction in a signal based on a modeling and identification method using recursive least-square error algorithm. The parameters of the harmonics are estimated from the uniformly sampled signal. A Hopfield neural network is used in Ref. [12] to determine simultaneously the supply-frequency variation, the fundamental-amplitude/phase variation as well as the harmonics-amplitude/phase variation. The nonlinear least error squares algorithm based on Taylor’s series expansion is used to compare with the Neural Network based Hopfield algorithm. A set of data taken on site was used. An enhanced FFT-based parametric (E-FFT) algorithm suitable for on-line harmonic analysis of electrical power systems is presented in Ref. [13]. The E-FFT algorithm exploits its iteration loops in combination with the characteristic of steep-descent gradient search strategy. The E-FFT algorithm performs reasonably well with short data record length. Ref. [14] offers the optimization of spectrum analysis for harmonics signal, by means of scale fine-tuning to overcome the drawbacks of FFT. These include the picked fence effect and the leakage effect. Simulated Annealing (SA) is a powerful technique for solving many optimization problems [16 –22]. It has the ability of escaping local minimal by incorporating a probability function in accepting or rejecting new solutions. It has been successfully used to solve combinatorial optimization problems, where the cost function is defined in a discrete domain. Recently, the SA algorithm has been successfully used to solve many important problems with functions of continuous variables, similar the problem of this study. The application of SA algorithm to such type of problems requires an efficient strategy to generate feasible trial solutions [16,22]. The other important point is the selection of an efficient cooling schedule that adapts itself according to statistics of the trial moves (acceptance/ rejection) during the search. This article presents a new application of SA algorithm for harmonics and frequency evaluation, for power quality analysis and frequency relaying. The proposed algorithm minimizes an objective function based on the sum of square of the errors to identify the amplitude, phase angle of each harmonic as well as the fundamental frequency. The proposed algorithm is tested on simulated and actual recorded data. Effects of different parameters on the estimated states are studied.

2. Problem formulation Consider the Fourier expansion of a current or voltage wave with N harmonic terms: yðtÞ ¼

N X

ðan sin 2pnft þ bn cos 2pnftÞ

n¼1

¼

N X

ð1Þ

An sinð2pnft þ un Þ

n¼1

Assume m samples of yðti Þ; i ¼ 1; 2; …; m are available. The objective is to estimate the values of An ; un and f so as to minimize the sum of the square of the error E given by: " #2 m N X X yðti Þ 2 ðAn sinð2pnfti þ un Þ ð2Þ E¼ i¼1

n¼1

In the above equation: n is the harmonic order, N is the largest harmonics order to be expected in the signal, An is the harmonic amplitude and un is the harmonic phase angle. Furthermore, f is the fundamental frequency of the signal. The objective function in Eq. (2) is a nonlinear function in the variables. In this article, two cases are studied. In the first case we assume that the signal fundamental frequency is known in advance and equals to nominal frequency 50 or 60 Hz. In the second case the signal frequency is assumed to be unknown. In Section 3, the SA algorithm to find the best estimate for the parameters in Eq. (2) is described.

3. The proposed SA algorithm The proposed algorithm aims to identify the amplitude and phase angle of each harmonic component as well as the fundamental frequency of a given signal. To obtain the optimal parameters estimates, the problem is formulated as one of nonlinear optimization problem. The objective criterion (function) (2), is chosen to minimize the sum of squares of the error between the sampled signal and the estimated one at all sampling time periods. The implementation details of the SAA are given in the following sections. 3.1. SA algorithm The major steps of the algorithm are summarized as follows: Step (0): Set iteration counter ITR ¼ 0 Set the initial temperature of the cooling schedule those results in high probability of accepting new solutions. Initialize the step size vector, Gstep ðiÞ for all variables, GðiÞ; i ¼ 1; 2; 3; …; N: Step (1): Find, randomly, initial values for the estimated parameters and set it as the current and best solution.

S.A. Soliman et al. / Electrical Power and Energy Systems 26 (2004) 31–36

Step (2): Determine the error (objective function) for the current estimated parameters. Step (3): Generate randomly a new estimate (new trial solution), as a neighbor to the current solution. Step (4): Calculate the objective function at the new estimate. Step (5): Perform the SA acceptance test; to accept or reject the trial solution (see Section 3.2). Step (6): Check for equilibrium at this temperature (see Section 3.3). Step (7): If the prespecified maximum number of iterations is reached then stop. Else go to step (8). Step (8): If the step size vector values for all variables ðGstep Þ are less than a prespecified value then stop, else go to step (9). Step (9): Update the step size vector values. (See Section 3.4). Decrease the temperature according to the polynomial time cooling schedule (see Section 3.5). Go to step (3). 3.2. SA test The implementation steps of the SA test as applied to each iteration in the algorithm are described as follows [16 and 22]. Step (1): At the same calculated temperature, Cpk apply the following acceptance test for the new trial solution. Step (2): Acceptance test: If Ej # Ei ; or If exp½Ei 2 Ej Þ=Cp $ R ð0; 1Þ; then accept the trial solution, set Xi to Xj and Ei ¼ Ej : Otherwise reject the trial solution. Where Xi ; Xj ; Ei ; Ej are the SA current solution, the trial solution and their corresponding cost, respectively. Step (3): Go to the next step in the algorithm. 3.3. Equilibrium test The sequence of trial solutions generated in the SAA at a fixed temperature is stopped as soon as thermodynamic equilibrium, detected by some adequate condition, is reached. Then the temperature and step vectors are suitably adjusted [16,22]. The test of equilibrium is done as follows: If the NTRACP ðTÞ , N1 £ n and NTR ðTÞ , N2 £ n then continue at the same temperature, otherwise end of temperature stage.

33

Where NTRACP ðTÞ; NTR ðTÞ are the number of trial accepted and attempted at temperature T, respectively. n is the number of variables in the problem. N1 and N2 are end temperature stage parameters (taken as 12 and 100) [22].

3.4. Step size vector adjustment In this work the step vector is updated jointly with cooling schedule temperature, according to the acceptance rate of the attempted moves at the previous temperature stage [16,22]. All its components are updates simultaneously. The following steps explain how the step vector is adjusted mathematically: Step (1): Calculate PðiÞ ¼ NTRACPðiÞ=NTRðiÞ; i ¼ 1; …; N: Where NTRACP ðiÞ is the number of trial accepted then variable i is changed. NTR ðiÞ is the number of trials attempted by changing variable i: Step (2): If PðiÞ . PMAX, then STEP(i) ¼ STEP(i) £ STEPMAX. If PðiÞ , PMIN, then STEP(i) ¼ STEP(i) £ STEPMIN. Where PMAX, PMIN, STEPMAX and STEPMIN are parameters taken in our implementation as 0.05, 0.5, 0.8 and 1.2, respectively [22]. 3.5. Cooling schedule A finite-time implementation of the SAA can be realized by generating homogenous Markov chains of finite length for a finite sequence of descending values of the control parameter. To achieve this, one must specify a set of parameters that governs the convergence of the algorithm. These parameters form a cooling schedule. The parameters of the cooling schedules are: an initial value of the control parameter decrement function for decreasing the control parameter and a final value of the control parameter specified by the stopping criterion, and a finite length of each homogenous Markov chain. In this work a polynomial-time cooling schedule is used in which the temperature is decreased based on the statistics of the trial solutions acceptance or rejection during the search. Details of the implemented cooling schedule are described in Refs. [17,18].

4. Numerical results 4.1. Signal with known frequency Example 1. In this example we assume the fundamental frequency of the signal equals 50 Hz, and the signal

34

S.A. Soliman et al. / Electrical Power and Energy Systems 26 (2004) 31–36

waveform contains the third, fifth, seventh and ninth harmonics. Thus the signal waveform can be written as: pﬃﬃ vðtÞ ¼ 2½sinðvt þ 308Þ þ 0:5 sinð3vt þ 208Þ þ 0:25 sinð5vt þ 158Þ þ 0:15 sinð7vt þ 108Þ þ 0:08 sinð9vt þ 108Þ The signal is sampled at a sampling frequency of 2000 Hz ðDT ¼ 0:5 msÞ; and a number of samples equals 200 samples are used to estimate the harmonics content of the signal using the proposed algorithm. It has been found that the proposed algorithm is succeeded in estimating the harmonics spectrum accurately. Also it has been shown, through extensive runs, that the number of samples as well as the sampling frequency has no effects on the estimated parameters of the voltage signal. Example 2. This example is solved in Ref. [9] using the fuzzy linear regression. The voltage signal waveform is contaminated with harmonics up to 19 harmonics and is given by: vðtÞ ¼ 1:0 sinðvt þ 108Þ þ 0:1 sinð3vt þ 208Þ þ 0:08 sinð5vt þ 308Þ þ 0:08 sinð9vt þ 408Þ þ 0:06 sinð11vt þ 508Þ þ 0:05 sinð13vt þ 608Þ þ 0:03 sinð19vt þ 708Þ The voltage signal is sampled at a sampling frequency of 2000 Hz and a 200-sample data window size is used. Table 1 gives the results obtained for the voltage harmonics amplitude as well as the phase angle. Examining this table reveals that: † The estimates of the proposed algorithm are accurate enough. Example 3. Actual recorded data The proposed algorithm is tested on actual recorded data, sampled at 1500 Hz, and 90 samples data window size is used. This signal is obtained by simulating a fault using EMTP. Fig. 1 gives the actual and estimated signals, while Table 2 gives the harmonic contents of the signal. Examining Fig. 1 and Table 2 reveals the following.

Fig. 1. Actual and estimated signals.

† Although the voltage signal is a highly distorted signal, the proposed algorithm is succeeded in producing good estimates for the harmonics contents of the signal. The two curves are coinciding. † During the fault period we recommend, for such a highly distorted data window, an exponential decayed dc term in the model of the voltage waveform, since the voltage signal is for a faulted phase. † We can consider harmonics of amplitude less than 3% of fundamental are not exist. 4.2. Signal with unknown frequency In this section, we assume that the signal frequency is unknown. In this case Eq. (2) will be a highly nonlinear cost function in the frequency, and harmonics parameters. Two examples are solved in this section: Example 3: The first example is the same as example 2, but the signal frequency is assumed to be unknown. We use for simulation 50 Hz, nominal frequency, one time, and another time we use 49.8 Hz, near nominal frequency. Fig. 2 shows the estimated and simulated waveforms. Examining this figure reveals that:

Table 1 Harmonics and phase angle estimate Order of harmonic

Magnitude (p.u.)

w (Degree)

1 3 5 9 11 13 19

0.99995 0.099958 0.8 0.8 0.06 0.05 0.03

10 19.976 30 40 50.0 60.0 70.00

† The proposed algorithm estimates exactly the simulated results, and gives harmonics parameters mentioned in Section 4.1 example 2. † Since the two curves are coinciding, the algorithm estimates exactly the signal frequency, nominal frequency and near nominal frequency (50 or 49.8 Hz). It has been shown through extensive runs that the proposed algorithm estimates exactly the signal frequency regardless the value of this frequency.

S.A. Soliman et al. / Electrical Power and Energy Systems 26 (2004) 31–36

35

Table 2 Estimated harmonics contents for the actual recorded data Harmonic order

Amplitude (p.u.)

Phase angle (deg.)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30

0.92937 0 0.00785 0.00017 0.00612 0 0.00121 0.00027 0 0.02976 0.0224 0.00139 0.01191 0.00022 0 0.00164 0.00003 0.02564 0.00035 0.00012 0.00055 0 0.0002 0.00462 0.00025 0.01513 0 0.04939 0.00037 0.00357

88.92156 0.07784 1.32428 0.35137 1.31303 0.02897 0.76314 0.26009 1.37055 0 12.25509 0.18803 2.8034 0.06027 0.6211 0.50068 0.26156 0 0.44847 0.06775 0.30928 0.26512 0 0.38181 0.18268 0.18616 0.49983 0.24213 0.27831 0

Example 4 While the second example, the proposed algorithm is tested on the same actual recorded data, stated in Section 4.1, but we assume that the signal frequency is unknown, and has to be estimated. The proposed algorithm is accurately estimated the frequency. This frequency has been found to be 50.15973 Hz. This value is very close to the value estimated by a frequency estimation algorithm of

Fig. 2. Estimated and simulated signals.

Fig. 3. Actual and estimated signals.

Ref. [15]. Fig. 3 compares the actual signal with the estimated signal. Table 3 gives the harmonic contents of the signal. Examining the figure and table reveals the following results: † The proposed algorithm estimates the signal with a good accuracy. Table 3 Estimated harmonics contents for the actual recorded data Harmonic order

Amplitude (p.u.)

Phase angle (deg.)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30

0.89927 0 0.00738 0 0.00628 0 0.00027 0.00067 0.0005 0.02979 0.02426 0 0.01005 0.00011 0.00074 0 0.0001 0.0214 0.00003 0 0.00148 0.00024 0.0001 0.00051 0.00029 0.00798 0.00021 0.004571 0.04999 0.00123

85.42126 1.03193 0.15929 0.39935 0.35735 0.56926 0.08843 0.23173 0.28249 0.10307 11.10537 0.24713 3.08373 0.31609 0.1531 0.48818 0 0.15004 0.29182 0.40445 1.13516 0.37999 0 0 0 0.22898 0.59784 0 0 0.39088

36

S.A. Soliman et al. / Electrical Power and Energy Systems 26 (2004) 31–36

† The harmonics contents is slightly different from those obtained in Table 2, since in this estimate, we assume the frequency of the signal is unknown, which produces a different cost function.

5. Conclusions In this article the SA based optimization algorithm is implemented to estimate harmonic contents of a signal as well as the frequency or harmonics as those available in the literature. The algorithm uses only samples of the voltage or current signal under study as well as the number of harmonics to be expected in the waveform. It has been shown that either the sampling frequency or the number of samples, providing that the sampling frequency satisfies the sampling theorem, does not affect the harmonics parameters. The proposed algorithm is suitable for off-line and at steady-state estimation.

References [1] Miller AJV, Dewe MB. The application of multi-rate digital signal processing techniques to the measurement of power system harmonic levels. IEEE Trans Power Deliv 1993;18(2):531–9. [2] Moo CS, Chang YN, Mok PP. A digital measurement scheme for time-varying transient harmonics. IEEE Trans Power Deliv 1995; 10(2):588– 94. [3] Smith BC, Waaston NR, Wood AR, Arrillaga J. A Newton solution for the harmonic phasor analysis of AC/dC. IEEE Trans Power Deliv 1996;11(2):965–71. [4] Rios S, Castaneda R, Veas D. Harmonic distortion and power factor assessment in City Street Gas discharge lamps. IEEE Trans Power Deliv 1996;11(2):1013–20. [5] Moreno Saiz VM, Barros Guadalupe J. Application of Kalman filtering for continuous real-time tracking of power system harmonics. IEE Proc Gen, Trans, Distrib (C) 1997;144(1):13 –20. [6] Probabilistic Aspects Task Force of the Harmonics Working Group Subcommittee, Time-varying harmonics: part I—characterizing measured data. IEEE Trans Power Deliv 1998;13(3):938– 44.

[7] Exposito AG, Macias JAR. Fast harmonic computation for digital relaying. IEEE Trans Power Deliv 1999;14(4):1263–8. [8] Nagndui E, Olivier G, April GE. Harmonics analysis in multipulse thyristor converters under unbalanced voltage supply using switching functions. Can J Electr Comput Engng 1999;24(4):137–47. [9] Soliman SA, Helal I, Al-Kandari AM. Fuzzy linear regression for measurement of harmonic components in a power system. Electr Power Syst Res 1999;50:99–105. [10] Bucci G, Landi C. On-line digital measurement for the quality analysis of power systems under non-sinusoidal conditions. IEEE Trans Inst Meas 1999;48(4):853–7. [11] Krim F, Benbaouche L, Chaoui A. A novel identification-based technique for harmonics estimation. Can J Electr Comput Engng 1999;24(4):149– 54. [12] Lai LL, Chan WL, Tse CT, So ATP. Real-time frequency and harmonic evaluation using artificial neural network. IEEE Trans Power Deliv 1999;14(1):52–9. [13] Lin HC, Lee CS. Enhanced FFT-based parametric algorithm for simultaneous multiple harmonics analysis. IEE Proc Gen Trans Distrib 2001;148(3):209 –14. [14] Tsao T-P, Wu R-C, Ning C-C. The optimization of spectral analysis for signal harmonics. IEEE Trans Power Deliv 2001;16(2):149–53. [15] Michie WC, Cruden A, Niewczas P, Madden WI, McDonald JR, Gauduin M. Harmonic analysis of current waveforms using optical current sensor. IEEE Instrumentation and Measurement Technology Conference, Budapest, Hungary; 2001. p. 1863–65. [16] Mantaway AH, Soliman SA, El-Hawary ME. An innovative simulated approach to the long-term hydro scheduling problem, Inter J. Elec Power and Energy Systems 2002;25(1):41–46. [17] Mantawy AH, Abdel-Magid YL, Selim SZ. A simulated annealing algorithm for unit commitment. IEEE Trans Power Syst 1997;13(1): 197 –204. [18] Aarts E, Korst J. Simulated annealing and Boltzman machines. A stochastic approach to combinatorial optimization and neural computing, New York: Wiley; 1989. [19] Cerny V. Thermo dynamical approach to the travelling salesman problem: an efficient simulation algorithm. J Optimization Theory Appl 1985;45(1). [20] Selim SZ, Alsultan K. A simulated annealing algorithm for the clustering problem. Pattern Recognit 1991;24(10):1003–8. [21] Metropolis N, Rosenbluth A, Rosenbluth M, Teller A, Teller E. Equations of state calculations by fast computing machines. J Chem Phys 1953;21:1087 –92. [22] Siarry P, Berthiau G, Haussy J. Enhanced simulated annealing for globally minimizing functions of many continuous variables. ACM Trans Math Software 1997;23(2):209–22.