An economic/emission dispatch based on a new multi-objective artificial bee colony optimization algorithm and NSGA-II

Abstract

The conventional energy resources have limited reserves and their utilization is adversely affecting the environment. Hence, it is necessary to generate electricity with the least cost and emission. The studies in the past reveal that the combined economic emission dispatch (CEED) problem has been solved by evolutionary and swarm intelligence-based optimization algorithms. However, the methodology to identify the best compromising solution for the CEED problem has not been studied much in the literature. In this paper, a multi-objective optimization algorithm, which reduces fuel cost of power generation as well as emission simultaneously, is presented. The algorithm is a combination of an artificial bee colony algorithm (ABC) and a non-dominated sorting genetic algorithm (NSGA-II) with a new constraint handling feature. To validate the effectiveness of the proposed algorithm it is applied to three benchmark systems, commonly used to study the effectiveness of optimization algorithms. Moreover, the best compromising solution, obtained by the proposed algorithm, is identified by using sixteen multi-attribute decision-making (MADM) methods. The non-dominated solutions reported in the past literature, for different test systems, are also analysed using MADM methods. It has been shown that the proposed algorithm gives better results without violating the constraints and with the minimum number of iterations.

Appendix 1

1.1 Analysis of the time consumed by the algorithm depending on size of population

To test the effect of population size on the time consumed by the proposed algorithm, the algorithm was run multiple times on 2 separate test systems for both ED and CEED problem. For every run, the size of population (foodnumber) was increased. 12 observations were taken which can be seen in the graphs below as blue lines with asterisks (Actual Values). Then, the graph of approximate equation was plot on the same figure, which can be seen as dashed orange lines (Approximation).

This experiment was conducted on a personal laptop with intel core i3-7020U processor and 8 GB of ram.

For 6-unit CEED problem, calculated approximated equation is:

$${\text{T}} = 1.9747*10^{ - 5} *{\text{N}}^{2} + 0.0604*{\text{N}} + 0.2387$$

In each run, the max_iter parameter was set to 300 iterations. The coefficient of N² is small and its effect can be neglected for population size of less than 1000. Usually in optimization problems, the population size is less than or equal to 1000, so increase in time is almost linear. This behavior can be observed in Fig. 12, which shows change in time consumed, with respect to varying population (foodnumber) size. 

Fig. 12

Time v/s Foodnumber graph for CEED problem on 6-unit test system

For 6-unit economic load dispatch (ELD) problem (considering a single objective of minimization of fuel cost), approximated equation is:

$$T = 1.4611*10^{ - 6} *{\text{N}}^{2} + 0.0101*{\text{N}} + 0.2051$$

In each run, the max_iter parameter was set to 200 iterations, which was sufficient to give best results. Here the coefficient of N² is even smaller than that of CEED problem. In the Fig. 13, the growth in time consumed for increase in foodnumber is not strictly linear (at foodnumber 600,700,800 etc.), rather its oscillating. The overall approximation line is almost linear.

Fig. 13

Time v/s Foodnumber graph for ELD problem on 6-unit test system

The same experiment was repeated on a larger test system, to check whether the linear growth in time can also be seen here. Keeping rest of the parameters same, the algorithm was run 12 times on 40-unit system for both CEED and ELD problems.

For 40-unit CEED problem, approximated equation is:

$$T = 8.9034*10^{ - 5} *{\text{N}}^{2} + 0.0788*{\text{N}}{-}0.9146$$

In each run, the max_iter parameter was set to 600 iterations. The coefficient of N² is higher than those in case of 6-unit test system. From the Fig. 14, it can be seen that the approximated line is not linear, it follows a parabolic path.

Fig. 14

Time v/s Foodnumber graph for CEED problem on 40-unit test system

For 40-unit ELD problem, approximated equation is:

$${\text{T}} = - 1.2202*10^{ - 6} *{\text{N}}^{2} + 0.0418*{\text{N}} + 0.971$$

In each run, the max_iter parameter was set to 600 iterations, which was sufficient to give best results. In the Fig. 15, it can be seen that the growth in time taken is almost linear, for both actual values and approximated lines.

After looking at these real time performances of the algorithm for given test conditions, it can be concluded that:

