Complexity / 2018 / Article

Research Article | Open Access

Volume 2018 |Article ID 7267593 | 23 pages |

Improved Firefly Algorithm: A Novel Method for Optimal Operation of Thermal Generating Units

Academic Editor: Michele Scarpiniti
Received08 Nov 2017
Revised20 Apr 2018
Accepted24 May 2018
Published09 Jul 2018


This paper presents a novel improved firefly algorithm (IFA) to deal the problem of the optimal operation of thermal generating units (OOTGU) with the purpose of reducing the total electricity generation fuel cost. The proposed IFA is developed based on combining three improvements. The first is to be based on the radius between two solutions, the second is updated step size for each considered solution based on different new equations, and the third is to slightly modify a formula producing new solutions by using normally distributed random numbers and canceling uniform random numbers of conventional firefly algorithm (FA). The effect of each proposed improvement on IFA is investigated by executing five benchmark functions and two different systems. The performance of IFA is investigated on six other study cases consisting of different types of objective function and complex level of constraints. The objective function considers single fuel with quadratic form and nonconvex form, and multifuels with the sum of several quadratic and nonconvex functions while a set of constraints taken into account are power loss, prohibited zone, ramp rate limit, spinning reserve, and all constraints in transmission power networks. The obtained results indicate the proposed improvements in terms of high optimal solution quality, stabilization of search ability, and fast convergence compared with FA. In addition, the comparisons with other methods also lead to a conclusion that the proposed method is a very promising optimization tool for systems with quadratic fuel cost function and with complicated constraints.

1. Introduction

The optimal operation of thermal generating units (OOTGU) has been widely concerned in the power system operation field due to its significant important role. In fact, thermal units are using a huge amount of fossil fuel for generating electricity while primary fuels are expensive and will be exhausted in the near future. The objective of the OOTGU problem is to determine the generated active power of thermal generating units so that the total fuel cost can be reduced as much as possible; while all the constraints of thermal generating units and of transmission power networks are taken into account. The problem becomes more and more complicated since more and more complicated constraints of generator are considered such as ramp rate, prohibited zone, voltage limitations, and generation limits. In addition, constraints related to transmission power network are not simple conditions to be exactly met such as the active power spinning reserve, active and reactive power balance, limitation of transmission lines, limitations of transformer tap, and limitation of capacitor bank’s reactive power.

There have been a huge number of studies considering OOTGU as a main problem. The studies have applied conventional methods such as deterministic algorithms and new methods. Most deterministic methods are based on gradient, derivative, Lagrange optimization function while the new methods are based on metaheuristic algorithms, modified versions of metaheuristic algorithms, and combination of two different metaheuristic algorithms. In addition, the metaheuristic methods integrated to other methods dealing with constraints of problem are also effectively combined optimization tools.

Deterministic methods have been applied for the problem such as Lagrange optimization function-based method (LM) [13], dynamic programming (DP) [1], linear programming method (LPM) [4], Hierarchical method (HM) [5], Hopfield neural networks (HNN) [3, 6], improved Hopfield neural network method (IHNNM) [7], augmented Lagrange Hopfield network (ALHN) [8], and enhanced augmented Lagrange Hopfield network (EALHN) [9, 10]. Among the methods, LM was the first optimization tool for solving the problem of OOTGU with quadratic fuel cost functions and the balance constraint of generated power and required power such as load demand and power losses. The obtained results from the method were acceptable since it could deal with the constraint of active power balance constraint and maximum error was very small. DP also achieved results as good as those from LM. However, DP required a huge number of iterations for search process, especially for large-scale systems with large number of units. LPM and HM have applied approximation techniques to simplify the complex of constraints and nonsmooth fuel cost functions and then LM was acted as an optimization search algorithm. Consequently, the applicability of the two methods was restricted for complicated systems with many nonlinear constraints and nonsmooth fuel cost functions. HNN and its improved versions in [79] consisting of IHNNM, ALHN, and EALHN have become more effective than LM, DP, LPM, and HM in dealing with large-scale systems with faster search process and better optimal solutions. The method group is mainly based on energy function and Hopfield neural network. However, HNN variants did not use Lagrange optimization function but they established energy function at the beginning and then Hopfield neural network was applied. On the contrary, ALHN variants have constructed an augmented Lagrange optimization function first and then the function was converted into energy function with the present of inverse sigmoid function. The ALHN variants could overcome several drawbacks of HNN variants such as high oscillation, high error, and low optimal solutions for complicated systems. In general, these methods in the deterministic group have the same limitations of application for systems considering nonconvex fuel cost function and nondifferentiable functions.

