Box-Muller harmony search algorithm for optimal coordination of directional overcurrent relays in power system

Directional overcurrent relays are used to protect interconnected power systems and looped distribution systems. The problem is modeled as a constrained nonlinear optimization problem in which the decision variables are the devices that control the act of isolation of faulty lines from the system without disturbing the healthy lines. Time dial setting (TDS) and plug setting (PS) are considered as optimal parameters. IEEE 3-bus model and IEEE 4-bus model are two models of this problem. This paper presents a new efficient and reliable approach based on a constrained harmony search algorithm (HSA) to solve the problems. The authors’ proposed technique is called Box-Muller harmony search (BMHS) algorithm. BMHS has been tested by applying it to two models. The simulation result of this technique is compared with those reported in the literature. The outcome is very encouraging and proves that the new technique outperforms them in terms of reaching a more optimal solution and speed.


INTRODUCTION
The problem of determining the time dial setting of directional overcurrent relays using the optimization techniques was first stated in 1987 (Perez and Urdaneta, 1999;Urdaneta et al., 1998).It has been shown that the linear programming technique can be successfully applied to this problem, guaranteeing the minimum possible settings of the relays that satisfy the time coordination constraints, that is, the optimal settings (Perez and Urdaneta, 1999;Urdaneta et al., 1998;Urdaneta et al., 1996;Chattopadhay et al. 1996;Urdaneta et al., 1997).The problem of re-setting the relays after a permanent topology change in the system reducing the number of relays to be reset was presented and solved using the concepts of multiple objective optimizations (Chattopadhay et al., 1996).The consideration of the dynamic changes of the system topology that take place during the fault clearing process *Corresponding author.E-mail: pp_shafipour@yahoo.com.using linear programming was also recently addressed (Urdaneta et al., 1997).However, the optimization problem has been formulated and solved for the systems including only directional overcurrent relays (DOCR) without considering the coordination of DOCR with other relay types.

OVERVIEW
As the dimension of the problem increases for the modern interconnected power systems, the complexity of the problem increases.Also, due to the complexities of non-linear programming techniques, most of the researchers have solved the problem in a linear environment by assuming the values of decision variables (all plug settings), which make the problem nonlinear.This assumption is made on the basis of engineering experience (Chattopadhay et al., 1996), (Irving and Elrafie, 1993;Urdaneta et al., 1996;Urdaneta et al., 2001), etc.The use of optimization techniques in relay coordination was first suggested by Urdaneta et al. (1998).Irving and Elrafie (1993) used the Sparse Dual Revised Simplex method of linear programming suggested by Irving and Sterling (1983) to optimize TDS settings for assumed non-linear PS settings.Laway and Gupta (1993) applied Simplex and Rosenbrock-Hillclimb methods to optimize TDS and PS settings respectively; in a similar way, as used by Urdaneta et al. (1998).These approaches were further followed by the simplex-based approaches with more and more sophistications about finer aspects of the relays (Urdaneta et al., 1996), Chattopadhyay et al. (1996), Perez and Urdaneta (1999), Abdelaziz et al. (2002) and So and Li (2000) used evolutionary programming.A survey of all coordination philosophies used by the various researchers in the past has been presented recently by Birla et al. (2005).Deep et al. (2006) used random search technique 2(RST2) of Shanker and Mohan (1987) to solve the relay coordination problem for IEEE 3-bus and IEEE 4-bus models.Dipti (2007), applied genetic algorithm (GA), self organizing migrating algorithm (SOMA), self organizing migrating genetic algorithm (SOMGA) which is a genetic algorithm hybridized with self organizing migrating algorithm to solve the problem.Bansal and Deep (2008) used the five different versions of Particle Swarm Optimization (PSOG, PSOGC, PSOL, PSOLC, CPSO) to solve the relay coordination problem.Barzegari et al. (2010) applied Harmony Search Algorithm (HSA) to solve the problem.

