Dynamic economic dispatch using hybrid metaheuristics

Dynamic economic dispatch problem or DED is an extension of static economic dispatch problem or SED which is used to determine the generation schedule of the committed units so as to meet the predicted load demand over a time horizon at minimum operating cost under ramp rate constraints and other constraints. This work presents an efficient hybrid method based on particle swarm optimization (PSO) and termite colony optimization (TCO) for solving DED problem. The hybrid method employs PSO for global search and TCO for local search in an interleaved mode towards finding the optimal solution. After the first round iteration of local search by TCO, the best local solutions are considered by PSO to update the schedules globally. In the next round, TCO performs local search around each solution found by PSO. This paper reports the methodology and result of application of PSO–TCO hybrid to 5-unit, 10-unit and 30-unit power dispatch problems; the result shows that the PSO–TCO (HPSTCO) gives improved solution compared to PSO or TCO (when applied separately) and also other hybrid methods.


Introduction
Usually the economic load dispatch problem (ELD) implies static economic dispatch problem or SED where the objective is to determine the optimal schedule of online generating units' outputs so as to meet the load demand at a certain time at the minimum operating cost under various system and operational constraints. In contrast, the objective of the dynamic economic dispatch (DED) problem is to schedule the generator outputs over a certain period of time economically. The DED problem takes into consideration the limits on the generator ramping rate coupled with real time intervals to keep the thermal stress on the generation equipment like the turbines and boilers within the safe limits and thus protect their life [1]. The DED problem divides the dispatch period into a number of small time intervals, and a SED is employed to solve the problem in each interval.
Since the DED problem was introduced in 1980s, several optimization techniques and procedures have been used for solving the DED problem with complex objective functions or constraints. There were a number of classical methods that have been applied to solve this problem such as the lambda iterative method [2], gradient projection method  [3], Lagrange relaxation [4], linear programming [5], dynamic programming [6] and interior point method [7]. Most of these methods are not applicable for non-smooth or non-convex cost functions. To overcome this problem, many heuristic optimization methods have been employed to solve the DED problem; such methods include ant colony optimization (ACO) [8], particle swarm optimization (PSO) [9][10][11], Levenberg-Marquardt back-propagation algorithm (LMBPA) [12], differential evolution (DE) [13], artificial immune system (AIS) algorithm [14], harmony search (HS) [15] and bee swarm optimization (BSO) algorithm [16] among others. Many of these techniques have proved their effectiveness in solving the DED problem without any or fewer restrictions on the shape of the cost function curves. These approaches solve the DED by employing an initial population of individuals each of which represents a candidate solution for the problem. Then, they evolve the initial population by successively applying a set of operators on the old solutions to transform them into new solutions. In recent years, the trend of solving DED problems has changed from single-heuristic techniques to hybrid metaheuristics-a combination of two or more techniques like PSO-ACO, DE-SQP, PSO-SQP, etc. It is proved that these hybrid techniques have capability to solve the DED problems better than the single-heuristic problems as the hybridization causes the individual techniques to mitigate their limitations and complement each other with their characteristic strength.