New methods based on metaheuristic algorithms have been widely and successfully applied for the considered OOTGU problem even for very complicated systems with nonconvex fuel cost functions and constraints taking prohibited zones, active power spinning reserve, and ramp rate into account. Differential evolution variants have been applied for solving the problem in [1113] consisting of conventional differential evolution (DE) [11], hybrid dynamic programming integer-coded (HDP-ICDE) [12], and colonial competitive differential evolution (CCDE) [13]. DE together with GA and PSO has been implemented for different cases with different constraints with the superiority of DE over genetic algorithm (GA) and particle swarm optimization (PSO) as a result. HDP-ICDE has used integer-coded differential evolution (ICDE) as an optimization tool for searching solutions meanwhile dynamic programming has been used as an evaluation tool for computing fitness function value of found solutions. The method could not overcome drawbacks of ICDE such as high number of solution evaluation and falling into local search. CCDE has changed the mutation operation by using different models based on so-far best solutions and replacement of worst solutions by more promising solutions that aim to improve the balance between local search and global search as well as to retain the best solutions. This method could show its superiority over DE and other versions of DE, but the test systems have been used and the result comparisons from the tests could not show its better performance than other methods because its cost and other methods’ cost were equal. The comparisons of search speed have not been carried out. The performance of cuckoo search algorithm (CSA) [14, 15] and improved cuckoo search algorithm (ICSA) [16] have been shown clearly since different types of constraints were considered and wide comparisons were implemented. However, the speed comparisons have been made only via execution time, leading to not good comparison criteria for conclusion. Different test systems considering constraints of transmission power network have been used to test the efficiency of PSO variants in [1719]. PSO with pseudo-gradient and constriction factor (PG-CF-PSO) [17] has determined the best direction for updating velocity by using pseudo-gradient theory and used constriction factor to find good search zone. It has been demonstrated to be more effective than PSO via three IEEE systems with 30, 57, and 118 buses. New adaptive particle swarm optimization (NAPSO) [18] has integrated both mutation operation and adaptive PSO for avoiding local optimal solution. The method has used different fuzzy and self-adaptive techniques for updating parameter and solutions. Each proposed technique has been investigated by testing different PSO methods with each modification such as NAPSO without mutation (NAPSO1), NAPSO without fuzzy (NAPSO2), and NAPSO without self-adaptive (NAPSO3). These methods have had worse results than NAPSO but better results than standard PSO. Combination of particle swarm optimization and tabu search algorithm (PSO-TSA) [19] has used the ability of PSO for local search based on three different phases and TSA for global search by adjusting obtained solutions of PSO. The method has been compared to three difference implemented methods such as GA, PSO, and TSA, and it can be shown better objective function. In spite of the superiority over the three methods, PSO-TSA has been more complicated, included more control parameters, and used high number of produced new solutions. Different GA variants have been proposed for dealing with different set of constraints of the considered problem such as real-coded genetic algorithm (RCGA) [20], RCGA with arithmetic-average-bound crossover and wavelet (NRCGA) [20], hybrid real-coded genetic algorithm (HRCGA) [21], improved GA (IGA) [22], and IGA with updated multiplier (IGAMU) [22]. These methods have been developed by improving GA, using RCGA and other improved crossover operation, and combining improved version of GA and multiplier updating technique. Generally, more complicated modifications lead to more improvement of GA variants. However, there have been more drawbacks such as higher number of control parameters, difficult task of selection of such factors, and complicated implementation when combining two or more methods. The different versions of krill herd algorithm (KHA) [23] have been used for solving different cases of the considered problem such as KHA without using genetic operation (KHA1), KHA using crossover operation (KHA2), KHA using mutation operation (KHA3), and KHA using crossover operation and mutation operation (KHA4). Several disadvantages of these KHA methods have been pointed out in [24] such as low convergence speed and low success rate. Thus, opposition-based krill herd algorithm (OBKHA) [24] has been developed by applying learning technique based on opposition and it has been compared to KHA methods. As derived from the comparisons, OBKHA could be considered better than KHA methods about optimal solutions reflected via less minimum cost and stable search ability reflected via less average cost, less maximum cost, and less standard deviation; however, the convergence speed and the success rate have not been investigated since there were no comparisons of used population and used iterations. Biogeography-based optimization algorithm (BBOA) has been employed for minimizing cost and emission in [25] and for minimizing cost in [26]. BBOA has some characteristics in common with other metaheuristic algorithms consisting of mutation and selection like GA, DE, and evolution programming, but this method has much more number of control parameters such as the probability of habitat modification and mutation, the rate of mutation and immigration, and lower and upper boundaries of probability. The results reported in [26] could show good performance of BBOA via comparisons with several existing methods on four test cases, but the keywords for determining the best control parameters were missed. Clearly, the application of BBOA to real power system with complicated constraints in transmission power network was limited and studies about BBOA for complicated systems have not been carried so far. Grey wolf optimization algorithm (GWOA) [27, 28] has been applied for simple systems without complicated constraints and small number of units. Thus, the demonstration of the methods’ performance has not been persuasive. On the contrary, crisscross optimization algorithm (CSOA) [29] and exchange market algorithm (EMA) [30] have had more result comparisons in the implementation of large-scale systems and complicated constraints. The ant lion optimization algorithm (ALOA) [31] was also a new method with few applications, and the method in the study has not shown its potential search ability persuasively due to simple test employment. In addition, many new methods have been developed for solving the problem such as modified symbiotic organisms search algorithm (MSOSA) [32], mine blast algorithm (MBA) [33], clonal algorithm (CA) [34], mathematical programming algorithm (MPA) [35], improved quantum-inspired evolutionary algorithm (IQIEA) [36], cuckoo optimization algorithm (COA) [37], improved colliding bodies optimization algorithm (ICBOA) [38], flower pollination algorithm (FPA) [39], natural updated harmony search (NUHS) [40], lightning flash algorithm (LFA) [41, 42], moth swarm algorithm (MSA) [43], and orthogonal learning competitive swarm optimization algorithm (OLCSOA) [44]. Among the remaining methods, LFA, MSA, ICBOA, and OLCSOA were developed in the last two years and they were considered better than most existing methods for considered test cases. LFA and OLCSOA have been implemented for systems with different complex levels of thermal generating units such as nonsmooth fuel cost functions with valve-point loading effects, multifuel types, and prohibited zone meanwhile challenges of transmission power network constraints have not been taken into account. On the contrary, MSA has been focused on dealing with the constraint complex rather than the challenges of fuel cost function of thermal generating units. In general, metaheuristic algorithms have been widely and successfully applied for different test systems and their superiorities have been pointed out mainly based on fuel cost function comparisons. The performance of ICBOA has been measured by testing on different study cases with different types of fuel cost function and the consideration of all constraints of transmission power networks. ICBOA together with other methods such as artificial bee colony (ABC), DE, PSO, BBO, standard colliding bodies optimization algorithm (CBO), and enhanced CBO have been measured by testing on different study cases with different types of fuel cost function and the consideration of all constraints of transmission power networks. ICBO has been pointed out obtaining better optimal solutions than all methods, but the further comparisons of convergence speed have not been carried out and discussed.