• 1. Although theoretically the complexity of proposed approach is O(DN²), the growth in time consumed is actually pretty linear, provided that the size of population (foodnumber) is less than 1000.

• 2. The algorithm is scalable in term of both population size (to get a greater number of solutions to choose from) and number of units in test systems (to apply for larger practical systems). The algorithm gives better results in any situation while consuming

The foodnumber was set to 50 for the actual test systems, since it gave good results and in acceptable time.

Fig. 15

Time v/s Foodnumber graph for ELD problem on 40-unit test system

1.2 Analysis of time consumed by repair algorithm

The while loop inside repair algorithm operates on a single solution until its constraint violation (error) becomes less than or equal to tolerance value. Since every solution is unique, each of them will take different amounts of iteration of the while loop. To determine time consumed in such case, the algorithm was run on CEED problem on different test systems and number of iterations took by repair operator were noted. Each iteration takes constant amount of time to run, thus the overall time taken by repair operation is directly proportional to the number of iterations.

Table

Table 15 Results of the experiment to analyse time complexity of repair algorithm

15 shows in brief the results of above experiment.

The above procedure was repeated 5 times on each test system and similar results were obtained as shown in the Table 15.

In case of 6-unit system, in which loss matrix were considered, there is only 1 solution which took 4490 iterations to complete the repair process, which is like an outlier for this experiment. 95% percentile is at 50 iterations, which is good. Average number of iterations is less than 30.

In case of 13-unit system, in which loss matrix is not considered, only 1 solution took 97 iterations to complete the repair process. Only 3 solutions took greater than 80 iterations to complete the process. Average number of iterations is less than 2, and standard deviation is less than 3. 95% of the solutions took less than 3 iterations to complete the repair process.

In case of 40-unit system, considering the transmission losses, the iteration taken for repair process are even less than those in 6-unit test system. Here an outlier solution is present which took 99 iterations to complete the repair process, which is still acceptable. Only 2 solutions took more than 80 iterations to complete the repair process. Average number of iterations is less than 4, with standard deviation less than 3. 95% of the solutions required less than 7 iterations to complete the repair process.

From the experiment it is clear that due to randomness present in the repair algorithm, the iterations performed on a single solution and thus the time consumed by repair algorithm is unpredictable beforehand. But, considering the average case performance over all runs and overall time took by whole algorithm (sub 5 seconds in every test case as discussed in section 5), repair algorithm can be used for large scale application without much effect on the performance.

1.3 Reproducibility of proposed algorithm

The Tables 16, 17 and 18 give the overall analysis of the results produced by proposed algorithm for 6 unit, 13 unit and 40 unit test systems respectively. It can be observed that the results produced in each run are almost equal. Also, due to repair algorithm, none of the solution in any test case or any test run have a constraint violation error. This shows that the results by proposed algorithm are reproducible

Table 16 Results for Best Compromising solutions for 6-unit system for a load of 1200 MW by MOABC-NSGA-II method for multiple test runs

Table 17 Results for Best Compromising solutions for 13-unit system for a load of 1800 MW by MOABC-NSGA-II method multiple test runs

Table 18 Results for Best Compromising solutions for 40-unit system for a load of 10,500 MW by MOABC-NSGA-II method multiple test runs

1.4 Demonstration of working of algorithm

Here, a multi-objective constrained benchmark function, Binh and Korn function is solved using the proposed algorithm. In this problem, 2 conflicting objectives functions f1 and f2 are to be minimized, subject to 2 inequality constraints c1 and c2. There are 2 decision variables, x and y. The problem is formulated as:

Minimize:

$$\begin{aligned} & f_{1} (x,y) = 4x^{2} + 4y^{2} ; \\ & f_{2} (x,y) = (x - 5)^{2} + (y - 5)^{2} ; \\ \end{aligned}$$

Subject to constraints:

$$\begin{aligned} & c_{1} (x,y) = (x - 5)^{2} + y^{2} \le 25; \\ & c_{2} (x,y) = (x - 8)^{2} + (y + 3)^{2} \ge 7.7; \\ \end{aligned}$$

Search domain:

$$\begin{aligned} & 0 \le x \le 5 \\ & 0 \le y \le 3 \\ \end{aligned}$$