Earlier hybrids
In 2005, Victoire et al. proposed hybrid EP-SQP made up of evolutionary programming and sequential quadratic programming technique to solve DED problems [17]. They also experimented with a hybrid of PSO and SQP techniques to solve DED problem with valve point loading (VPL) effect [18]. In 2009, Yuan et al. [19] hybridized PSO with differential evolution method for solving DED with VPL. In 2010, hybrid SOA-SQP method that combined seeker optimization algorithm (SOA) with SQP was used by Sivasubramani and Swarup [20] to solve DED with VPL. In 2012, Cai et al. [21] reported application of hybrid CPSO-SQP in DED with VPL. Again in 2012, Swain et al. [22] hybridized gravitational search algorithm (GSA) with SQP (GSA-SQP) to solve DED with VPL.
Elaiw et al. [23] compared, in 2013, the efficacy of hybrids DE-SQP and PSO-SQP in solving DED with VPL effect. Chen et al. [24] used a combination of three methods, namely low-discrepancy sequences (LDS), improved shuffled frog leaping algorithm (ISFLA) and SQP, to solve DED problem. In this hybrid (termed as LDISS), LDS is used to generate initial population, ISFLA is liable for global search and SQP is used for local search. In 2013, Mohammadi-Ivatloo et al. [25] introduced hybrid immune genetic algorithm to solve DED considering VPL and prohibited operating zone and ramp rate constraints along with transmission losses. Zhang et al. [26] proposed hybrid bare bones (BB)-PSO or BBPSO in 2014 to solve DED with VPL only.
Among recent developments are two hybrid techniques, BBO-PSOTVAC and FA-PSOTVAC, developed in 2018 by Hamed et al. [27], combining firefly algorithm (FA) and biogeography-based optimization (BBO) with time-varying acceleration-based particle swarm optimization (PSOTVAC) to improve the solution of DED. In the same year, Pan et al. [28] solved the DED problem with VPL using a hybrid technique MILP-IPM involving mixed-integer linear programming (MILP) and interior point method (IPM). Another very recent development (in 2018) by Xiong and Shi [29] is a hybrid of BBO and brain storm optimization (BSO) to get a better solution for DED with VPL.
In the present work, the authors have applied for the first time a hybrid computational approach HPSTCO that combines PSO and TCO to solve the DED problem. The hybrid method employs PSO iterations for global search and TCO iterations for exploring the locality near the global solutions, interleaving both the search processes to overcome the drawback of fast convergence to (selection of ) global optimal solution in the original PSO method. The interleaving process requires PSO and TCO pass their solutions to each other. The solutions of TCO are updated in PSO iterations by considering the global best solution. Similarly, the solutions of PSO are adjusted by considering the locally observed information by TCO. A solution point searched by the PSO method can be used as an initial condition in the TCO method. The hybrid HPSTCO model has been programmed in MATLAB and simulation run executed for 5-unit, 10-unit and 30-unit DED system with parameters referred from the literature. Performance of the HPSTCO, as compared using a benchmark function, is quite encouraging.

Motivation and contribution
The novelty of the present study lies in the fact that PSO and TCO together, i.e., their hybrid combination (HPSTCO) has never been tried before to optimize small to largescale economic load dispatch problem. However, PSO and TCO individually and in combination with other metaheuristics have been applied earlier in different ELD problems. The contribution of the paper lies in the fact that it has established by reporting four distinct test cases and comparing the result with other hybrid methods for each of these four test cases that the HPSTCO hybrid is quite effective and advantageous in dealing with small-scale (5-and 10-unit) as well as medium-scale (30-unit) DED problem. In medium to large-scale systems with higher-capacity turbines, the fuel cost function is highly non-smooth and non-convex and contains discontinuous values at each boundary, forming multiple local optima. The complexity of the problem also increases significantly with the increase in the number of generating units because of their combinatorial nature. The present work has tackled this challenge nicely having no earlier precedence of application of this particular (HPSTCO) hybrid optimization mechanism. Therefore, it can be said that this paper introduces a new metaheuristics in DED with significant results.
The current study is done on four different test cases of DED involving 5, 10 and 30 generating units: Case 1 5-unit system with valve point effects, ramp rate constraints, prohibited operating zones and transmission losses Case 2 10-unit system with valve point effects, ramp rate constraints and transmission losses Case 3 10-unit system with valve point effects and ramp rate constraints without transmission losses Case 4 30-unit system with valve point effects and ramp rate constraints without transmission losses.

Organization
The rest of this paper is organized as follows: The "Model" section presents the DED problem formulation with all constraints and limitations. The "Method" section explains the basic concepts and searching principle of each of PSO and TCO methods. In the same section, the proposed HPSTCO algorithm is introduced and described with flowchart and details of the stages. Case studies, simulation results, result analysis and comparison are presented in "Result and discussion" section. Finally, "Conclusion" section draws some concluding remarks on the limitation of the present work and future scope of research.

Model
A comprehensive study of basic DED problem is done here. A non-smooth, non-convex, non-differentiable single-and multi-objective multi-constraint model of ED problem is formulized in this section.

Objective function
The objective function of DED problem, which is to minimize the total production cost over the operating horizon, can be written as: where C T (in$/h) is the total generation cost, C i,t is the generation cost of ith unit at time t, n is the number of dispatch-able power generation units; here, n = 5, 10 and 30, and P i,t (in MW) is the power output of ith unit at time t. T is the total number of hours from operational point of view. The basic ELD objective function is represented by a non-smooth curve (quadratic polynomial) with VPL effect (ripple effect) modeled with a sinusoidal function as shown in Eq. (2).
where a i (in $/h), b i (in $/MWh) and c i (in $/MW 2 h) are the cost coefficients of the ith unit, and e i (in $/h) and f i (in 1/MW) are the VPL coefficients of the ith unit. P min , i (in MW) is the minimum generation capacity limit of unit i. In the generation cost function, the term e i sin(f i (P min i − P i )) represents the VPL effect.
The objective function (Eq. 1) of the DED problem should be minimized subject to the following constraints.