The firefly algorithm (FA) is also a population-based metaheuristic algorithm similar to PSO, DE, CSA, and so on built by Yang in 2008 for solving optimization problems [45]. The configuration of FA consists of three procedures of updating the distance between two considered fireflies, updating a step size, and updating new solutions. FA has shown its superiority over other traditional algorithms consisting of GA and PSO [46] by the comparisons of different benchmark optimization functions. In 2011, Yang and his coworkers applied the algorithm for solving the considered OOTGU problem [47] and they have stated that the method has been efficient for the problem by carrying out result comparisons via testing on three different systems. However, the performance of the method has not been demonstrated for large-scale systems with large number of generators and highly complicated constraints in transmission networks. Thus, there were several improvement versions of the method for dealing with the problem and the comparisons with it were also done for evaluation. An improved firefly algorithm (IFA) has been proposed by Kazemzadeh-Parsi in [48] by applying three different modifications. The first modification was to use the better solutions of the previous iteration for replacing the worst solutions in the current iteration. The second modification was to replace low solutions by randomly produced solutions, and the third modification was to move the worse solutions only to one representative solution instead of all better solutions. The representative solution is determined by using the average solution of all better solutions. The memetic firefly algorithm (MFA) in [49] has focused on the balance of exploration acting as global search and exploitation acting as local search by using adaptive attractiveness β and adaptive step parameter α with respect to the change of iteration in formula updating new solutions for each lower quality solution. Another improved FA with using adaptive step parameter α (ASPFA) was proposed in [50] by suggesting an adaptive formula for updating step parameter α based on current iteration and the maximum number of iterations. The chaotic firefly algorithm (CFA) was developed in [51] for solving the considered OOTGU problem with different test cases of fuel cost function and constraints. The method has used chaotic distribution to produce the values of attractiveness β and step parameter α instead of using random distribution in FA. Moustafa et al. [52] have implemented conventional FA, IFA, MFA, and ASPFA on 3-unit system and 6-unit system for comparison. The result comparisons have indicated that IFA with three modifications was the best one and FA was the worst one among four FA variants. However, there was no comparison between these methods with others. On the contrary, CFA in [51] has been tested on many tests with more complicated objective functions considering valve-point effects and multifuels, and more complicated constraints considering POZ, ramp rate, and active power spinning reserve. CFA has been demonstrated to be more effective than many methods such as PSO, ABC, DE, and GA variants. Another modified firefly algorithm (α-MFA) based on the modification on distance between two fireflies and the modification on generation of step parameter α was proposed in [53] for dealing with OOTGU problem with 3-, 6-, 13-, and 15-unit systems considering valve-point loading effects, ramp rate constraints, and POZ constraints. The five improved versions of FA have been demonstrated to be better than conventional FA and other methods. However, the superiority of the methods was not shown clearly since it has not been implemented for complicated systems taking all constraints of transmission power networks and generators into consideration.

In the paper, we propose three modifications on the FA in order to tackle several disadvantages of FA such as premature convergence to local optimum solution and impossibility of jumping out of the search zone with many local optimum solutions. In the first modification, we propose a new formula to update radius between a considered firefly (a solution) to another firefly (another solution) with lower fitness function than the considered solution. The proposed radius based on and the best solution become more effective better than that based on and used in FA. In the second modification, we propose a new algorithm for producing new solutions of an old solution by suggesting two models for the updated step size. A large or a smaller updated step size will be used to diversify new solutions and to avoid converging to a local optimum and trapping into search zone with many local optimums. In the third modification, uniform distribution is replaced with normal distribution aiming at diversifying the search zone. As a result, the improvement of the new algorithm is highly considerable compared to FA. The application of each modification and three modifications is evaluated by testing on difference of four systems with other six cases. These test cases consider thermal generating units using single fuel and multifuels, convex and nonconvex objective function, and all constraints in real power systems. The main contributions of the paper can be summarized as follows: (i)Carry out three modifications on the new solution produced phase of conventional FA and point out the specific impact of each modification on the obtained results.(ii)The proposed method has few control parameters, and the parameters are easily tuned. The proposed method is more effective than FA, but the implementation of the proposed method is as simple as that of FA.(iii)Thoroughly present the application of the proposed IFA method for solving different systems with different constraints in optimal operation of power systems. The exact selection of control variables can help the proposed method satisfy all considered constraints.(iv)Develop basic application of the proposed method for optimization problem and lead to conclusion that the proposed method is effective for systems with complicated constraints but it works less effectively for systems considering nonconvex objective functions.

The rest of the paper is organized as follows. Section 2 describes the problem formulation. Section 3 presents the firefly algorithms. Section 4 shows the implementation of IFA for OOTGU problem. Section 5 presents the results of numerical simulations. Finally, the conclusions and future works are given in Section 6.

2. Problem Formulation

2.1. Objective Function

In OOTGU with single fuel, the fuel cost of each generating unit is expressed as a quadratic function of its power output. The objective of the problem is to minimize the total fuel cost of available units and can be written as follows [1]: where is the fuel cost function of thermal unit and it can be represented in a single quadratic form or piecewise form depending on the number of fuel cost options that thermal plants use. If the thermal units use only one fuel, its form is the second order equation as shown in (2) while the form of when using multiple fuels is a piecewise curve consisting of at least two second order equations. For the case of multifuel options, is shown in (3) [1].

In addition, once valve-point loading effects (VPLE) of thermal units are considered during the process of increase or decrease of power output, (2) and (3) will become more complicated models as follows [20].

For easy understanding, the case of single fuel option and the case of multifuel options neglecting and considering VPLE are depicted in Figures 1 and 2. It is clear that the curves with VPLE are much more complicated.

2.2. Set of Constraints
2.2.1. Generator Constraints

Active and reactive power and voltage of generators are restricted by the following inequality [17]. where and are the lower and upper limitations of reactive power of generator , respectively; and are the lower and upper limitations of voltage of generator , respectively; and and are the operating values of reactive power and voltage of generator , respectively.

It is noted that the constraint (5) is always taken into account in the OOTGU problem while the constraints (6) and (7) are considered only when real transmission power networks are being considered.

2.2.2. Prohibited Operating Zone (POZ) Constraints

In addition to the active power limitation constraint in (5), thermal units are also constrained by prohibited operating zones due to the security requirement. Thus, power output of each thermal unit should be within its lower limit and upper limit and outside the POZ. An example with two prohibited operating zones of thermal units is shown in Figure 3, and a typical thermal unit with prohibited operating zones is mathematically modeled as follows [20].

2.2.3. Spinning Reserve Constraint

In order to avoid lack of energy in case that the largest power unit shutdowns, total active power reserve of all thermal units should be required to be equal to or higher than the largest power unit. The requirement is as below [20]. where is the total active power that the power systems require for spinning reserve and is the active power spinning reserve of thermal generating unit .

2.2.4. Ramp Rate Constraint

During the operation of thermal units, generation increase and generation decrease are constrained by the ramp up and down the rate limit, respectively. Thus, initial power output of units and the ramp rate limits will determine the lowest and the highest power output of unit at the moment. The description can be presented in the following inequalities [20].

2.2.5. Real Power Balance

The total real power output of generating units satisfies total load demand plus system power losses [1]. and the total power loss is calculated using Kron’s (1):