PROBLEM FORMULATION
The operating time (T) of a DOCR is the non-linear function of the relay settings (time dial setting (TDS) and plug setting (PS) and the fault current (I) seen by the relay).Therefore, relay operating-time equation for a directional overcurrent relay is given by a non-linear equation as follows: (1) *denotes the multiplication.Only TDS and PS are the unknown variables in the equation.These are the "decision variables" of the problem, and are the constants representing the behavior of characteristic in a mathematical way, in which operating time of the DOCR varies and are given as 0.14, 0.02 and 1.0 respectively as per IEEE standard (1997).Value of I is also known, as it is a system dependent parameter and continuously measured by measuring instruments (Bansal and Deep, 2008).
The * _ variable is equal to .Characteristics of standards of overcurrent relay are shown in Table 1 (Kavehnia et al., 2006).
The relay, which is supposed to operate first to clear the fault, is called the primary relay.A fault close to relay is known as the close-in fault for the relay and a fault at the other end of the line is known as a far-bus fault for this relay.Conventionally, objective function in coordination studies is constituted as the summation of operating-times of all primary relays, responding to clear all close-in and far-bus faults.These times have limits.
As what is expressed in Abyaneh et al. (2003) in order to coordinate two overcurrent relays, one as the main relay (m) and the other as the backup relay (b), (it is shown in Figure 1), the difference between the operation time of the backup relay and the operation time of the main relay for faults and should be more than CTI and are short circuits at near bus and far bus of the main relay respectively.CTI is the coordination time interval for main and backup relays.
So the constraints for coordination of relays m and b will be in the form of Inequality (2).
(2) is operating time of the backup relay and is operating time of the primary relay.Value of CTI is known.
So that the formulation problem is as follows (Bansal and Deep, 2008): Minimize Subject to: (3) , is number of relays responding for close-in fault, is number of relays responding for far-bus fault, − − is primary relay operating-time for close-in fault.Here: (4) − − is primary relay operating-time for far-bus fault.Here: (5) The values of constants i a , i b , i c and i d are given in the two Tables , 2 and 5 (Bansal and Deep, 2008).
is operating time of the backup relay.Here: is operating time of primary relay.
Here: The values of constants i e , i f , j g and j h are given in Tables 3 and 6 (Bansal and Deep, 2008).
According to Table 1 and Equation (1

The harmony search algorithm
The harmony search algorithm (HSA) is a music-inspired evolutionary algorithm, mimicking the improvisation process of music players (Geem et al., 2001;Geem et al., 2008).The HSA is simple in concept, few in parameters, and easy in implementation, with a theoretical background of stochastic derivative (Geem et al., 2008).The algorithm was originally developed for the discrete optimization and later expanded for the continuous optimization (Lee and Geem, 2005).For continuous variables, Das et al. (2011) provided theoretical background of exploratory power.It has been successfully applied to various benchmark and real-world problems including travelling salesman problem (Geem et al., 2005), parameter optimization of river flood model (Kim et al., 2001), design of pipeline network (Geem et al., 2002;Geem, 2006), and design of truss structures (Lee and Geem, 2004).Recently, parameter-setting-free (AKA adaptive) HS was developed (Geem and Sim, 2010).
The steps in the procedure of harmony search are as follows (Fesanghary et al., 2008): Step 1: Initialize the problem and algorithm parameters Step 2: Initialize the harmony memory Step 3: Improvise a new harmony Step 4: Update the harmony memory Step 5: Check the stopping criterion These steps are described as follows:

Initialize the problem and algorithm parameters
In Step 1, the optimization problem is specified as follows: Minimize r is the objective function, M is the number of inequality constraints and P is the number of equality constraints.x is the set of each decision variable i x .N is the number of decision variables.
The lower and upper bounds for each decision variable are i L x and i U x respectively.The HSA parameters are also specified in this step.These are the harmony memory size (HMS), or the number of solution vectors in the harmony memory, harmony memory considering rate (HMCR), pitch adjusting rate (PAR), and the number of improvisations (NI), or stopping criterion.The harmony memory (HM) is a memory location where all the solution vectors (sets of decision variables) are stored.The HM is similar to the genetic pool in the genetic algorithms (GAs) (Fesanghary et al., 2008).
Here, HMCR and PAR are parameters that are used to improve the solution vector.Both are defined in Step 3.

Initialize the harmony memory
In Step 2, the HM matrix is generated from a uniform distribution in the ranges ] , [ Infeasible solutions that violate the constraints have a chance to be included in the HM with the hope of forcing the search towards the feasible solution area.Static penalty functions are used to calculate the penalty cost for an infeasible solution.The total cost for each solution vector is evaluated using: Where i α and j β are the penalty coefficients.Generally, it is difficult to find a specific rule to determine the values of the penalty coefficients and normally these parameters remain problemdependent.Selecting a suitable value of penalty coefficients, however, proves to be very difficult.If they are chosen to be too large, the search terrain of ) x fitness( r may become too rugged to be searched by gradient-based methods.If they are too small, the solution to unconstrained problem Equation (11) may not be a CLM (constrained local minima) to the original problem or even may not be a feasible solution.Usually the selection of penalty coefficients is randomly done (Wang, 2001;Wang et al., 2005).

Improvise a new harmony
, is generated based on three rules: (1) Memory consideration (2) Pitch adjustment (3) Random selection Generating a new harmony is called "improvisation" (Fesanghary et al., 2008).In the memory consideration, the value of the first decision variable x ′ for the new vector is chosen from any of the values in the specified HM range are chosen in the same Fetanat et al. 4083 manner.The HMCR, which varies between 0 and 1, is the rate of choosing one value from the historical values stored in the HM, while (1-HMCR) is the rate of randomly selecting one value from the possible range of values, as shown in Equation ( 12). end.
Where rand () is a uniform random number between 0 and 1 and i X is the set of the possible range of values for each decision variable, that is For example, a HMCR of 0.85 indicates that the HSA will choose the decision variable value from historically stored values in the HM with an 85% probability or from the entire possible range with a (100 to 85)% probability (Fesanghary et al., 2008).Every component obtained by the memory consideration is examined to determine whether it should be pitch adjusted.This operation uses the PAR parameter, which is the rate of pitch adjustment as follows: . end Where bw is an arbitrary distance bandwidth.To improve the performance of the HSA and eliminate the drawbacks associated with fixed values of PAR and bw, Mahdavi et al. (2008) proposed an improved harmony search (IHS) algorithm that uses variable PAR and bw in improvisation step.In their method PAR and bw change dynamically with generation number as expressed below: Recently other variants of harmony search have been proposed.Omran and Mahdavi (2008) proposed a new variant of harmony search, called the global best harmony search (GHS), in which concepts from swarm intelligence are borrowed to enhance the performance of HSA such that the new harmony can mimic the best harmony in the HM.Geem (2009) claimed that HS is better than GHS for a large-scale problem with 454 discrete variables while GHS is better than HS for small or medium sized problems.Also, Geem (2008) proposed a new stochastic derivative for the discrete variables based on a harmony search algorithm to optimize the problems with discrete variables and the problems in which the mathematical derivative of the function cannot be analytically obtained.

Update harmony memory
If the new harmony vector, ) , , , ( , has a better fitness function than the worst harmony in the HM, the new harmony is included in the HM and the existing worst harmony is excluded from the HM.

Check stopping criterion
The HSA is terminated when the stopping criterion (example, maximum number of improvisations) has been met.Otherwise, Steps 3 and 4 are repeated.

Box-Muller harmony search algorithm (BMHS)
Inspired by the concept of Box-Muller transform, a new variation of HSA is proposed in this study.A Box-Muller transform is a method of generating pairs of independent standard normally distributed (zero expectation, unit variance) random numbers, given a source of uniformly distributed random numbers (Box and Muller, 1958).
It is commonly expressed in two forms.The basic form as given by Box and Muller takes two samples from the uniform distribution on the interval (0, 1) and maps them to two normally distributed samples.The polar form takes two samples from a different interval, (−1, +1), and maps them to two normally distributed samples without the use of sine or cosine functions.One could use the inverse transform sampling method to generate normally-distributed random numbers instead; the Box-Muller transform was developed to be more computationally efficient.Suppose 1 U and 2 U are independent random variables that are uniformly distributed in the interval (0, 1): Z and 2 Z are independent random variables with a normal distribution of standard deviation 1.
The derivation is based on the fact that, in a two-dimensional Cartesian system where X and Y coordinates are described by two independent and normally distributed random variables, the random variables for 2 R and θ (shown earlier) in the corresponding polar coordinates are also independent and can be expressed as: The Box-Muller harmony search (BMHS) algorithm generates random numbers by means of Box-Muller method.BMHS has exactly the same steps as the HSA but the only difference lies in the Steps 2 and 3 which are modified as presented in Figure 2. In Figure 2, Z-rand generates random number in range 0.0~1.0by.method.

RESULTS
To examine the applicability and accuracy of the proposed method, BMHS has been applied to two models.The obtained results are compared with the results of the other methods.Programming commands used in this simulation are taken from VB.Net software.

3-bus system
For the coordination problem of IEEE 3-bus model, value of each of and is 6 (equal to number of relays or twice the lines).Accordingly, there are 12 decision variables (two for each relay) in this problem, that is, 1 to 6 and 1 to 6 .The 3-bus system can be visualized as shown in Figure 3, also Table 2 shows i a , i b , i c and i d and Table 3 shows operating time of backup relay and operating time of primary relay.
The problem formulated in problem formulation is solved using BMHS technique.The parameters used for the simulation are as follows: HMCR = 0.85, min PAR = 0.05, max PAR = 0.95, HMS (The population size) = 10 and NI=10000.Also, In Equation( 11), i α and j β respectively are 80 and 1000.Figure 4 illustrates the convergence of BMHS where it clearly shows that the optimal solution is obtained after 9500 iterations.
According to Equation (1), if PS is constant then the time based on TSM changes linearly.Thus, objective function and constraints become linear based on TSMs.So the optimization problem becomes a linear programming (LP) form (Perez and Urdaneta, 1999;Zeineldin et al., 2005;Yue et al., 2006;Braga and Saraiva, 1996).In this model employed in this paper, obtained PSs via BMHS were put into objective function and constraints.Then the problem, using general algebraic modeling system (GAMS) software and MATLAB Optimization Toolbox, was solved.
Table 4, shows a comparison of the optimal values of decision variables (DV) and objective function (OBJ) obtained by BMHS, LP and other algorithms for IEEE 3bus model.

4-bus system
The next coordination problem is of IEEE 4-bus model, value of each of cl N and far N is 8 (equal to number of relays or twice the lines).Accordingly, there are 16 decision variables (two for each relay) in this problem, that is, 1 to 8 and 1 to 8 .The 4-bus system can be visualized as shown in Figure 5.The objective   ).The values of constants i a , i b , i c , i d and i e , i f , j g , j h for IEEE 4-bus model are given in Tables 5  and 6 respectively.β respectively are 80 and 1000.Figure 6 shows the convergence of BMHS.This solution is  obtained after 19000 iterations.It is important to mention that this solution was obtained using HSA after 100,000 iterations.Again, in the same manner as previous section, the LP problem is solved with GAMS and MATLAB Optimization Toolbox.Table 7 shows the results for IEEE 3-bus model obtained by BMHS, LP and other algorithms.