Real power balance constraint
In Eq. (3), P i (in MW) is the power generated by the ith unit, P D (in MW) is the total load demand and P L (in MW) is the total transmission loss of the system at time t. P L is computed using B-coefficients that can be calculated by using Kron's loss formula known as B-matrix coefficients. In this work, B-matrix coefficients method is used to calculate system loss, as follows:

Generator capacity constraint
The generator power output (P i ) of ith generator is within minimum power P min i and maximum power P max i (in MW).

Ramp rate limit (RRL)
A production unit, which is used for generating power P i0 , can increase or decrease its active power output (P i,t ) within upper ramp rate (UR i ) limit (in MW/h) and down ramp rate (DR i ) limits (in MW/h) as shown in Eqs. (6) and (7).

Prohibited operating zone (POZ)
Prohibited operating zone demarcates the scope of active power output of a generator which is otherwise affected due to the technical operation of shaft (unreasonable vibrations of bearing). Usually, modification of power is not allowed in the prohibited spans. The allowable operating range of a generator is given as in Eq. (8).
here j is the number of POZs, P upper i,j−1 is the 'upper boundary' and P lower i,j is the 'lower boundary' of the jth POZ of the ith unit. n i is the number of prohibited operation zones of unit i. The main objective of DED is to minimize the generation cost C T and optimize the power generation schedule (P i,t ) as in Eqs. (1) or (2) subject to satisfying the constraints in Eqs. (3) to (9) used with different combinations in different test cases.
For increase in output power:

Method
The hybrid approach HPSTCO taken up in this study comprises two basic metaheuristics, namely particle swarm optimization (PSO) and termite colony optimization (TCO). A brief outline and working principle of these two optimization techniques are first discussed in this section.

Particle swarm optimization (PSO)
PSO is a swarm intelligence technique inspired from social behavior of bird flocking and fish schooling. Birds and fish follow the neighbor that is nearest to the food, when they search for food. Each individual solution in PSO is named as 'particle' and represents a bird or a fish in the search space. Each particle has a position, velocity and fitness value. While they move in the solution space of fitness function, the particles aim to improve their next position based on their past experience and the best position in the swarm. Therefore, every individual is gravitated toward a stochastically weighted average of the previous best position of its own and that of its neighborhood companions [30]. In every iteration of PSO, the position and velocity of every particle is updated and the value of fitness function at its current location is evaluated.
Mathematically, given a swarm of particles, each particle i is associated with a position vector X i = {X i1 , X i2 , … X iD }, which is a feasible solution for an optimal problem in the D-dimensional search space S. Let the previous best position or pbest of a particle i be denoted by X pi and the best position that has ever been found by any particle or gbest be denoted by X gi . At the start of search, all the positions and velocities are initialized randomly. At each iteration, the position vector of each particle i is updated by adding an increment vector or velocity V i = {V i1 , V i2 , …, V iD } as per Eq. (11). The velocity is updated according to Eq. (10): V i k and V i {k + 1} represent velocity vectors for particle i in the previous and current iterations, c 1 and c 2 are two positive constants, and r 1 and r 2 are two random parameters of uniform distribution in range of [0, 1], which limit the velocity of the particle in the coordinate direction. The new location of each particle should be compared with the pbest value. If the new location of the particle is better than the pbest value, then the pbest is updated for the new location. Otherwise the original value of pbest is stored unchanged. The new global optimum solution gbest is updated according to the gbest of the new particle swarm. This iterative process will continue until a stop criterion is satisfied or maximum number of iterations has been done. Eventually, the particle swarm will converge to the global optimum solution.