In addition to (11), real power balance is also presented by the following equation since real transmission power networks are considered [17]. where and are the real part and the imaginary part of admittance of a transmission line connecting buses and .

Similarly, reactive power balance is also modeled as follows [17]: where is the generated reactive power of capacitors at bus and constrained by the following inequality [17]: where and are, respectively, the lowest and the highest reactive power of the capacitor bank at bus .

2.2.6. Transformer Tap Constraints

In a real transmission power network, transformers are usually used to step down voltage and supply electricity to lower voltage transmission lines. The voltage of the lower voltage lines can be adjusted by setting the tap of the transformers, so that the following rule is always exactly met [17]. where and are the lower limit and the upper limit of transformer’s tap at bus .

2.2.7. Stably Working Power Network Constraints

In order to operate electric power network continuously, load bus voltage and apparent power flow of transmission lines must follow the following limitations [17]. where and are, respectively, lower voltage and upper voltage of load at bus ; is the apparent power flow capacity of line ; and is the current apparent power flow of line .

3. Firefly Algorithms

3.1. The Conventional Firefly Algorithm

Firefly algorithm, a metaheuristic algorithm, was first built by Yang in 2008 for solving optimization problems [45]. The configuration of FA was developed by the three following ideas: (i)Each firefly attracts another one by its brightness.(ii)The fireflies with higher brightness have higher attractiveness level to other ones.(iii)The fireflies with lower brightness level move to other ones with higher brightness level.

The three behaviors of nature fireflies have inspired Yang in building an optimization algorithm called firefly algorithm. There is a mutual connection between the behaviors of fireflies and the construction of FA. In fact, each firefly corresponding to each optimal solution will own its brightness corresponding to the fitness function of the optimal solutions. The action that the fireflies with darker brightness will look for and get to other fireflies producing higher brightness level is similar to the newly produced solutions based on old solutions with a better fitness function. Consequently, in FA, each old solution can be newly produced several times depending on the comparison of its brightness with other ones. As a result, only one new solution of each old solution is kept based on the comparison of fitness function.

Suppose that each solution is a position of firefly at the current iteration. When the fitness function of solution is higher than that of another solution , the distance between the firefly and is obtained by using the following model.

Then, the updated distance is employed to be substituted into another (19) to calculate a new attractiveness. Then, the new position for the ith considered firefly can be determined corresponding to the generation of a new solution of the ith solution. The procedure of generating a new solution is carried out as (20). where is a random number of solution and β0 is the attractiveness at zero distance and normally set to 1. is a solution having lower fitness function than ; and is an updated step size calculated by the following model.

Equations (18), (19), and (20) continues to be determined for the ith solution until there are no more solutions that have lower fitness function.

In summary, for each solution , we can have either one or higher than one new solution or no new solution depending on the fitness comparison between solution and other solutions among the current population. The statement can be described according to the following term.

For the first term in (22), if the considered solution is also the global best solution, there will be no any new solution generated. For the second case, there will be only one new solution, is generated if the considered solution is the second best solution, and is the global best solution among the population. For other cases, it means is the third best solution or worse than the third best solution and even it is the worst solution, there will be from two new solutions to new solutions . In this regard, the set of new solutions of solution will be evaluated by fitness function value comparison and the best one with the lowest fitness is retained while others are discarded.

The whole description of FA is shown in detail in flowchart of Figure 4.

3.2. The Proposed Improved Firefly Algorithm

FA has been applied for solving different optimization problems such as OOTGU problem and benchmark optimization functions, and its results have been better than those from other popular methods like PSO, GA, and DE. But the method has not met the demand of dealing with large-scale systems and complicated constraints of OOTGU problem, and many studies have worked on improving it [4951, 53]. In fact, FA still suffers from several drawbacks such as premature convergence to local optimum or convergence to near global optimum with high number of iterations [54]. In order to solve the difficulties of using FA, we suggest three improvements of the attractiveness β (namely, that is the radius), updated step size , and random number (rand) in (20). The detail of the three improvements is as follows.

3.2.1. The First Improvement

Instead of using the distance between the considered solution i and another better solution to determine the radius, the best solution is recommended to be used for calculating the radius as shown in the expression.

We use to replace in (18) because the change can produce a more effective scaling factor, which is the attractiveness β. The attractiveness β is a function with respect to the radius as seen in (18). The usefulness of the change is demonstrated in the numerical results where IFA1 (FA with the first improvement) is better than FA in terms of optimal solution reflected via the best solution comparison and the stable solution search ability reflected via the standard deviation comparison.

3.2.2. The Second Proposed Improvement

In the second improvement, we propose a new technique to generate updated step size and obtain lower solution fitness than those of FA. It is clear that the manner of producing an updated step size by using (21) is similar to mutation operation of differential evolution algorithm (DEA) in which β acts as mutation factor ranging from 0 to 2. Some studies before [55, 56] have pointed out the disadvantages of DEA such as low convergence to global optimum or trapping into local optimum easily. Consequently, the proposed improvement aims to tackle the constrictions of FA by using the following expressions.

It is clear that the first updated step size is less than the second updated step size . Thus, the use purpose of each one is different, in which the former narrows the search zone nearby old solutions while the latter can expand exploitation to avoid falling into the same solutions. Derived from the difference, a condition of using either or is suggested based on the fitness ratios as follows:

and are not fixed at constant values, but they are not control parameters and their values can be exactly calculated by

It is clear that there are two possibilities for comparison between and . is either less or higher than . Since the first assumption occurs, it means that solution is better than the average solution of the whole population and if the second assumption happens, solution is worse than the mentioned average solution. In the case that solution is better than the mentioned average solution, solution is very close to current good solutions, which are better than the mentioned average solution. So if we use the first updated step size , the change from current solution is small and the new solution can be very close to or the same as the current good solutions. Clearly, the second updated step size is much larger than and it can push the obtained new solution far away from the current good solutions, avoiding jumping to local optimum zone and premature convergence. On the contrary, since solution is worse than the mentioned average solution, it is not close to the current good solutions and the search around it with a not large step size is better. The impact of the proposed improvement is significant for obtaining good result of the proposed method via the observation of comparison between FA and IFA2 (FA with the second improvement) in numerical results.

3.2.3. The Third Proposed Improvement

In the third proposed improvement, we suggest using normal distribution (randn) instead of using uniform distribution (rand) like the conventional FA. The suggestion aim is to diversify the search zone as used in Levy flights of Cuckoo search algorithm [46]. As a result, the equation of producing new solutions is as follows:

