EP Based Optimization for Estimating Synchronizing and Damping Torque Coefficients

This paper presents Evolutionary Programming (EP) based optimization technique for estimating synchronizing torque coefficients, Ks and damping torque coefficients, Kd of a synchronous machine. These coefficients are used to identify the angle stability of a system. Initially, a Simulink model was utilized to generate the time domain response of rotor angle under various loading conditions. EP was then implemented to optimize the values of Ks and Kd within the same loading conditions. Result obtained from the experiment are very promising and revealed that it outperformed the Least Square (LS) method and Artificial Immune System (AIS) during the comparative studies. Validation with respect to eigenvalues determination confirmed that the proposed technique is feasible to solve the angle stability problems.


INTRODUCTION
Small signal stability analysis of power systems becomes more important nowadays. Under small perturbations, this analysis is to predict the low frequency electromechanical oscillations resulting from poorly damped rotor oscillations. These oscillations stability becomes a very important issue as reported in [Feilat et al. (1999), Rouco et al. (1993), Abu-AlFeilat (2000), El Naggar et al. (2002), El Naggar et al. (2004]. The operating conditions of the power system are change with time due to the dynamic nature, so it is need to track the system stability on-line. To track the system, some stability indicators will be estimated from given data and these indicators will be updated as new data received. Synchronizing torque coefficients K s and damping torque coefficients K d are used as stability indicators. To achieve stable operation of the machine, both K s and K d must be positive [Peter et al. (1998), Glover et al. (2002), Kundur ( 2004), Hadi (2004)].
Certain techniques have been proposed to estimate the value of K s and K d which involved optimization technique. Some techniques have been explored by means of frequency response analysis [Shepperd (1961), DeMello et al. (1968), deOliveria et al., (1994, Padiyar et al. (1991)]. [Glover et al. (2002)] decomposes the change in electromagnetic torque into two orthogonal components in the frequency domain. The two equations are expressed in terms of the load angle deviation then solved directly. Static and dynamic time domain estimation methods were also proposed in this study.
Least Square (LS) method can be one of the possible techniques in addressing this phenomenon. It has been used as static parameter estimation [Alden et al. (1979)]. However, several disadvantages have been identified with LS. Amongst them are the long computation time and the requirement of data updating. It also requires monitoring the entire period of oscillation. Recently, evolutionary algorithms such as Evolutionary Programming (EP) and Artificial Intelligent System (AIS) have received much attention for global optimization problems. These evolutionary algorithms are heuristic population-based search methods that used both random variation and selection. The search for an optimal solution is based on the natural process of biological evolution and is accomplished in a parallel method in the parameter search space. EP-based method has been applied in various researches in static and dynamic system stability [Dobson et al. (1992), Lai et al. (1998), Jason et al. (1999), Fogel et al. (2000, Abido et al. (2002), Musirin et al. (2002), Rajan (2004), Hassim et al. (2006, Talib et al. (2007)]. On the other hand, Artificial Immune System (AIS) approach to optimization is more recent exploitation of natural phenomena in power system than EP. EP and AIS share many common aspects; EP tries to model the natural evolution while AIS tries to benefit from the characteristics of a human immune system [Rahman et al. (2004), Dasgupta (2006), Liu et al. (2006, Hunjan et al. (2007), Wei et al. (2008]. This paper presents an efficient online estimation technique of synchronizing and damping torque coefficients in solving angle stability problems. The method is based upon the population-based search methods that use both random variation and selection. The method is used to estimate synchronizing torque coefficients, K s and damping torque coefficients, K d from the machine time responses of the change in rotor angle ǻį(t), the change in rotor speed ǻȦ(t) and the change in electromechanical torque ǻT e (t). The goal is to minimize the estimated coefficient error and the time consumed. The proposed EP technique is used to find the best solution of the formulated problem. Results obtained from the experiment using EP have been compared with AIS and LS; resulting significant merit.

The System Model:
A simplified block diagram model of the small-signal performance is shown in Fig. 1. In this representation, the dynamic characteristics of the system are expressed in terms of K constants with linearized single machine infinite bus (SMIB) system, this model is represented with some variables such as electrical torque, rotor speed, rotor angle and exciter output voltage. From the transfer function block diagram the following state-space form is developed. The system matrix A is a function of the system parameters, which depends on the opening conditions. The perturbation matrix B depends on the system parameters only. Details of matrix A and matrix B are explained in Appendix.
The interaction among these variables is expressed in terms of the 4 constants K 1 , K 2 , K 3 and K 4 . These constants with the exception of K 3 , which are only a function of the ratio of impedance, are function of the operating real and reactive loading as well as the excitation levels in the generator.

Concept of Synchronizing and Damping Torque:
A single machine connected to infinite bus system is considered. The system comprises a steam generator connected via a tie line to a large system represented as infinite bus. The dynamic stability study is performed by linearizing the power system under consideration around an operating point to represent the system is state space model. The machine differential equations, the exciter equation and the block diagram can be found in [Kundur (1994)].
The change of electromagnetic torque ǻT e (t) can be broken down into two components namely the synchronizing torque K s and damping torque K d . The synchronizing torque component is in phase and proportional with the change in rotor angle ǻį(t), and the damping torque is in phase and proportional with the change in rotor speed ǻȦ(t). The estimated torque can be written as: Where: ǻį(t) : Change in rotor angle ǻȦ(t) : Change in rotor speed K s : Synchronizing torque coefficients K d : Damping torque coefficients Evolutionary Programming: The Evolutionary Programming (EP) is one of the evolutionary computing which uses the models of biological evolutionary process for the solution of complex engineering problems. The search for an optimal solution is based on the natural process of biological evolution and is accomplished in a parallel method in the parameter search space. EP belongs to the generic fields of the simulated evolution and artificial life. It is robust, flexible and adaptable and it can yield global solutions to any problem, whatever the form of the objective function.
The advantages of EP over other conventional optimization techniques can be summarized as follows [Dobson et al. (1992), Lai et al. (1998), Jason et al. (1999), Fogel et al. (2000, Abido et al. (2002), Musirin et al. (2002), Rajan (2004), Hassim et al. (2006, Talib et al. (2007)]: 1. EP searches the problem space using a population of trials representing possible solutions to the problem and not a single point. This will ensure that EP is less possibility getting trapped on local minima. Therefore, EP can reach to a global optimal solution. 2. EP uses performance index or objective function information to guide the search for solution. Therefore, EP can easily deal with non-smooth and non-continuous objective functions. 3. EP uses probabilistic transition rules instead of non-deterministic rules to make decisions. Moreover, EP is a kind of stochastic optimization algorithm that can search a complicated and uncertain area to find the global minimum. This makes EP more flexible and robust than conventional methods.

Algorithm for EP:
In the EP algorithm, the population has 2n candidate solutions with each candidate solution is an mdimensional vector, where m is the number of optimized parameters. The EP algorithm can be described as: Step 1 (Initialization): Generation counter i is set to 0, and generate n random solutions . In this initial population, minimum value of objective function J min will be searched, the target is to find the best solution x best with objective function J best Step 2 (Mutation): Each parent x k produces one offspring x k+n . Each optimized parameter P l is perturbed by a Gaussian random variable . The standard deviation specifies the range of the optimized 2 0, l N V l V parameter perturbation in the offspring. equation is as follows: where ȕ is a scaling factor, and is the objective function of the trial solution x k .

k J x
The value of optimized parameter will be set at certain limit if any value violates its specified range. The offspring can be described as: Step 5 (Combination): All members in the population are combined with all members from the offspring k x to become 4n candidates. These individuals are then ranked in descending order, based on their k n x objective function as their weight.
Step 6 (Selection): The first n individuals with higher weights are selected are selected along with their objective functions as parents of the next generation. The generation counter will be set to and 1 i i algorithm will start again from Step 2.
Step 7 (Stopping criterion): The search process will be terminated if one of the followings is satisfied: 1. It reaches the maximum number of generations 2. The value of is very close to 0. max min

J J
The flow chart of EP is shown in Fig. 2.

Application of EP:
The EP algorithm described before has been applied to search for optimal or near optimal values of K s and K d . In our implementation, the search will terminate if the following occur: It reaches the 1000 number of generations The value of . max min 0.0001 J J t

Least Square Method:
All the data of ǻį(t), ǻȦ(t) and ǻT e (t) can be obtained from either offline simulation or online measurements. Following a small disturbance, the time responses of these three items are recorded. The least square (LS) technique is then used to minimize the sum of the square of the differences between the electric torque ǻT e (t) and the estimated torque . The error is defined as: e e

E t T t T t ' '
The torque coefficients K s and K d are calculated to minimize the sum of the error squared over the interval of oscillation t, as given in Eq.(4), where t=NT (N is the number of samples and T is the sampling period). For correct estimation of K s and K d , the interval t should be chosen adequately. The suitable value of t which makes K s and K d constants during the oscillation period was found to be the entire period of oscillation. In matrix notation, the above problem can be described by an over-determined system of linear equations as follows: where , and . The estimated vector x is such that the function In this case the estimated vector x will be given by: where is the left pseudo inverse matrix. Solving Eq. (7) gives the values of K s and K s for the t A corresponding operating point.

Artificial Immune System:
Artificial immune system (AIS) approach to optimization is more recent exploitation of natural phenomena in power system than EP. EP and AIS share many common aspects, whereas EP tries to model the natural evolution, AIS try to benefit from the characteristics of a human immune system. Basic algorithm for AISbased optimization is called the Clonal Selection Algorithm (CSA) and it works as follows [Rahman et al. N affinity values. The better the affinity the higher the number of clones generated for each of the selected antibodies. The clones form a group C. 5. The clone group C undergoes an affinity maturation, in which the clones are mutated inversely proportionally to their affinities: the better the affinity the smaller the mutation rate. 6. The affinities of the affinity-matured clones are calculated. 7. If an affinity-matured clone has a better affinity value than the parent antibody, replace the parent antibody with the affinity-matured clone. 8. Replace the d lowest affinity antibodies with randomly created new antibodies. 9. Go to c) if run time constraints have not been met, otherwise exit.
A flow chart of AIS is shown in Fig. 3.

Test System:
In this study, performance evaluation of the EP for the estimation of K s and K d is compared with LS and AIS estimation method. The evaluation is carried out by conducting several offline simulation cases on the Linearized model of SIMB. In this study, block diagram as shown in Fig. 1 is used for offline simulation to generate the required ǻį(t), ǻȦ(t) and samples in MATLAB Simulink environment. The parameters e T t ' of the SMIB system are given in Appendix. Stable and unstable study cases are simulated using different types of disturbances. Data size is set to 20 second, while number of samples is set to 400 samples. Using ǻį(t), ǻȦ(t) and as generated sample e T t ' data, 3 sets of MATLAB files: LS, EP and AIS based simulation are developed. The simulation diagram is show in Fig. 4. During simulation, all parameters are adjusted until an optimal solution is obtained. The results of EP and AIS are compared with the LS solution.  In this study, 8 sets of ǻį(t), ǻȦ(t) and ǻT e (t) samples are generated using offline simulation of block diagram implemented in MATLAB Simulink. 8 samples of ǻį(t) data in graph forms are shown in Fig. 5 until Fig. 12. Different value of P and Q are used to simulate these cases. For verification, the eigenvalues for all cases have been calculated and written below of each graph. For cases which all eigenvalues are negative, it is stable. For cases which have positive eigenvalues, it is unstable. From the result of eigenvalues, the first 5 cases (Fig. 5~9) are stable and the other 3 cases (Fig. 10~12) are unstable. Table 1 shows the results obtained for eight different study cases. The estimated constants obtained using EP, AIS and LS methods are shown as well as the eigenvalues for each case. Results shown in this table for the proposed method are the steady state values obtained at the end of the simulation. It is found that 400 samples within a data size of 20 seconds after the disturbance are sufficient to reach steady state solution.
As both values of K s and K d are positive, the result indicate that case 1, 2, 3, 4 and 5 are stable cases. On the other hand, case 6, 7 and 8 are unstable cases as the value of K d is negative. Eigenvalues shown in the last column verify the result obtained.
In all cases, all 3 methods give accurate and close results. Although results using EP and AIS method are close, difference of value between EP method and LS method is closer than different of value between AIS method and LS method. These show that simulation results from EP method are more accurate and closer than simulation results using AIS method.
Except for LS method, times consume to calculate the value of K s and K d until it reaches steady state solution has been recorded for EP and AIS method. As the value of K s and K d calculated is same from first iteration, the time consume for LS method is not been recorded. Comparing EP and AIS for all cases, the average of time consume to calculate using AIS is about 48 seconds, while EP is about 30 seconds, the time almost 40% longer compared to EP. As a result it shows that calculation method using EP is faster than AIS.
For the effect of error-contaminated data on the accuracy of the estimated value of K s and K d , simulation has been done by introducing about 10% of bad data (zeros) at different locations of ǻį(t), ǻȦ(t) and ǻT e (t). For comparison, case 5 and 8 has been selected.  Table 2 shows the result of EP, AIS and LS method calculation with bad data. It shows that although bad data were injected into the system, estimates for K s and K d using EP and AIS method are not affected as the value is identical with the results obtained in Table 1. On the other hand, LS is affected with the 10% of bad data which has been implemented.     More than that, it also gives false result for case 8. Using LS method, the simulation gives positive value of damping coefficient K d that indicates the system is stable for case 8, but the results given by EP and AIS method give negative value which indicates the system for case 8 is unstable. The result also can be confirmed by the positive value of eigenvalues which verify that case 8 is unstable. As a result, estimates for K s and K d using EP and AIS method are more accurate compare to LS method. Fig. 13 shows convergence process for case 8 without bad data, while Fig. 14 shows convergence process for case 8 with bad data.

Conclusion:
In this paper, three methods for accurate estimation of the synchronizing and damping torque coefficients, K s and K d are presented. The performance of Evolutionary Programming (EP) is compared with the Artificial Immune System (AIS) and Least Square (LS) method. Compared with AIS and LS, EP gives several advantages. This includes better data accuracy and 60% shorter computing time compared to AIS. EP also never affected with bad data consumed in the system compared to LS which give false decision on the stability. The proposed method can be considered as a reliable and efficient tool in the area of power system stability analysis.