Termite colony optimization
TCO is an optimization method inspired from intelligent behaviors of termites during their mound structure building process. Initially the termites arbitrarily search for soil pallets, and after finding it, they deposit it on the mound. Later on, the termites move on the basis of observed trail of pheromone (a chemical) that they deposit on the path on returning after depositing soil pallets on the termite mound. The pheromone acts as attractive stimulus to other members of the colony to follow smaller paths with higher intensities as it is a volatile chemical that evaporates with time. Termites that travel the shortest path reinforce this path with more amount of pheromone, thereby helping others to follow them.
Assuming that the size of the termite population M is within the D-dimensional search space, the position of the ith termite is denoted by � . . , X iD } which indicates a possible solution of an optimization problem. The cost/fitness function value for each position X i is fit ( X i ) which represents amount of pheromone deposited on a hill. The basic steps of TCO can be summarized as follows: 1. Initialize the population as weights, position of termites and the number of iterations. (Every termite has its distinct random position, velocity, desirability and rate of evaporation of pheromone). 2. Evaluate the fitness function value for each termite. 3. Determine the best position and evaporation rate of pheromone of each termite. 4. Determine the position of the best termite. 5. Update the evaporation rate of pheromone, velocity and position of each termite. 6. Stop if the condition of optimization is satisfied. If not, repeat from step 2.
If τ i {t − 1} and τ i {t} stand for the pheromone level at the current and previous locations, respectively, of ith termite, then the pheromone updates rule states: where ρ is the evaporation rate of pheromone taken in the range of [0-1].
After updating the pheromone level, each termite adjusts its route and moves to a new location. Therefore, termite movement is a function of pheromone level at the visited location and the distance between a termite location and the visited locations. Now there are two possible directions of movement: if there is no previously visited location (by the swarm) in the neighborhood of a termite, it moves randomly; if there are one or more visited locations, then the termite selects the location with highest level of pheromone and moves to that position. When the termite moves randomly to search a new gainful position, then position is updated as: here X i (t − 1) and X i (t) represent the current and new position of the termite, respectively; R w is a random walk function of current position and radius of search. When the termite moves toward a gainful position or best local position B i having higher level of pheromone compared to the current position, then the position is updated as: here 1 < w b ≤ 2 and 0 < r b < 1 probabilistically controls the attraction of the termite toward local best position.

Hybrid of PSO and TCO (HPSTCO)
The present work adopts a hybrid of PSO and TCO algorithm (call it as HPSTCO), expecting their usefulness in solving DED problems would be enhanced when used as a combination in complementary mode. The HPSTCO exploits the global search potential of the PSO along with the local search potential of TCO in a given search space. While PSO iterations produce globally distributed solutions (overlooking the localized search space around each global solution), its hybrid partner TCO complements PSO by exploring in more detail any potential localized solution. The solutions obtained by the PSO iterations are fed to the TCO iterations in order to gravitate more termites toward gainful positions. Again the solution found by the termites in TCO updates the positions of the corresponding particles, thereby giving a good starting point of the particles in the global search space.
The basic input parameters of HPSTCO are: maximum number of iterations (max_ iter), population size (s), number of PSO iterations (n 1 ), number of TCO iterations (n 2 ) and number of solutions which are fed from PSO (TCO) to the TCO (PSO) at the switching time (η). The parameter n 1 (n 2 ), respectively, shows how many times PSO iterations (TCO iterations) should be executed before a switching time, implying n 1 iterations of PSO are followed by n 2 iterations of TCO.
The pseudocode of the hybrid algorithm is given in Fig. 1.

Initialization
The HPSTCO algorithm starts with PSO iterations with n number of particles placed in random position in the solution space. A position is a candidate for the priority list P = (p 1 , p 2 , … p n ). Each element of the list represents an activity, and its corresponding value shows the priority of that activity. Hence, the position vector X i = {X i1 , X i2 , …, X iD } of each individual i represents the priority values of n activities. A solution space of priorities will be created where the lower and upper bounds will be defined as Lb = 0:0 and Ub = 1.0. The value of each element must be limited to [Lb, Ub].

Global search
Each particle presents a possible schedule for the DED problem. The velocity of the particles is updated by Eq. (13) which is a modified form of Eq. (10) of original PSO: where 0 < γ < 1 is a constriction factor that improves the convergence speed. The position of each particle is updated by considering its current position and pbest and gbest values are determined by calculating the fitness of the proposed schedules.

Switching from global to local search
The HPSTCO method simply switches from PSO to TCO, while switching a part of solutions found by the PSO that is passed to the TCO. Each solution determines the start position of a termite in the next iteration of TCO. Basically each particle switches its type as termite.