The effectiveness of the proposed improvement is demonstrated by comparing IFA3 (FA with the third modification) with FA in numerical results.

The implementation of the proposed IFA for a general optimization problem is similar to the flowchart shown in Figure 4 of FA; however, the difference between the two considered algorithms is the way to produce new solutions (shown in Step 5). The whole search procedure of the proposed IFA is described in Algorithm 1.

Set control parameters to values
Randomly initialize population
Calculate fitness of the whole population
Determine the best solution and the worst solution
Set ;
Set ;
Calculate , and
    Determine and
    Determine attractiveness β and new solution
    Calculate fitness function of
    End    % if
  End      % for
Compare each old solution and each new solution to keep better one.
Compare these fitness function to choose the best fitness function
Set with the lowest fitness to
End        % for
Determine the best solution and the worst solution
End        % while

4. The Implementation of IFA for OOTGU Problem

4.1. Dealing with Load Demand-Supply Balance and Other Constraints

In order to deal all constraints, decision variables and dependent variables are chosen thoroughly. For the case that constraints of real transmission power network are not taken into account, power output of the first () thermal units is included in the position of each firefly as shown in where

As a result, the load demand-supply balance constraint can be dealt successfully by using the dependent variable obtained by the below model [15].

For complicated case with the consideration of all transmission power network constraints, active power output of all generators excluding at slack bus, voltage of all generators, tap setting of all transformers, and reactive power output of all capacitors are chosen as decision variables. Thus, , , and are not the same as (28), (30), and (31) and can be obtained as follows:

Note that is not included in the set of decision variable because it is supposed to be the active power of generator located at slack bus. Other dependent variables are obtained by running optimal power flow program, Mathpower.

4.2. Penalizing the Violation of Dependent Variables

There is a possibility that can violate its limitation, lower than the lowest generation or higher than the highest generation. Therefore, the violation must be controlled and considered in the quality evaluation of solutions. The task is done by calculating the penalty term as indicated in the following equation.

Similarly, when the real transmission power network constraints are running, the third penalty term calculation is carried out by the following model. where the limits related to dependent variables are determined by

4.3. Fitness Function

The fitness function of all solution should be determined to arrange the effectiveness of all the solutions. Consequently, the fitness function considering objective function and the sum of penalty terms is shown in the following model. where is penalty factor used to amplify the violation of the dependent variable [16].

4.4. The Procedure of IFA for OOTGU Problem

The whole search process of the proposed method when solving the considered problem is described in detail as follows.

Step 1. Initialize a population with solutions satisfying the condition of (29).

Step 2. Calculate dependent variables for solution by using (32) or running Mathpower optimal power flow program.

Step 3. Determine penalty terms by using (34), (35), and (36).

Step 4. Find the fuel cost coefficients corresponding to power output of TGU in case of multifuel options.

Step 5. Calculate fitness function by using (38).

Step 6. Assign the solution with lowest fitness into the best one and start the first iteration ().

Step 7. Consider solution .

Step 8. Produce new solutions for solution by using three proposed modifications.

Step 9. Check limits for and perform Steps (2)–(4) above.

Step 10. Calculate fitness function for all new solutions and assign the solution with the lowest fitness into .

Step 11. If , go to the next step. Otherwise, set and back to Step 8.

Step 12. Compare each new solution to each old solution to keep better one.

Step 13. Determine the best solution.

Step 14. If Iter is equal to the highest iteration number , stop the iterative algorithm. Otherwise, set and back to Step 7.

5. Numerical Results

In the paper, the efficiency of proposed IFA method is verified by testing on six main cases with different characteristics of fuel cost function, thermal unit constraints, and power network constraints. The detail of test systems is described as follows.

Case 1. Six thermal generating unit test system with POZ constraint and power loss constraints [57].

Case 2. Fifteen thermal generating unit test system with POZ constraint and without power loss constraints [22].

Case 3. Twenty thermal generating unit test system with power loss constraints [3].

Case 4. Forty thermal generating unit test system with valve-point loading effects [58].

Case 5. Ten thermal generating units with multifuels. Subcase 5.1.Neglect valve-point loading effects and power losses [7].Subcase 5.2.Consider valve-point loading effects and neglect power losses [20].Subcase 5.3.Consider valve-point loading effects and the constraints of POZ, ramp rate limits, and active power spinning reserve [20].

Case 6. The IEEE-118 bus system with all constraints of real transmission power networks [59].

The population size and the highest iteration number selected for implementation of IFA are reported in Table 1. For each study case, each method is run one hundred independent trials by using MATLAB and computer with 4 GB of Ram and 2.4 GHz processor.

CaseNumber of units


5.1. Impact Analysis of Each Proposed Improvement

In this section, the contribution of each proposed modification to the performance of the proposed method is investigated by running conventional FA, FA with the first proposed improvement (IFA1), FA with the second proposed improvement (IFA2), FA with the third proposed improvement (IFA3), and the proposed IFA method on five benchmark functions and two power systems. Such benchmark function is given in Table 2 [17, 60] while the two power systems with 6 units and 15 units considering POZ constraints are taken from Cases 1 and 2.







5.1.1. Testing on Benchmark Functions

In order to implement FA, IFA1, IFA2, IFA3, and the proposed method, the population size and number of iterations are set to 50 and 1000, respectively. In addition, we also employ Wilcoxon signed-rank test [61] and set the significance improvement level to 0.01 for calculating value, value, and value for better comparison between the proposed method and others. As a result, the obtained results in terms of minimum fitness, standard deviation of 50 trial runs, value, value, and value are reported in Table 3. Among the result numbers, the comparison of the best cost is to reflect the best optimal solution and comparison of the standard deviation cost is to reflect the stabilization of search ability. As pointed in [62], value can be either equal to or smaller than 0.01 or higher than 0.01. In case that value of a compared method is higher than 0.01, the significance improvement level of the proposed method over the compared method is not adopted but the significance level is accepted in case such value is less than 0.01 or equal to 0.01 [62]. Table 3 shows that all minimum values and standard deviation of the proposed method are smaller those of other ones meanwhile value, value, and value are, respectively, equal to −6.154, 0, and 0. value = 0 indicates that the proposed method reaches the significance improvement level over other FA variants while value = 0 determines the accuracy of value and value. For other comparisons among FA, IFA1, IFA2, and IFA3, the minimum fitness and standard deviation reports can indicate that IFA2 is the best modification; IFA3 is the second best modification while IFA1 is only better than conventional FA. The further investigation of performance of each FA methods can be seen clearly in the following section since two power systems with 6 units and 15 units considering POZ constraints are employed for testing.