DISCUSSION
Depending on the results presented on Tables 4 and 7,   for the two models, we find that there is big change on problem parameters before and after applying BMHS.
Using the two random functions in Box-Muller transform, we can make more changes in solutions; hence, the efficiency and flexibility of the problem increase.Therefore, Box-Muller transform, when applying the proposed approach on the optimization problem, affects the optimization parameters values.From Figures 4 and  6, it is clear that the BMHS method is fast and reaches the optimal solution with a relatively low number of iterations.This behavior represents a powerful feature of this new method.In this study, the proposed method shows better performance than the other algorithms.

Conclusion
The problem is to determine the optimal values of time dial setting (TDS) and plug setting (PS) so that the sum of the relays operation time can be minimized.Up to this time, linear programming (LP) and other methods are applied in this problem; each of the aforementioned methods has advantages and disadvantages.The harmony search algorithm (HSA) is one of the metaheuristic methods that has a high degree of flexibility, fast convergence, and implement ability on any multi-bus power systems.In this paper, using Box-Muller transform previous harmony algorithms speed has increased that authors called it BMHS here.Two models of this problem namely IEEE 3-bus and IEEE 4-bus are solved using BMHS.Two models have been tested and the obtained results using BMHS in Tables 4 and 7 have been compared with other optimization methods; the result indicates its superiority.Based on this study, BMHS can be applied satisfactorily to different types of optimization problems in the area of power systems.
gn) is the pitch adjusting rate for each generation, min PAR is the minimum pitch adjusting rate, max PAR is the maximum pitch adjusting rate and gn is the generation number.

Figure 4 .
Figure 4.The convergence of BMHS for IEEE 3-bus model.

Table 1 .
Characteristics of standards of overcurrent relay.

Table 2 .
Values of constants i

Table 3 .
Values of constants i

Table 4 .
Optimal decision variables and optimal values of objective function for IEEE 3-bus model.

Table 5 .
Values of constants i Figure 6.The convergence of BMHS for IEEE 4-bus model.

Table 6 .
Values of constants i e , i f , j g and j h for IEEE 4-bus model.

Table 7 .
Optimal decision variables and objective function optimal values for IEEE 4-bus model.