Local search
The TCO uses the solutions which are passed from PSO as the start positions of its termites. Next, TCO tries to find improved solutions in the local neighborhoods of those solutions (now the termites). To determine the neighborhood for each termite, the Euclidian distances of all termites from the candidate termite are calculated. If the distance is smaller than a threshold, the corresponding termite is considered as a neighbor of the candidate termite. The threshold value is dynamically adjusted, gradually decreasing as the algorithm proceeds. The termite with no neighbor moves randomly following Eq. (13); the termite having one or more neighbors selects one of them randomly as its neighbor and updates its position following Eq. (14).

Switching from local to global search
In this phase, each termite switches its type as particle. The solution found by termites updates the positions of the corresponding particles in the PSO. The earlier best position of each particle (pbest) and the global best position of the entire swarm is updated accordingly. The updated fitness of the new solution for a particle is compared to its previous fitness value; if the new fitness value is better, it will be considered as the new pbest. Similarly, the gbest position is compared with this new pbest position, and if the later has better fitness compared to the gbest, then the gbest value is updated with the current pbest.

Constraint handling
In each cycle of HPSTCO, a new population of feasible and infeasible solutions is generated. An infeasible solution is the one which violates the constraints of the problem. After detection of an infeasible solution, it is recovered as a feasible solution. The activity which violates the constraints is changed with the next activity (in the activity list) with lesser priority, and the constraint handling process is applied on the new activity list. This process is iterated until the infeasible solution is converted to a feasible solution.

Result and discussion
In order to review the effectiveness of HPSTCO, it is applied to solve the DED problem on three test systems having 5, 10 and 30 generators, considering valve point loading effect. The algorithm has been coded using MATLAB and implemented on a 64-bit PC with the detailed settings as follows: Hardware CPU: Intel ® Core ™ i5-6200U, frequency: 2.30 GHz, RAM: 8.0 GB, hard drive: 500 GB Software Operating system: Windows 10, package: MATLAB 8.1 (R2014a).
The values of the input parameters of the algorithm are depicted in Table 1.
The simulation in MATLAB is done on four different test cases of DED involving 5, 10 and 30 generating units: Case 1 5-unit system with valve point effects, ramp rate constraints, prohibited operating zones and transmission losses Case 2 10-unit system with valve point effects, ramp rate constraints and transmission losses Case 3 10-unit system with valve point effects and ramp rate constraints without transmission losses Case 4 30-unit system with valve point effects and ramp rate constraints without transmission losses

Test case 1: 5-unit system
In this test system, the valve point loading effects, ramp rate constraints, prohibited operating zones, transmission and generation limits have been considered. The essential input data of the 5-unit system are enlisted in Table 2 [22] that includes prohibited zones of units 1 to unit 5. These zones result in two disjoint subregions for each of units 1, 2, 3, 4 and 5. The B-coefficients matrix used for calculating power loss is given in Table 3. The load demand of the system is separated into 24 dispatch intervals of a day as shown in Table 4. The population size is 50. The fuel cost and transmission losses obtained by the HPSTCO technique are 42,151.3377 $/day and 194.3182 MW, respectively, as shown in Table 5. The graphical representation of Table 5 is shown in Fig. 2. Table 6 shows the comparison results for the fuel cost obtained for 5-unit DED system by HPSTCO with other hybrid methods as reported in the literature. Table 6 shows that the minimum cost yielded by CMIWO [36], MGDE [37], BBOSB [29], MILP-IPM [28], HIGA [25], BBPSO [26], LDISS-2 [24]    43,213$/day, respectively, whereas the cost for HPSTCO is 42,151.3377$ only. The average execution time required for one complete solution was 0.98 min till eighth iteration, and thereafter, the convergence curve becomes a straight line, which is acceptable for DED solutions, though it is not the least in comparison to the time taken by other methods. The convergence characteristic of the proposed algorithm is depicted in Fig. 3.