To keep things simple, the population size is declared as 3. Rest parameters are same. Here only 1 iteration of the algorithm is explained. The repair algorithm proposed in this research is only for economic/emission dispatch problems and thus will not work for other benchmark problems. This problem is solved simultaneously on MATLAB, so all the matrices shown here are the actual values given by the program.

In the initialization phase, population X is initiated using Eq. (7) as:

$${\text{X}} = \left[ {\begin{array}{*{20}c} {1.643661} & {0.492924} \\ {4.660542} & {0.789429} \\ {1.44897} & {1.639382} \\ \end{array} } \right]$$

Each row represents a different solution / individual in the population. First column represents first decision variable (x), and 2^nd column represents 2^nd decision variable (y). All values are initialized within the search domain for respective decision variable. Trial counter for each solution is initiated to zero.

$${\text{trial}}\_{\text{counter}} = \left[ {\begin{array}{*{20}c} 0 \\ 0 \\ 0 \\ \end{array} } \right]$$

Constraint violation for all solutions is calculated, which are as follows:

$${\text{c}}1 = \left[ {\begin{array}{*{20}c} {11.50799} \\ {0.73843} \\ {15.32633} \\ \end{array} } \right]\quad {\text{c}}2 = \left[ {\begin{array}{*{20}c} {52.60356} \\ {25.51175} \\ {64.49324} \\ \end{array} } \right]$$

Each row represents the c1 and c2 values of respective solution. For all solutions, c1 is less than 25 and c2 is greater than 7.7. So, all constraints are satisfied. Thus, error matrix has all values zero.

$${\text{error}} = \left[ {\begin{array}{*{20}c} 0 & 0 \\ 0 & 0 \\ 0 & 0 \\ \end{array} } \right]$$

Now, objective function values for all solutions are calculated and stored in matrix Z as:

$$Z = \left[ {\begin{array}{*{20}c} {11.77838} &\quad {31.57875} \\ {89.37539} &\quad {17.84414} \\ {19.1012} &\quad {23.93251} \\ \end{array} } \right]$$

First column stores value of objective f1 and 2nd column stores value of objective f2.

Now main loop of the algorithm starts with employed bee phase.

Consider the first solution in the population. The decision variable to change is selected randomly by MATLAB to be 2. Neighbouring solution is also selected randomly by MATLAB, which is solution at index 3.

Using Eq. (8), algorithm creates a new temporary solution (xx) as

$${\text{xx}} = \left[ {\begin{array}{*{20}c} {1.507249} & {0.492924} \\ \end{array} } \right]$$

cosntraint violation error value and objective function cost value of the new solution are given as:

$${\text{xx}}\_{\text{error}} = \left[ {\begin{array}{*{20}c} 0 & 0 \\ \end{array} } \right]$$

$${\text{xx}}\_{\text{cost}} = \left[ {\begin{array}{*{20}c} {10.05909} & {32.51305} \\ \end{array} } \right]$$

Now to check if new solution dominates old solution, first error values are compared. In this case, for both variables it is zero, so now objective function values are compared. Comparing the first row of Z matrix and xx_cost for each objective, it can be observed that neither solution clearly dominates the other. Since the new solution is not dominating the old solution, new solution is discarded. Trial counter of old solution is increased by 1.

Same procedure is repeated for solutions 2 and 3 completing the employed bee phase. No improvement was done in either solution 2 or 3. Trial counter of all solution is now equal to 1.

Now, fitness value and probability of each solution is calculated. For which non-dominated sorting [159] is performed on the whole population. First, feasible and non-feasible solutions are separated. Here all solutions are feasible solutions. For feasible solutions, dominance is decided based on the value of objective function. In current iteration, all solutions are non-dominant to each other. Thus, all solutions are given rank value of 1. Rank matrix can be represented as:

$${\text{rank}} = \left[ {\begin{array}{*{20}c} 1 \\ 1 \\ 1 \\ \end{array} } \right]$$

After this, crowding distance is assigned to each solution. For solutions at both ends each rank (front), crowding distance is infinity (Inf). Crowding distance for rest solutions is calculated using the formula discussed in [159]. In this case, crowding distance is assigned as follow:

$${\text{crowding}}\_{\text{distance}} = \left[ {\begin{array}{*{20}c} {{\text{Inf}}} \\ {{\text{Inf}}} \\ 0 \\ \end{array} } \right]$$

After this, the whole population and its related matrices are sorted based on values of rank and crowding distance. First whole population is sorted in descending order using crowding distance, then it is sorted in ascending order using rank. This sorting satisfied the criteria discussed in section 4.4.

For this example, the sorting doesn’t affect the order of elements.

After this, fitness value is assigned to the solutions based on their position in the population. Since the population is sorted, the number of solutions dominated by any solution is equal to the number of solutions in the population present below it. Fitness value and then probability of selection are calculated using Eq. (10) and Eq. (11) respectively.

$${\text{fitness}}\_{\text{value}} = \left[ {\begin{array}{*{20}c} {0.6667} \\ {0.3333} \\ 0 \\ \end{array} } \right]$$

$${\text{probabilty}}\_{\text{value}} = \left[ {\begin{array}{*{20}c} {0.6667} \\ {0.3333} \\ 0 \\ \end{array} } \right]$$

The onlooker bee phase is like employed bee phase. The only difference is in the selection of neighbour solution. In onlooker bee phase, a solution with higher probability value has better chance of getting selected as a neighbour. Also, neighbour should not be the same as current element. Roulette wheel selection can be used to select the neighbour based on their probability.

After selection of neighbour solution, just like employed bee phase, one of the decision variables is randomly selected. Then using Eq. (9), a new solution xx is created. If this new solution dominates old solution, substitution takes place, else trial counter of old solution is increased by 1.

In current problem, no new dominating solution was created, so the population and rest of the matrices remained unchanged. Trial counter for all solutions is now 2.

The scout bee phase is initialized only when trial_counter for any solution in the population is greater than trail_limit. The old solution is discarded, and new solution is initialized randomly using Eq. (7). Its objective function values and constraint violation error values are calculated. All these values are replaced in the original matrices in position of old element.

In our problem, no solution entered scout bee phase, as trial_limit for crossover was set to 60.

After this, crossover is done to generate new solutions using already present solutions. Here, probability of crossover (pc) was set to 0.9. So total number of crossovers done is given by

$$Number\,of\,crossovers = round\left( {\frac{food\,number*pc}{2}} \right) = round\left( {\frac{3*0.9}{2}} \right) = 1$$

Each crossover produces 2 offspring, so total number of new solutions formed = 2.

For each crossover, 2 parents from the current population are selected using tournament_selection method. After that SBX crossover operation [159] is performed, which creates 2 new offsprings (Xcross). In current problem, the values got are as follows:

$${\text{X}}\_{\text{cross}} = \left[ {\begin{array}{*{20}c} {1.591839} & {0.509774} \\ {4.712364} & {0.772579} \\ \end{array} } \right]$$

Each row represents a newly formed solution. Columns represent x and y decision variables respectively.

In mutation, a single parent is selected using tournament selection method. Then, polynomial mutation [159] is performed to create a new offspring (X_mutate). Each mutation operation creates a single new offspring.

For mutation, probability of mutation (mu) was set to 0.1. So total number of mutations done is given by

$$Number\,of\,mutations = round\left( {foodnumber * mu} \right) = round\left( {3*0.1} \right) = 0$$

Due to small size of population no mutation operation was performed.

All these newly from solutions are added to the population, at the end. So now the population matrix X is

$${\text{X}} = \left[ {\begin{array}{*{20}c} {1.643661} &\quad {0.492924} \\ {4.660542} &\quad {0.789429} \\ {1.44897} &\quad {1.639382} \\ {1.591839} &\quad {0.509774} \\ {4.712364} &\quad {0.772579} \\ \end{array} } \right]$$

Similarly, objective function cost values and constraint violation error values for new solutions are calculated. These values are appended in the respective matrices. The trial_counter for new solution is set to zero, which is appended to the previous matrix.