FunctionMethodMin.Std. dev. value value value




IFA0.000203445.1E − 05


5.1.2. Testing on Power Systems with POZ Constraints

In this section, the performance of FA, IFA1, IFA2, IFA3, and the proposed IFA method is illustrated by running on power systems of Study Cases 1 and 2. Furthermore, four other popular methods consisting of particle swarm optimization (PSO) [63], differential evolution (DE) [64], flower pollination algorithm (FPA) [65], and crow search algorithm (CSA) [66] are also implemented for solving two test systems. For running the methods, we set the population and the number of iteration to those values of the proposed method shown in Table 1. In addition, we also tune other control parameters of these methods to their predetermined ranges. For the implementation of PSO, acceleration factors and are set to the range of [0.2, 2.0] with a step size of 0.2. For the implementation of DE, mutation factor and crossover factor are, respectively, set to range of [0.2, 1.2] and [0.2, 1.0] with a step size of 0.2. For the implementation of FPA, probability of using either local search or global search is set to the range of [0.1, 1.0] with a step size of 0.1. For the implementation of CSA, awareness probability and flight length are, respectively, set to the range of [0.2, 1.0] and [0.2, 2.0] with a step size of 0.2. The obtained results from the proposed method and these compared methods such as minimum cost and standard deviation and the results from Wilcoxon signed-rank test such as value, value, and value are reported in Tables 4 and 5 for Study Cases 1 and 2.

MethodBest cost ($)Std. dev. ($) value value value


MethodBest cost ($)Std. dev. ($) value value value


For Case 1, the proposed method has obtained the best results with less best cost than FA, IFA1, IFA2, and IFA3, respectively, by $2.14, $1.843, $0.001, and $0.182. Similarly, the standard deviation cost of IFA is, respectively, less than that of FA, IFA1, IFA2, and IFA3 by $26.765, $18.05, $1.322, and $17.113. Furthermore, value, value, and value also indicate the significance improvement of the proposed method over FA variants. Comparison among IFA1, IFA2, and IFA3 can see that such three improved versions has better performance than FA; meanwhile, IFA2 is the best improvement and IFA1 is the worst improvement. Analysis of Case 2 also gives the same evaluation for FA methods.i

For comparison with PSO, FPA, DE, and CSA, all the costs yielded by the proposed method are smaller than those from PSO, FPA, DE, and CSA. In particular, the standard deviation of the proposed method is approximately zero while this value of other methods is much higher than zero. Furthermore, value, value, and value calculated by setting the significance improvement level to 0.01 also indicate that the proposed method can improve working performance significantly when compared to PSO, FPA, DE, and CSA. Among IFA1, IFA2, and IFA3, only IFA2 outperforms PSO, FPA, DE, and CSA since the best cost and standard deviation cost of IFA2 are smaller than those of PSO, FPA, DE, and CSA. Consequently, it can conclude that the proposed method has high performance for searching optimal solutions while the second modification has the most impact on results.

Besides, the more accurate comparison of the best cost improvement and the standard deviation improvement can be quantitative by using the following equation.

In (39), result of IFA or result of compared method can be the best cost, mean cost, worst cost, or standard deviation cost dependent on decision of evaluation but the two most important evaluations must be used including the best cost improvement and the standard deviation improvement. As result is the best cost, IFA result improvement is the improvement of optimal solution quality of IFA over the compared method. In addition, as the result is standard deviation, the IFA result improvement is the improvement of stabilization of IFA over the compared method. As a result, the improvement of IFA over others in terms of the best optimal solution quality and the stabilization is reported in Table 6 for Cases 1 and 2. As observed from Table 6, the best cost improvement of IFA over FA and other methods is not high but the improvement of standard deviation cost is approximately 100% for the two cases.

CaseCompared methodBest cost improvement (%)Std. dev. improvement (%)



The best cost of fifty independent trials for Cases 1 and 2 obtained by the proposed IFA and FA is, respectively, plotted in Figures 5 and 6. The figures can evaluate the standout of IFA over FA clearly once all minimum costs of IFA are nearly on a flat line and much lower than those of FA, especially for Case 2 while fifty cost values of FA have high fluctuations for the two cases.

5.2. Comparison and Discussion

In order to investigate the further performance of IFA, we carry out the comparisons of results obtained by IFA and many optimization algorithms. In addition to the comparison of the best costs, another comparison criterion is also considered to be the number of fitness evaluations , which is shown in the following equation. where is the number of generations in each iteration. For some optimization algorithms possessing two new solution generations such as CSA [15] and ICSA [16], is 2 while for others possessing one new solution generation like PSO [11] and DE [11], is only 1. For the proposed IFA, there is only one new solution generation in each iteration; thus, is also equal to 1. As a result, the value of is added into each table for the comparison and the indication is that a method with lower is considered to be more efficient if its best cost is also lower or equal.

5.2.1. Comparisons of Test Case 1

Table 7 reports comparisons of the best cost, mean cost, worst cost, and from IFA and other methods for Case 1, and especially the best cost improvement of IFA compared to each one is also reported at the final column. It notes that if the best cost improvement (BCI) in % has plus sign (+), the result from IFA is better than the compared method, and if BCI has minus sign (−), the compared one has better optimal solution than IFA. Besides, if BCI equals zero, the proposed method and compared methods have the same best cost and the same quality of optimal solution. The observation of BCI can lead to an evaluation that the proposed method is slightly better than most methods (the improvement is from 0.0002% to +0.1030%). KHA4 [23], OBKHA [24], and EMA [30] have the same cost as the proposed method, but α-MFA [53] has slightly lower cost than IFA. However, the recalculated cost of α-MFA [53] is $15445.46, which is higher than the cost of the proposed method. The comparison of can indicate that the proposed method is at least 4 times faster than other ones since of IFA is 300 but that of other ones is from 1200 (α-MFA [53]) to 50,000 (EMA [30]). In spite of using lower value of , the mean cost and the highest cost of the proposed method are also much less or approximately equal to those from other ones. As a result, it can result in a conclusion that the proposed method is much effective than other methods when solving test Case 1 with quadratic fuel cost function, POZ constraint, and power loss constraint.