Test case 2: 10-unit system
In this test case, the valve point effect, ramp rate constraints, transmission and generation limits are considered. The basic input data of the 10-unit system are listed in Table 7.
The B-matrix coefficients (per MW) for calculating power loss are given in Table 8. The load demand of the system is divided into 24 dispatch intervals as shown in Table 9.
Results obtained by MATLAB simulation are presented in Table 10. The graphical representation of Table 10 is shown in Fig. 4. From Table 11 that compares the output of HPSTCO with that of other recently published hybrid methods such as hybrid MILP-IPM [28], hybrid BBOSB [29], HIGA [25], hybrid LDISS [24] and hybrid EP-SQP [32], it is found that the cost (in $/day) yielded by these methods for 10-unit system is 1,040,676, 1,038,362.014, 1,041,087.802, 1,039,083 and 1,035,748, respectively. In comparison, cost and loss yielded by HPSTCO are 1,035,730.203 $ and 811.6073 MW only which is the least among all. The average execution time required for one complete solution was 1.85 min, which is not the least of all but less than many DED solutions. The convergence characteristic of the HPSTCO is depicted in Fig. 5 which shows that the result converged after 50 generations and 1.85 min.

Test case 3: 10-unit system without transmission loss
Unlike test case 2, here a 10-unit system is considered without transmission loss. Like test case 2, valve point effect, ramp rate constraint and generation limits are considered. Input data or DED parameters of the 10-unit system are sane as listed in Table 7. The load demand of the system is divided into 24 dispatch intervals same as shown in Table 9. Results of best generation schedule at each hourly interval as obtained through MATLAB simulation are presented in Table 12. The fuel cost yielded by the HPSTCO method is 1,015,438.967 $/day. The graphical representation of Table 12 is shown in Fig. 6.

Test case 4: 30-unit system
In this test case, input data [15] are obtained by tripling the data of 10-unit system given in Tables 7 and 8. The load demand of the system as divided into 24 dispatch intervals is given in Table 14. In this case, the VPL effects, ramp rate constraints and generation limits are considered. DED results obtained by MATLAB simulation are presented in      Table 15, the graphical representation of which is shown in Fig. 8. The fuel cost obtained by the proposed method is 1,051,964.4$/day. In Table 16, the simulation result of proposed HPSTCO is compared with other recently published hybrid methods, namely BBOSB [29], FA-PSOTVAC [27], BBPSO [26], HIGA [25], LDISS [24], HHS [15] and hybrid EP-SQP [17] Figure 9 shows that the technique takes only nine iterations to reach to steady state which is a good convergence characteristics compared to other methods.

Conclusion
In this paper, HPSTCO has been taken up as a cost minimization and schedule optimization method for 24-h time interval in four test cases representing small-to mediumscale thermal power generation system. Such a hybrid method was never implemented before for dynamic emission dispatch. A synergistic combination of two popular techniques for optimization has been able to mitigate the limitations of the individual techniques. Besides improving the convergence rate, the exploration of neighborhood area for finding local optima has bettered. The hybrid has overcome the problem of convergence to local optima and yields a good globally optimal solution. From the trial runs of the test cases, it can be concluded that HPSTCO is reliable, robust and can consistently provide high-quality solutions of DED considering practical operational constraints, such us valve point effects and multiple fuel changes. The convergence characteristics of HPSTCO are also quite acceptable though not one of the best.
The performance of HPSTCO in terms of cost minimization and dispatch schedule optimization when compared with different other hybrids is found to be quite competitive and can be safely used as an effective metaheuristic for small to medium scale, simple to complex DED problems. In future, this hybrid method can be used to solve the problem of dynamic economic emission dispatch (DEED) problem, multi-objective economic dispatch (MOED) and multi-objective economic emission dispatch (MEED) problem and multi-area economic dispatch (MAED) problem of large-scale power generation system. The HPSTCO method can also be applied to find the impact on optimum dispatch problem of renewable energy like solar and wind energy.       Optimization algorithm: particle swarm optimization (PSO) X i : position of a particle i; Vi k: velocity vectors for particle i in the previous iterations; � Vi {k + 1}: velocity vectors for particle i in the current iterations; D: dimension; S: d-dimensional search space; c 1 , c 2 : positive constants; r 1 , r 2 : random parameters of uniform distribution; X pi : local best (pbest) position of a particle i; X gi : global best (gbest) position of a particle i; V i : velocity of particle i; w: inertia weight.
Termite colony optimization (TCO) X i : cost/fitness function value for each position of the termite; M: size of the termite population; D: dimension; fit(X i ): fitness function value for each position of termite; R w : random walk function of current position radius of search; B i : best local position of termite; τ i {t − 1}: pheromone level of ith termite at the current locations; τ i {t}: pheromone level of ith termite at the previous locations; ρ: evaporation rate of pheromone; w b , r b : probabilistically controls parameters for attracting the termite toward local best position; X max , X min : maximum and minimum limit of search space along a dimension.