$$Z \, = \left[ {\begin{array}{*{20}c} {11.77838} &\quad {31.57875} \\ {89.37539} &\quad {17.84414} \\ {19.1012} &\quad {23.93251} \\ {11.17528} &\quad {31.7777} \\ {91.213} &\quad {17.95382} \\ \end{array} } \right]\quad {\text{error}} = \left[ {\begin{array}{*{20}c} 0 & 0 \\ 0 & 0 \\ 0 & 0 \\ 0 & 0 \\ 0 & 0 \\ \end{array} } \right]\quad {\text{trial}}\_{\text{counter }} = \left[ {\begin{array}{*{20}c} 0 \\ 0 \\ 0 \\ 0 \\ 0 \\ \end{array} } \right]$$

After this, a final non-dominated sorting operation is performed on the whole operation. Since all solutions are feasible, ranks are assigned according to the objective function values. Then crowding distance is assigned. The result of these operation is as follow:

$${\text{rank}} = \left[ {\begin{array}{*{20}c} 1 \\ 1 \\ 1 \\ 1 \\ 2 \\ \end{array} } \right]\quad {\text{crowding}}\_{\text{distance}} = \left[ {\begin{array}{*{20}c} {0.642406} \\ {Inf} \\ {0.642406} \\ {Inf} \\ {Inf} \\ \end{array} } \right]$$

Then, the whole population and its related matrices are sorted based on values of rank and crowding distance according to the criteria discussed in section 4.4. After sorting, the matrices look like this:

$${\text{X}} = \left[ {\begin{array}{*{20}c} {4.660542} &\quad {0.789429} \\ {1.591839} &\quad {0.509774} \\ {1.643661} &\quad {0.492924} \\ {1.444897} &\quad {1.639382} \\ {4.712364} &\quad {0.772579} \\ \end{array} } \right]\quad {\text{Z}} = \left[ {\begin{array}{*{20}c} {89.37539} &\quad {17.84414} \\ {11.17528} &\quad {31.7777} \\ {11.77838} &\quad {31.57875} \\ {19.1012} &\quad {23.93251} \\ {91.213} &\quad {17.95382} \\ \end{array} } \right]$$

$${\text{trial}}\_{\text{counter}} = \left[ {\begin{array}{*{20}c} 2 \\ 0 \\ 2 \\ 2 \\ 0 \\ \end{array} } \right]\quad {\text{crowding}}\_{\text{distance}} = \left[ {\begin{array}{*{20}c} {Inf} \\ {Inf} \\ {0.642406} \\ {0.642406} \\ {Inf} \\ \end{array} } \right]\quad {\text{ranks}} = \left[ {\begin{array}{*{20}c} 1 \\ 1 \\ 1 \\ 1 \\ 2 \\ \end{array} } \right]$$

Finally, extra solutions from the population are removed using truncation operation. Here, the two elements at the bottom of the matrix are removed because they are the weakest solutions. Similarly, last 2 rows are removed from the rest of the data matrices. This keeps the population size equal to the food number.

After completing one iteration, the matrices are as follow:

$${\text{X}} = \left[ {\begin{array}{*{20}c} {4.660542} &\quad {0.789429} \\ {1.591839} &\quad {0.509774} \\ {1.643661} &\quad {0.492924} \\ \end{array} } \right]\quad {\text{Z}} = \left[ {\begin{array}{*{20}c} {89.37539} &\quad {17.84414} \\ {11.17528} &\quad {31.7777} \\ {11.77838} &\quad {31.57875} \\ \end{array} } \right]$$

$${\text{trial}}\_{\text{counter}} = \left[ {\begin{array}{*{20}c} 2 \\ 0 \\ 2 \\ \end{array} } \right]\quad {\text{crowding}}\_{\text{distance}} = \left[ {\begin{array}{*{20}c} {Inf} \\ {Inf} \\ {0.642406} \\ \end{array} } \right]\quad {\text{ranks}} = \left[ {\begin{array}{*{20}c} 1 \\ 1 \\ 1 \\ \end{array} } \right]$$

This whole procedure repeats for multiple iterations, until maximum number of iterations is reached. Since the population size was too low in this case, the results obtained were not good enough.

When the same problem is solved with population size of 50 and for 100 iterations, the algorithm took 1.3925 seconds and the pareto front obtained was as follow (Fig.

Fig. 16

Pareto Optimal front for Binh and Korn function

16).