MethodBest cost ($)Mean cost ($)Worst cost ($)BCI (%)

DE [11]15,449.7715,449.8715,449.783600+0.0433
GA [11]15,45915,52415,4693600+0.1030
PSO [11]15,45015,49215,4543600+0.0448
KHA1 [23]15,450.7515,452.8215,455.5010,000+0.0497
KHA2 [23]15,448.2115,450.8315,453.4010,000+0.0333
KHA3 [23]15,445.3615,447.2215,449.6010,000+0.0148
KHA4 [23]15,443.0815,443.1915,443.3010,0000.0000
OBKHA [24]15,443.0815,443.9215,443.330.0000
BBO [26]15,443.1010,000+0.0002
EMA [30]15,443.0715,443.0850,0000.0000
AIS [34]15,4484000+0.0319
NAPSO [18]15,443.7715,443.7715,443.773000+0.0045
NAPSO1 [18]15,449.1415,451.0415,452.923000+0.0393
NAPSO 2 [18]15,444.7715,446.7715,449.773000+0.0110
NAPSO 3 [18]15,449.0815,450.3315,452.783000+0.0389
α-MFA [53]15442.91200−0.0011

5.2.2. Comparisons of Test Case 2

This section compares the performance of the proposed method with other ones when dealing with a larger scale system compared to Case 1. The obtained results from the proposed method and other compared methods are given in Table 8. It is obviously seen that ICSA [16] has the best optimal solution; however, the confirmation of a valid solution is not done for ICSA because there was no optimal solution reported for the case in [16]. In comparison with mean cost and maximum cost of other ones, the proposed method has higher values because we run the proposed method with of 1000 while others were run by a very high value. In case that we increase of the proposed method to 6000 and higher than 6000, the result is that the mean cost and the highest cost are approximately equal to the best cost. The illustration aims to demonstrate the fast convergence of the proposed method compared to others, and it is clear that the proposed method is better than this compared method for Case 2 with 15 units considering quadratic fuel cost function and POZ constraint.

MethodBest cost ($)Mean cost ($)Worst cost ($)

CSA [15]32,544.970432,545.00632,546.67348000
ICSA [16]32,542.559332,543.1688432,546.65746000
IGA [22]32,04.736250,000
IGAMU [22]32,544.99050,000

5.2.3. Comparisons of Test Case 3

In the section, a larger system with 20 units considering quadratic fuel cost function and power loss constraint is employed for comparison. The information of the best cost and is shown in Table 9 for comparison with deterministic methods and metaheuristic methods. The best can see that the proposed method can obtain less cost than FA, Lamda, and HM methods in [3] and BBOA in [26] but slightly higher than cost of CSA and ICSA in [16]. The deviation implies that the proposed method is more effective than lambda, HM, and BBOA, but it is not less effective than methods in [16]. Besides, the proposed method has used much less than other ones including FA. As we use a higher number of equal to 10,000, the best cost obtained by the proposed method is also as good as that of methods in [16]. The numerical comparison can show a good performance of the proposed method when dealing with a system with 20 units considering quadratic fuel cost function and power loss constraint.

MethodCost ($)

Lambda [3]624656.639
HM [3]62456.6341
CSA [16]62456.63310,000
ICSA [16]62456.63310,000
BBOA [26]62456.792620,000

5.2.4. Comparisons of Test Case 4

Test Case 4 is more complicated than above cases due to the consideration of 40 units and nonconvex fuel cost functions. The result comparisons of the case are shown in Table 10. The values of BCI can see that the proposed method has better cost than methods with plus sign “+” of BCI such as BBO, OLCSO, KHA1, KHA2, KHA3, FA [47], and CFA [51]. In comparison with remaining methods, the proposed method has lower quality of solution with higher cost excluding FPA, which did not report optimal solution and had the surprised lowest value. The comparison of also goes against the proposed method in the evaluation of performance. In fact, the proposed method uses a higher number of compared to CCDE, FA, and KHA methods but it uses smaller than the rest of methods including CFA, which is an improved version of FA. The analysis of best cost and can point out that the application of the proposed method is not as effective as the applications with quadratic fuel cost function in Cases 1, 2, and 3 above.

MethodBest cost ($)Mean cost ($)Worst cost ($)BCI (%)

CCDE [13]121,412.58121,413.56121,416.5370,000−0.002
CSA [16]121,412.54121,520.41121810.25160,000−0.002
BBO [26]121,426.95121,508.03121,688.66200,000+0.010
FPA [39]121,074.50121,095.70121,196.302000−0.281
OLCSO [44]121,415.82121,504.05121,460.78400,000+0.001
FA [47]121,415.05121,416.57121,424.5650,000+0.000
KHA1 [23]121,460.42121,468.98121,477.4610,000+0.038
KHA2 [23]121,448.36121,453.68121,461.3910,000+0.028
KHA3 [23]121,423.46121,428.23121,433.5610,000+0.007
KHA4 [23]121,412.60121,413.15121,415.0010,000−0.002
CFA [51]121,415.02121,416.23121,420.67336,000+0.00035

5.2.5. Comparisons of Test Case 5

This study case considers different fuel cost functions with nonsmooth forms such as the sum of several quadratic functions for Subcase 5.1, the sum of several nonconvex functions for Subcases 5.2 and 5.3, and different complicated constraints such as POZ, ramp rate limits, active power spinning reserve, and power losses for Subcase 5.3. Subcase 5.1 is the simplest among the three subcases considering four load values, 2400, 2500, 2600, and 2700 MW. Table 11 reports the comparison of the best cost and the value of from methods. It can see that the proposed method can yield optimal solutions as good as others excluding FA with the worst results and the proposed method’s is also the lowest value excluding comparison to ALHN and EALHN, which are not metaheuristic methods and face to the difficulties of applying for nonconvex fuel cost function. The comparison can imply that the proposed method is very promising for dealing with ten-unit system considering multifuel options where nonconvex fuel cost function or valve-point loading effects are not taken into account.

MethodPD (MW)

ALHN [8]481.723526.239574.381623.809
EALHN [9]481.723526.239574.381623.809
DE [11]481.723526.239574.381623.80912,000
HDP-ICDE [12]481.7226526.2388574.3808623.8098000
HRCGA [21]481.7226526.2388574.3808623.809230,000
AIS [34]481.723526.24574.381623.80924000
CSA [15]623.80926000
CFA [51]623.8339156,000

Comparisons for Subcases 5.2 reported in Table 12 can see that KHA methods in [23] are the best ones with the surprised lowest BCI compared to the proposed method; however, the validation confirmation of solutions cannot be done due to the lack of optimal solution report in [23]. Compared to popular methods such as PSO, GA, TSA (Tabu search algorithm), and FA, the proposed method is better in terms of lower cost and lower . On the contrary, the proposed method is less effective than remaining methods in [13, 15, 16, 19]. As we set to higher values such as 3000, 7000, and 10,000, the results obtained are also worse than these methods. Thus, we conclude that the proposed method is not much effective for the subcase considering valve-point loading effects.

MethodBest cost ($)Mean cost ($)Worst cost ($)BCI (%)

ICSA [16]623.8684623.9495626.366612,000−0.0013
CCDE [13]623.8288623.8574623.89047000−0.0077
CSA [15]623.8684623.9495626.366610,000−0.0013
DSPSO-TSA [19]623.8375623.8625623.90013000−0.0063
PSO [19]624.3045624.5054625.925230000.0685
GA [19]624.505624.7419624.816930000.1006
TSA [19]624.3078635.0623624.828530000.0690
KHA1 [23]611.3276613.0895614.829310,000−2.0528
KHA2 [23]609.0768610.3271611.210510,000−2.4299
KHA3 [23]607.5437608.1164.608.543110,000−2.6884
KHA4 [23]605.7582605.8043605.942610,000−2.9911

The same conclusion is also seen for Subcase 5.3 reported in Table 13 since the proposed method can improve the results 0.026% compared to RCGA and 7.284% compared to FA, but the results are 0.022% worse than that of NRCGA. It is clear that the proposed method cannot be better than further improved method like NRCGA. Once again, the proposed method is not effective for system considering VPLE.

MethodBest cost ($)Mean cost ($)Worst cost ($)BCI (%)

RCGA [20]624.6605625.9201628.925330000.026
NRCGA [20]624.355624.5792624.75413000−0.022

5.2.6. Comparisons of Test Case 6

An IEEE-118 bus system is employed for the test case. The system is similar to a real system with the presence of nearly all electric components such as transformers, capacitors, transmission lines with impedance and admittance, active and reactive power of generators, voltage of load buses and generator buses, the transmission capacity of conductors, and tap setting of transformers. The system has 64 load buses, 186 transmission lines, 54 thermal generating units with quadratic fuel cost function, 14 capacitors, and 9 transformers. The results of minimum cost, mean cost, worst cost and , and the best cost improvement of the proposed method compared to others are reported in Table 14. It can be seen from BCI that MSA [43] is the best method with the best cost improvement of 0.93% compared to the proposed method. On the other hand, the proposed method can improve the best cost from 1.7% (compared to COA) to 10.08% (compared to PSO); meanwhile, its is equal to of PSO and PG-CF-PSO, and much lower than that of all remaining methods. In summary, the proposed method is a promising method for the system with 54 units, quadratic fuel cost function, and all constraints in transmission power networks.

MethodBest cost ($)Mean cost ($)Worst cost ($)BCI (%)

PSO [17]145,520.01158,596.173184,686.82510,00010.08
PG-CF-PSO [17]139,604.13152,204.261170,022.97310,0006.27
COA [39]133,110.43153,110.432138,260.4031.70
ICBOA [38]135,121.57450,0003.16
CBO [38]135,073.00450,0003.13
ECBO [38]135,172.27450,0003.20
DE [38]142,751.118225,0008.34
ABC [38]135,145.189450,0003.18
BBOA [38]135,272.196450,0003.27
MSA [43]129,640.719−0.93

The optimal solutions obtained by the proposed method for Cases 1, 2, 4, and 6 and Subcase 5.3 are given in Tables 1517.

Case 1Case 2Subcase 5.3
Unit (MW)Unit (MW)Unit (MW)


Unit (MW)Unit (MW)Unit (MW)Unit (MW)



(MW)0.8976 (MW)34.7545 (pu)1.0133
(MW)1.7628 (MW)34.628 (pu)1.0264
(MW)0.1324 (MW)19.5482 (pu)1.0033
(MW)12.2424 (MW)56.655 (pu)1.0053
(MW)395.304 (MW)6.0809 (pu)0.993
(MW)76.2481 (MW)29.5586 (pu)0.9881
(MW)10.092 (MW)15.4431 (pu)0.9933
(MW)40.2488 (MW)0.1302 (pu)0.9868
(MW)29.3122 (MW)30.0209 (pu)1.0214
(MW)0.6953 (pu)1.0489 (pu)1.0158
(MW)188.9729 (pu)1.0345 (pu)1.0278
(MW)269.9192 (pu)1.0035 (pu)1.0094
(MW)21.4182 (pu)0.9753 (pu)1.0106
(MW)7.1482 (pu)1.0454 (pu)1.0287
(MW)6.0145 (pu)1.0275 (pu)0.996
(MW)2.3381 (pu)1.06 (pu)1.0191
(MW)0.892 (pu)1.0618 (pu)0.9841
(MW)36.72 (pu)1.0581 (pu)1.0346
(MW)28.3862 (pu)1.0388 (pu)1.0382
(MW)20.5522 (pu)1.0586 (MVAr)−14.3
(MW)193.6641 (pu)1.0993 (MVAr)9.8
(MW)50.635 (pu)1.0322 (MVAr)−15.2
(MW)66.2727 (pu)1.035 (MVAr)3.9
(MW)70.6712 (pu)1.0407 (MVAr)5.5
(MW)131.0352 (pu)1.0137 (MVAr)4.7
(MW)151.3198 (pu)1.0087 (MVAr)8.5
(MW)0.9109 (pu)1.0247 (MVAr)11.3
(MW)341.1965 (pu)1.0525 (MVAr)11.8
(MW)295.7935 (pu)1.028 (MVAr)0
(MW)22.3211 (pu)1.0495 (MVAr)6.3
(MW)0.0596 (pu)1.0207 (MVAr)14.8
(MW)0.2632 (pu)1.0253 (MVAr)5.3
(MW)22.6539 (pu)1.0209 (MVAr)2.6
(MW)21.9574 (pu)1.0345 (pu)0.98