# Artificial Bee Colony Algorithm Based onK-Means Clustering for Multiobjective Optimal Power Flow Problem

Authors: Liling Sun; Jingtao Hu; Hanning Chen
Journal: Mathematical Problems in Engineering
 (2015)

---

## Abstract

An improved multiobjective ABC algorithm based onK-means clustering, called CMOABC, is proposed. To fasten the convergence rate of the canonical MOABC, the way of information communication in the employed bees’ phase is modified. For keeping the population diversity, the multiswarm technology based on K-means clustering is employed to decompose the population into many clusters. Due to each subcomponent evolving separately, after every specific iteration, the population will be reclustered to facilitate information exchange among different clusters. Application of the new CMOABC on several multiobjective benchmark functions shows a marked improvement in performance over the fast nondominated sorting genetic algorithm (NSGA-II), the multiobjective particle swarm optimizer (MOPSO), and the multiobjective ABC (MOABC). Finally, the CMOABC is applied to solve the real-world optimal power flow (OPF) problem that considers the cost, loss, and emission impacts as the objective functions. The 30-bus IEEE test system is presented to illustrate the application of the proposed algorithm. The simulation results demonstrate that, compared to NSGA-II, MOPSO, and MOABC, the proposed CMOABC is superior for solving OPF problem, in terms of optimization accuracy.


---

## Body

## 1. Introduction

Real-life optimization problems usually have several conflicting objectives. The solution of such problems is very difficult because their objectives tend to be in conflict with each other. So a decision maker (DM) needs to determine an optimal trade-off among them. These problems can be modeled asmultiobjective optimization problems (MOPs). Generally, the optimal solution of a multiobjective optimization problem is a set of optimal solutions (largely known as Pareto-optimal solutions). A Pareto-optimal solution to an MOP is a candidate for the optimal trade-off but it is not possible to find a better solution than others from the Pareto-optimal solutions, because one cannot be better than others without any further information. Therefore, it is necessary to find as many Pareto-optimal solutions as possible [1].Swarm intelligence (SI) is an innovative artificial intelligence technique for solving complex multiobjective optimization problems, such as Pareto-archived evolution strategy (PAES) [2], Pareto envelope-based selection algorithm- (PESA-) II [3], nondominated sorting genetic algorithm II (NSGA-II) [4], strength Pareto evolutionary algorithm (SPEA2) [5], indicator-based evolutionary algorithm (IBEA) [6], multiobjective particle swarm optimization (MOPSO) [7], multiobjective evolutionary algorithm based on decomposition (MOEA/D) [8], and two-lbest multiobjective particle swarm optimization (2LB-MOPSO) [9]. As the widely adopted intelligent search method, the artificial bee colony (ABC) algorithm is a powerful search technique that drew inspiration from the biological foraging behaviors observed in bee colony [10, 11]. Due to its simplicity and efficiency, ABC has become a high competitor in solving the MO problems. Many researchers have presented several existing multiobjective ABC (MOABC) algorithms [12, 13]. However, these proposed algorithms still suffer from low convergence rate and lack the diversity of swarm.To conquer the weakness of initial MOABC, an improved multiobjective ABC algorithm based onK-means clustering, named CMOABC, is proposed. In the CMOABC, for accelerating the convergence rate of canonical MOABC, the employed bees’ phase is modified to adjust the strength of information communication. Additionally, the population is partitioned into several subpopulations based on K-means clustering. Information communication between the subpopulations depends on reclustering the population after each specific iteration. To further enhance the population diversity, a number of individuals with worse performance regenerate in the reclustering process. With the test against several multiobjective benchmark functions, the results demonstrate that the proposed CMOABC outperforms the MOABC, MOPSO, and NSGA-II on convergence metric and diversity metric.In the real-world optimization problems, optimal power flow (OPF) is a classical multiobjective problem. Traditionally, the basic objective of optimal power flow (OPF) problem is to schedule the committed generating units to meet the system load demand at minimum operating cost while satisfying the various system equality and inequality constraints [14]. But the passage of clean air act amendments in 1990 forced the utilities to reduce the emission from fossil fuel fired thermal station [15]. Therefore, in addition to fuel cost, emission must also be considered as an objective. On the whole, OPF problem is a nonlinear, constrained optimization problem where many competing objectives are present.In the early research, many traditional optimization techniques such as linear programming (LP) [16], nonlinear programming (NLP) [17, 18], quadratic programming (QP) [19], Newton method, and interior point method (IPM) [20–22] have already been employed to solve the OPF problem. However, these methods suffer from severe limitations in handling the OPF problem, including discrete-continuous functions, and complex constraints [20]. In order to overcome the limitations of classical optimization techniques, recently, a wide variety of the multiobjective heuristic evolutionary algorithms have been proposed to solve OPF problem, mainly including genetic algorithm (GA) [23, 24], enhanced genetic algorithm (EGA) [25], evolutionary programming (EP) [26, 27], multihive multiobjective bee algorithm (M2OBA) [28], modified shuffle frog leaping algorithm (MSFLA) [29], differential evolution (DE) [30, 31], harmony search (HS) algorithm [32], biogeography based optimization (BBO) [33, 34], multiobjective evolutionary algorithm (MOEA) [35], multiobjective particle swarm optimization [36, 37], and fuzzy clustering-based particle swarm optimization (FCPSO) [38]. However, all the mentioned mathematical techniques have some drawbacks such as being trapped in local optima and they are only suitable for considering a specific objective function in the OPF problem. Due to its outstanding performance on the test against benchmark functions, CMOABC is suitable for solving the multiobjective problems. Therefore, CMOABC is utilized to solve OPF problem. Compared with MOABC, MOPSO, and NSGA-II, the proposed CMOABC algorithm can accommodate considerable potential for solving OPF problem.The rest of this paper is organized as follows. Section2 presents optimal power flow (OPF) formulation, including the objective function and the constraints of the OPF problem. The proposed algorithm CMOABC is detailed in Section 3, and also in this section we test the CMOABC on benchmark functions: four ZDT and two DTLZ. Afterward, in Section 4, the CMOABC is applied to solve the OPF problem. Finally, Section 5 outlines the conclusions.


## 2. Optimal Power Flow Problem Formulation

The multiobjective optimal power flow (OPF) is a nonlinear, nonconvex optimization problem. The main goal of OPF is to optimize the settings of control variables in terms of one or more objective functions while satisfying several equality and inequality constraints. In general, the OPF problem can be mathematically formulated as follows:(1)MinimizeFx,u=f1x,u,…,fmx,u(2)Subjected  togx,u=0,(3)hx,u≤0,where F(x,u) is the objective function to be optimized, g(x,u) is the equality constraints representing nonlinear power flow equations, and h(x,u) is the system operating constraints. Here, u is the vector of independent control variables including(1)
generator active power outputPG except at slack bus PG1;(2)
generator bus voltageVG;(3)
transformer tap settingT;(4)
shunt VAR compensationQC.Hence,u can be expressed as(4)uT=PG2⋯PGNg,VG1⋯VGNg,T1⋯TNt,QC1⋯QCNC,where Ng, Ni, and NC denote the number of generating units, number of regulating transformers, and number of shunt compensators, respectively. Generators active powers (except slack bus PG1) and generators bus voltages are continuous variables, while the tap settings of the tap changing transformers and VAR injections of the shunt capacitors are discrete variables.Andx is the vector of dependent variables including(1)
slack bus generated active powerPG1;(2)
load (PQ) bus voltageVL;(3)
generator reactive power outputQG;(4)
transmission line loading (line flow)SL.Hence,x can be expressed as(5)xT=PG1,VL1⋯VLNpq,QG1⋯QGNg,SL1⋯SLNl,where Npq is the number of PQ buses and Nl is the total number of transmission lines.


### 2.1. Objective Function

In this paper, the OPF problem is to minimize three competing objective functions, namely, the total fuel cost, the total emission cost minimization, and total power losses, while satisfying several equality and inequality constraints. Generally the problem is formulated as follows.


#### 2.1.1. Minimization of Total Fuel Cost

This objective function is to minimize the total fuel costfcost of the system. The fuel cost curves of the thermal generators are modeled as quadratic cost curves and can be represented as follows:(6)fcost=∑i=1NgfiaiPGi2+biPGi+ci,where ai, bi, and ci are the fuel cost coefficients of the ith generator and PGi is real power output of the ith generator.


#### 2.1.2. Minimization of Total Power Losses

The power flow solution gives all bus voltage magnitudes and angles. Then, the total MW active power loss in a transmission network can be described as follows:(7)flost=∑k=1NlgkVi2+Vj2-2ViVjcos⁡⁡δi-δj,where Nl is the number of transmission lines, Vi and Vj are the voltage magnitudes at the ith bus and jth bus, respectively, and δi and δj are the voltage angles at the ith bus and the jth bus, respectively.


#### 2.1.3. Total Emission Cost Minimization

In this paper, two important types of emission gasses, namely, sulphur oxidesSOx and nitrogen oxides NOx, are taken as the pollutant gasses. The emission gasses generated by each generating unit may be approximated by a combination of a quadratic and an exponential function of the generator active power output. Here, the total emission cost is defined as follows [26]:(8)femission=∑i=1Ngαi+βiPGi+γiPGi2,where femission is the total emission cost (ton/h) and αi, βi, and γi are the emission coefficients of the ith unit.


### 2.2. Equality Constraints

The equality constraintsg(x,u) represented by (2) are typical power flow equations which are defined as follows:(9)PGi-PDi-Vi∑j=1NVjGijcos⁡⁡δi-δj+Bijsin⁡δi-δj=0i=1,2,…,N,QGi-QDi-Vi∑j=1NVjGijcos⁡⁡δi-δj-Bijsin⁡δi-δj=0i=1,2,…,N,where N is the number of buses in the system. PGi and QGi are the active and reactive power generators at the ith bus; PDi and QDi are the active and reactive power loads at the ith bus, respectively; Gij and Bij are the real and imaginary parts of the ijth element of the bus admittance matrix, respectively; Vi and Vj are the voltage magnitudes at the ith bus and the jth bus, respectively; δi and δj are the voltage angles at the ith bus and the jth bus, respectively.The equality constraints of power balance equations shown in (9) are forced by unconstrained Newton-Raphson based power flow calculations; therefore, there is no need to integrate them into the objective function.


### 2.3. Inequality Constraints

The inequality constraintsh(x,u) represented by (3) are the power system operating limits including the following.(1) Generator Constraints. Generator active power PG, generator reactive power QG, and generator voltage magnitude VG are restricted by their lower and upper limits:(10)PGi,min⁡≤PGi≤PGi,max⁡i=1,…,Ng,QGi,min⁡≤QGi≤QGi,max⁡i=1,…,Ng,VGi,min⁡≤VGi≤VGi,max⁡i=1,…,Ng.(2) Transformer Constraints. Transformer taps have minimum and maximum setting limits:(11)TGi,min⁡≤Ti≤Ti,max⁡i=1,…,Nt.(3) Switchable VAR Sources. The switchable VAR sources have restrictions as follows:(12)QCi,min⁡≤QCi≤QCi,max⁡i=1,…,Nc.(4) Security Constraints. These include the limits on the load bus voltage magnitudes and transmission line flows limits:(13)VLi,min⁡≤VLi≤VLi,max⁡i=1,…,Npq,SLi≤SLi,max⁡i=1,…,Nl.In this work, the penalty factor method is utilized for handling the inequality constraints. In this regard, each control vector which violates constraints will be fined by these penalty factors. Therefore, in the next step, this control vector will be deleted automatically.


## 2.1. Objective Function

In this paper, the OPF problem is to minimize three competing objective functions, namely, the total fuel cost, the total emission cost minimization, and total power losses, while satisfying several equality and inequality constraints. Generally the problem is formulated as follows.


### 2.1.1. Minimization of Total Fuel Cost

This objective function is to minimize the total fuel costfcost of the system. The fuel cost curves of the thermal generators are modeled as quadratic cost curves and can be represented as follows:(6)fcost=∑i=1NgfiaiPGi2+biPGi+ci,where ai, bi, and ci are the fuel cost coefficients of the ith generator and PGi is real power output of the ith generator.


### 2.1.2. Minimization of Total Power Losses

The power flow solution gives all bus voltage magnitudes and angles. Then, the total MW active power loss in a transmission network can be described as follows:(7)flost=∑k=1NlgkVi2+Vj2-2ViVjcos⁡⁡δi-δj,where Nl is the number of transmission lines, Vi and Vj are the voltage magnitudes at the ith bus and jth bus, respectively, and δi and δj are the voltage angles at the ith bus and the jth bus, respectively.


### 2.1.3. Total Emission Cost Minimization

In this paper, two important types of emission gasses, namely, sulphur oxidesSOx and nitrogen oxides NOx, are taken as the pollutant gasses. The emission gasses generated by each generating unit may be approximated by a combination of a quadratic and an exponential function of the generator active power output. Here, the total emission cost is defined as follows [26]:(8)femission=∑i=1Ngαi+βiPGi+γiPGi2,where femission is the total emission cost (ton/h) and αi, βi, and γi are the emission coefficients of the ith unit.


## 2.1.1. Minimization of Total Fuel Cost

This objective function is to minimize the total fuel costfcost of the system. The fuel cost curves of the thermal generators are modeled as quadratic cost curves and can be represented as follows:(6)fcost=∑i=1NgfiaiPGi2+biPGi+ci,where ai, bi, and ci are the fuel cost coefficients of the ith generator and PGi is real power output of the ith generator.


## 2.1.2. Minimization of Total Power Losses

The power flow solution gives all bus voltage magnitudes and angles. Then, the total MW active power loss in a transmission network can be described as follows:(7)flost=∑k=1NlgkVi2+Vj2-2ViVjcos⁡⁡δi-δj,where Nl is the number of transmission lines, Vi and Vj are the voltage magnitudes at the ith bus and jth bus, respectively, and δi and δj are the voltage angles at the ith bus and the jth bus, respectively.


## 2.1.3. Total Emission Cost Minimization

In this paper, two important types of emission gasses, namely, sulphur oxidesSOx and nitrogen oxides NOx, are taken as the pollutant gasses. The emission gasses generated by each generating unit may be approximated by a combination of a quadratic and an exponential function of the generator active power output. Here, the total emission cost is defined as follows [26]:(8)femission=∑i=1Ngαi+βiPGi+γiPGi2,where femission is the total emission cost (ton/h) and αi, βi, and γi are the emission coefficients of the ith unit.


## 2.2. Equality Constraints

The equality constraintsg(x,u) represented by (2) are typical power flow equations which are defined as follows:(9)PGi-PDi-Vi∑j=1NVjGijcos⁡⁡δi-δj+Bijsin⁡δi-δj=0i=1,2,…,N,QGi-QDi-Vi∑j=1NVjGijcos⁡⁡δi-δj-Bijsin⁡δi-δj=0i=1,2,…,N,where N is the number of buses in the system. PGi and QGi are the active and reactive power generators at the ith bus; PDi and QDi are the active and reactive power loads at the ith bus, respectively; Gij and Bij are the real and imaginary parts of the ijth element of the bus admittance matrix, respectively; Vi and Vj are the voltage magnitudes at the ith bus and the jth bus, respectively; δi and δj are the voltage angles at the ith bus and the jth bus, respectively.The equality constraints of power balance equations shown in (9) are forced by unconstrained Newton-Raphson based power flow calculations; therefore, there is no need to integrate them into the objective function.


## 2.3. Inequality Constraints

The inequality constraintsh(x,u) represented by (3) are the power system operating limits including the following.(1) Generator Constraints. Generator active power PG, generator reactive power QG, and generator voltage magnitude VG are restricted by their lower and upper limits:(10)PGi,min⁡≤PGi≤PGi,max⁡i=1,…,Ng,QGi,min⁡≤QGi≤QGi,max⁡i=1,…,Ng,VGi,min⁡≤VGi≤VGi,max⁡i=1,…,Ng.(2) Transformer Constraints. Transformer taps have minimum and maximum setting limits:(11)TGi,min⁡≤Ti≤Ti,max⁡i=1,…,Nt.(3) Switchable VAR Sources. The switchable VAR sources have restrictions as follows:(12)QCi,min⁡≤QCi≤QCi,max⁡i=1,…,Nc.(4) Security Constraints. These include the limits on the load bus voltage magnitudes and transmission line flows limits:(13)VLi,min⁡≤VLi≤VLi,max⁡i=1,…,Npq,SLi≤SLi,max⁡i=1,…,Nl.In this work, the penalty factor method is utilized for handling the inequality constraints. In this regard, each control vector which violates constraints will be fined by these penalty factors. Therefore, in the next step, this control vector will be deleted automatically.


## 3. Multiobjective Artificial Bee Colony Algorithm Based onK-Means Clustering

### 3.1. Canonical ABC Algorithm

The artificial bee colony (ABC) algorithm, proposed by Karaboga and Basturk [39] and further developed by Basturk and Akay [11, 40] for real-parameter optimization, which simulates the intelligent foraging behavior of a honeybee swarm, is one of the most recently introduced swarm-based optimization techniques.The entire bee colony contains three groups of bees: employed bees, onlookers, and scouts. Employed bees explore the specific food sources; meanwhile they pass the food information to onlooker bees. The number of employed bees is equal to that of food sources; in other words, each food source owns only one employed bee. Then onlooker bees choose good food sources based on the received information and then further exploit the food near their selected food source. The food source with higher quality would have a larger opportunity to be selected by onlookers. There is a control parameter called “limit” in the canonical ABC algorithm. If a food source is not improved anymore whenlimit is exceeded, it is assumed to be abandoned by its employed bee and the employed bee associated with that food source becomes a scout to search for a new food source randomly. The fundamental mathematic representations are listed as follows [41].Step 1 (initialization phase).
In initialization phase, a group of food sources are generated randomly in the search space using the following equation:(14)xi,j=xjmin⁡+rand⁡0,1xjmax⁡-xjmin⁡,where i=1,2,…,SN, j=1,2,…,D. SN is the number of food sources. D is the number of variables, that is, problem dimension. xjmin⁡ and xjmax⁡ are the lower and upper bounds of the jth variable, respectively.Step 2 (employed bees’ phase).
In the employed bees’ phase, the neighbor food source (candidate solution) can be generated from the old food source of each employed bee in its memory using the following expression:(15)vi,j=xi,j+δxi,j-xk,j,where k is a randomly selected food source and must be different from i; j is randomly chosen indexes; δ is a random number in range [-1,1].Step 3 (onlooker bees’ phase).
In the onlooker bees’ phase, an onlooker bee selects a food source depending on the probability value associated with that food source andPi can be calculated as follows:(16)Pi=fitnessi∑j=1SNfitnessj,where fitnessj is the fitness value of the jth solution.Step 4 (scout bees’ phase).
In scout bees’ phase, if a food source cannot be improved further through a predetermined cycle (called “limit” in ABC), the food source is supposed to be abandoned. The employed bee subsequently becomes a scout. A new food source will be produced randomly in the search space using (14).The employed, onlooker, and scout bees’ phase will recycle until the termination condition is met.


### 3.2. CMOABC for Optimization

This section presents the detailed description of the proposed multiobjective artificial bee colony algorithm based onK-means clustering (named CMOABC). The proposed algorithm makes improvements on two aspects: (1) modifying the method of information communication in employed bees’ phase and (2) partitioning the population using K-means clustering. The number of the clusters will be chosen from the predefined set G. After each specific iteration, when the defined criteria are met, parts of individuals in clusters will be removed and an equal number of new individuals are regenerated randomly. Then the number of clusters will change into a new value of the set G and the population will be reclustered based on K-means clustering.


#### 3.2.1. Modified Employed Bees’ Phase of ABC

To increase the exploitation of ABC algorithm and fasten the convergence rate, the employed bees’ phase is modified. The new food source (candidate solution) is generated using the following way:(17)vij=xij+ϕxij-xkjrand≥Cr,vij=xij+φxbest,j-xijrand<Cr,(18)Cr=iteritermax⁡,where xbest,j is the position of the global best solution and the term xbest,j can drive the new candidate solution to the global best solution, iter is the current iteration, and itermax⁡ is the maximum iteration.The parameterCr in (18) plays an important role in balancing the exploration and exploitation of the candidate solution search. If Cr=0, (18) is identical to (15). With the increase of Cr, the probability of the candidate solution learning to the best solution increases correspondingly. In this way, in the beginning of optimization process, the proposed algorithm operates as the canonical artificial bee colony algorithm, which can well keep the population diversity; in the end of optimization process, the modified ABC algorithm has a considerable improvement on both convergence rate and local search.Furthermore, in (17), the different learning intensities, which are determined by ϕ and φ, are adopted. We adopt the beta probability distribution [42] to tune these two parameters, which is based on the similar strategy of tuning parameters in [43]. The beta distribution is flexible for modeling data that are measured in a continuous scale on a truncated interval in range [0,1], since its density is a versatile way to represent different shapes depending on the values of the two parameters that index the distribution. The probability density (Bd_pd) of the beta distribution in the range [0,1] is calculated as follows:(19)Bd_pd=xα-11-xβ-1Bα,β,where B is the beta function.The mean (E) and variance (Var) of the beta distribution are calculated as follows:(20)EX=αα+β,Var⁡X=αβα+β2α+β+1.The values ofϕ and φ are tuned by MATLAB∗ script (see (21)). An example of ϕ and φ values generated using beta distribution is shown in Figure 1. Consider(21)α=rand,β=0.7∗iteritermax⁡+0.2∗rand,ϕ=betarnda,b,11,α=rand,β=0.6∗iteritermax⁡+0.1∗rand,φ=betarnda,b,11,where itermax⁡ denotes the max iteration and iter is the current iteration.Figure 1
(a)F and (b) CR values for a run generated using beta distribution in the beta distribution.
(a)
(b)


#### 3.2.2. Cluster Setting

TheK-means clustering method is employed to partition the population into subpopulations. The basic concepts of K-means clustering are presented firstly, and then we will give a detailed description of its application in our proposed algorithm.(I) Basic Parameters of
K
-Means. The cluster centers are substituted for center positions of food sources and the formula of computing the centers is shown in (22). If the ith cluster contains ni members and the members are denoted by x1i,x2i,…,xnii, then the center (clustericenter) is determined as(22)clustericenter=∑i=1nixini.The radius (R) of a cluster is defined as the mean distance (Euclidean) of the members from the center of the cluster. Thus, R can be written as (23)Ri=∑p=1nixp-clustericenterni.(II) Clusters in the Proposed Algorithm CMOABC. In the proposed algorithm, the stochastically generated population is partitioned into n subpopulations based on the widely adopted K-means cluster method (Algorithm 1) [44]. The number of clusters is determined by the predefined set G={g1,g2,…,gm}, where g1>g2>⋯>gm. Every cluster operates as the modified ABC introduced in the above section. During optimization, it may happen that two or more clusters come close to each other or get overlapped to a high degree. Then, they will practically search the same domain of the functional landscape. To avoid this scenario, the distances between each two clusters are calculated as follows:(24)Dis_cluster=clustericenter-Nei_clustericenter,where Dis_cluster is the distance between one cluster and its neighbor and Nei_clustericenter is the center of the ith cluster’s neighbor. clustericenter is the center of the ith cluster.Algorithm 1:Pseudocode of CMOABC.
//Step  1. Initialization
Generalize a population of NP individuals in the search region randomly;
Set the numbers of clustersG; Create the external archive EA; initial the current iteration iter = 1;
//Step  2. Loop
while stopping criteria are not satisfied do
Calculate TI according to (26)
Decide the numbergi of clusters according to G
Partition the whole population based onK-means clustering
For all cluster
ComputeDis_cluster according to  (24)
Decide whether some clusters need to be removed according toDISm
Selectxup individuals based on non-domination.
End for
For iter = 1 : TI
Update the individuals’ position in each cluster according to  (17)
Update the external archive and updateNum_Update
iter = iter + 1;
End for
Calculate the number of individuals in each cluster needing to be regenerated, according to (27)
Generate a certain number of the new individuals in every cluster;
End whileIf the distance between one cluster and its neighbors is smaller than the specific distanceDISm, one of the clusters will be removed and its nondomination solutions are stored:(25)DISm=0.2∗min⁡⁡Ri,Ri_neighbor,where Ri is the radius of clusteri and Ri_neighbor is the radius of the neighbors of clusteri.During optimization, information communication is not an inconsiderable aspect. For keeping the good ability of exploration, in the proposed algorithm, there is no information communication amonggi clusters in specific iterations. Therefore, in order to exchange information among individuals, the whole population is repartitioned into gi+1 clusters based on K-means clustering after each TI iteration, where gi and gi+1 are orderly chosen from the predefined set G={g1,g2,…,gm}. That is, the individuals in a cluster may be distributed into different new clusters when the number of the clusters is changing. Moreover, to balance the exploration and exploitation, the value of TI is not a constant. To strengthen exploration, the rate of sharing information should keep a low value in the first half of search process. To improve the ability of exploitation, in the second half, the interval iteration TI changes to a small value. The formula of TI is listed as follows:(26)TI=floor0.03∗itermax⁡if  iter≤0.5∗itermax⁡,floor0.06∗itermax⁡if  iter>0.5∗itermax⁡,where itermax⁡ is the maximum iterations; iter is the current iteration.After eachTI iteration, a certain number of individuals in each cluster should be regenerated according to this cluster’s contribution to the external archive that is used to store the nondomination solutions. For the jth cluster, the number of solutions updating to the external archive during each TI iteration is recorded in Num_Update(j). Then, according to its position in the sort of Num_Update, the number of individuals needed to regenerate in the jth cluster is calculated in (27). The individuals that will be removed in cluster j are determined by nondomination sort:(27)Num_regeneratej=Sort_Updatejgi∗Num_indj2,where Num_indj is the number of individuals in the jth cluster; Sort_Updatej indicates the jth cluster’s position in the sort of Num_Update; and gi is the current number of clusters.Additionally, to prevent the size of a cluster (Num_ind) from becoming too large, we set a limit Cmax⁡ as the maximum of cluster size. If the number of individuals in a cluster exceeds the limit, the best Cmax⁡ solutions are reserved based on nondomination and others are removed to other clusters stochastically. Setting formula of Cmax⁡ is shown in (28)Cmax⁡=1.5∗NPgi,where NP is the number of individuals in the population and gi is the number of clusters in current.


### 3.3. External Archive and Crowding Distance

#### 3.3.1. External Archive

In multiobjective optimization, in the absence of preference information, none of the solutions can be said to be better than the others. Therefore, in order to keep the best solutions, a fixed-sized external archive is used by the CMOABC algorithm. The external archive is used to keep a historical record of the nondominated vectors found along the search process [6, 45].


#### 3.3.2. Crowding Distance

In order to prevent the number of nondominated solutions from exceeding the size of EA, the crowding distance [4] is used to remove the crowded solutions. But in the proposed algorithm, the rules of removing the crowded solutions are modified. In each iteration, the distances between each nondominant solution and adjacent nondominant solution are calculated first, whether the size of external archive reaches the defined size or not. The detailed improvement is described as follows.(1)
The size of external archive exceeds the defined size of EA: calculate the crowding distance and remove the crowded solutions [4].(2)
The size of external archive does not outnumber the defined size: the crowding distance will be compared with a metric DIS. DIS is defined in (29); those nondominant solutions in the EA whose distances are all less than DIS will be removed from the EA [44]. Consider(29)DIS=∑i=1n-1fi-fi-1n+1.


#### 3.3.3. Constraints Handling Rules

To solve constraints, this paper uses a modified way to handle constraints which is introduced by Mohamed and Sabry [45]. Based on the initial constraints handling methods [46], three modified rules are as follows.(1)
The trial vector is feasible and the target vector is not.(2)
The trial vector and target vector are both feasible, but the trial vector has fitness value smaller than or equal to the corresponding target vector.(3)
The trial vector and target vector are both infeasible, but the trial vector has overall constraint violation smaller than or equal to the corresponding target vector.The above selection procedure allows the trial vector to be entered in the new population if it has the same amount of constraint violation or objective function value as the target vector. Therefore, this simple modification can help the algorithm to spread out and pass through the search space, so the algorithm can escape from stagnation.


### 3.4. Benchmark Test

#### 3.4.1. Test Function

To fully evaluate the performance of the CMOABC algorithm without a biased conclusion towards some chosen problems, in this paper, six commonly recognized benchmark functions have been used. Four of these test problems ZDT1, ZDT2, ZDT3, and ZDT6 [32] are of two objectives, while the other two, DLTZ2 and DLTZ6 [33], are of three objectives.


#### 3.4.2. Performance Measures

In order to facilitate the quantitative assessment of the performance of a multiobjective optimization algorithm, two performance metrics are taken into consideration: (1) convergence metricλ and (2) diversity metric η [4].Convergence Metric. This metric measures the extent of convergence to a known set of Pareto-optimal solutions:(30)λ=∑i=1NdiN,where N is the number of nondominated solutions obtained with an algorithm and di is the Euclidean distance between each of the nondominated solutions and the nearest member of the true Pareto-optimal front. To calculate this metric, we find a set of H=500 uniformly spaced solutions from the true Pareto-optimal front in the objective space. For each solution obtained with an algorithm, we compute the minimum Euclidean distance of it from H chosen solutions on the Pareto-optimal front. The average of these distances is used as the convergence metric λ. Figure 2 shows the calculation procedure of this metric.Figure 2
Convergence metricλ.Diversity Metric. This metric measures the extent of spread achieved among the obtained solutions. Here, we are interested in getting a set of solutions that spans the entire Pareto-optimal region. This metric is defined as(31)η=df+dl+∑i=1N-1di-d¯df+dl+N-1d¯,where di is the Euclidean distance between consecutive solutions in the obtained nondominated set of solutions and N is the number of nondominated solutions obtained by an algorithm. d¯ is the average value of these distances. df and dl are the Euclidean distances between the extreme solutions and the boundary solutions of the obtained nondominated set, as depicted in Figure 3.Figure 3
Diversity metricη.


#### 3.4.3. Experimental Setting

The effectiveness of CMOABC will be demonstrated on various standard test problems in this section. In addition, the proposed algorithm will be compared with the nondominated sorting genetic algorithm II (NSGA-II) [4], the multiobjective particle swarm optimization (MOPSO) [47], and the multiobjective artificial bee colony algorithm (MOABC) [13]. In order to compare the different algorithms with a fair time measure, the number of function evaluations (Fes = 10,000) is used for the termination criterion.For two test examples presented in this study, the CMOABC algorithm parameters are set as follows: the population size is 500 and the setG is defined as {50,20,10,5}.For NSGA-II, the population size, crossover, and mutation probabilities are selected as 500, 0.8, and 0.2, respectively, for the two test examples. For MOPSO, the population size, mutation rate, and divisions for the adaptive grid are selected as 500, 0.5, and 30. For MOABC, the population size is selected as 500, and for other parameter settings refer to [13].


#### 3.4.4. Two Objective Functions’ Results

The results reported in terms of the best, worst, average, median, and standard deviation of the performance measures are listed in Tables1–4. Figures 4–7 show the optimal front obtained by four algorithms for each two objective functions, respectively. The continuous lines represent the Pareto-optimal front, while mark spots represent found nondominated solutions by the algorithms.Table 1
Comparison of performance on ZDT1.
ZDT1
CMOABC
MOABC
NSGA-II
MOPSO
Convergence metric
Average
8
.
4
93
2
e
-
0
0
4
1
.
3
87
9
e
-
0
0
2
2
.
3
01
2
e
-
00
1
1
.
6
89
2
e
-
00
3
Median
8
.
2
39
8
e
-
0
0
4
1
.
9
08
6
e
-
0
0
2
1
.
4
32
1
e
-
00
1
1
.
4
30
9
e
-
00
3
Min
7
.
5
65
6
e
-
0
0
4
9
.
8
86
7
e
-
0
0
3
7
.
2
36
2
e
-
00
2
8
.
4
79
8
e
-
00
4
Max
1
.
1
02
1
e
-
0
0
3
3
.
6
60
9
e
-
00
2
8
.
8
03
1
e
-
00
1
2
.
1
97
9
e
-
00
3
Std.
8
.
0
98
8
e
-
0
0
5
1
.
0
15
7
e
-
00
2
2
.
3
71
5
e
-
00
1
6
.
7
12
1
e
-
00
4
Diversity metric
Average
6
.
7
12
9
e
-
0
0
2
2
.
8
96
2
e
-
0
0
1
7
.
5
98
9
e
-
0
0
1
6
.
3
10
2
e
-
0
0
1
Median
6
.
2
48
3
e
-
0
0
1
3
.
3
66
3
e
-
0
0
1
5
.
1
98
6
e
-
0
0
1
6
.
8
96
3
e
-
0
0
1
Min
5
.
3
26
6
e
-
0
0
1
2
.
9
63
1
e
-
0
0
1
4
.
6
01
2
e
-
0
0
1
6
.
5
6
1
1
e
-
0
0
1
Max
6
.
2
30
9
e
-
0
0
1
3
.
7
35
4
e
-
0
0
1
9
.
6
12
3
e
-
0
0
1
7
.
2
02
2
e
-
0
0
1
Std.
4
.
4
15
6
e
-
0
0
2
2
.
4
21
2
e
-
0
0
2
1
.
5
30
1
e
-
0
0
1
2
.
4
26
2
e
-
0
0
2Table 2
Comparison of performance on ZDT2.
ZDT2
CMOABC
MOABC
NSGA-II
MOPSO
Convergence metric
Average
2
.
0
30
6
e
-
0
0
4
1
.
0
15
1
e
-
0
0
3
8
.
1
21
8
e
-
0
0
2
1
.
5
20
1
e
-
0
0
1
Median
3
.
3
72
7
e
_
00
4
9
.
4
90
1
e
-
0
0
4
9
.
5
23
1
e
-
0
0
2
7
.
9
93
1
e
-
0
0
2
Min
9
.
2
46
6
e
_
00
5
6
.
9
88
3
e
-
0
0
4
2
.
1
83
4
e
-
0
0
2
6
.
8
64
1
e
-
0
0
4
Max
1
.
8
27
6
e
-
0
0
3
1
.
0
92
4
e
-
0
0
3
1
.
1
63
2
e
-
0
0
1
6
.
1
07
3
e
-
0
0
1
Std.
6
.
0
91
9
e
-
0
0
4
2
.
5
06
1
e
-
0
0
4
3
.
3
87
6
e
-
0
0
2
6
.
5
66
3
e
-
0
0
1
Diversity metric
Average
6
.
4
83
2
e
-
0
0
2
2
.
8
95
7
e
-
0
0
1
6
.
0
09
7
e
-
0
0
1
6
.
2
90
7
e
-
0
0
1
Median
6
.
7
37
9
e
-
00
2
3
.
0
01
5
e
-
0
0
1
5
.
9
50
5
e
-
0
0
1
6
.
3
0
1
1
e
-
0
0
1
Min
6
.
1
05
8
e
-
00
2
2
.
7
63
1
e
-
0
0
1
5
.
4
03
8
e
-
0
0
1
5
.
1
87
6
e
-
0
0
1
Max
5
.
8
34
9
e
-
0
0
1
3
.
3
12
0
e
-
0
0
1
6
.
8
72
4
e
-
0
0
1
7
.
2
27
2
e
-
0
0
1
Std.
1
.
9
45
0
e
-
0
0
1
2
.
2
09
2
e
-
0
0
2
4
.
5
62
8
e
-
0
0
2
7
.
6
24
1
e
-
0
0
2Table 3
Comparison of performance on ZDT3.
ZDT3
CMOABC
MOABC
NSGA-II
MOPSO
Convergence metric
Average
4
.
3
25
6
e
-
0
0
4
1
.
0
01
2
e
-
0
0
3
2
.
9
47
6
e
+
0
0
0
1
.
6
42
1
e
-
0
0
2
Median
1
.
5
90
8
e
-
0
0
4
9
.
7
59
7
e
-
0
0
4
2
.
8
10
5
e
+
0
0
0
2
.
9
83
1
e
-
0
0
3
Min
9
.
3
66
4
e
-
0
0
5
8
.
9
18
7
e
-
0
0
4
5
.
0
66
8
e
-
0
0
1
1
.
8
27
7
e
-
0
0
3
Max
3
.
8
93
0
e
-
0
0
3
2
.
5
60
9
e
-
0
0
3
4
.
8
21
8
e
+
0
0
0
5
.
2
1
1
4
e
-
0
0
2
Std.
1
.
2
97
7
e
-
0
0
3
2
.
4
48
1
e
-
0
0
4
1
.
6
28
8
e
+
0
0
0
2
.
6
67
1
e
-
0
0
2
Diversity metric
Average
7
.
5
59
2
e
-
0
0
2
6
.
4
35
2
e
-
0
0
1
9
.
9
80
2
e
-
0
0
1
8
.
0
00
2
e
-
0
0
1
Median
6
.
9
82
5
e
-
0
0
2
6
.
4
18
5
e
-
0
0
1
8
.
9
97
9
e
-
0
0
1
6
.
9
87
1
e
-
0
0
1
Min
2
.
11
1
8
e
-
0
0
2
5
.
9
63
1
e
-
0
0
1
4
.
6
61
7
e
-
0
0
1
5
.
0
80
5
e
-
0
0
1
Max
6
.
8
03
2
e
-
0
0
1
6
.
8
11
1
e
-
0
0
1
1
.
0
54
2
e
+
0
0
0
8
.
9
39
2
e
-
0
0
1
Std.
2
.
2
67
7
e
-
0
0
1
1
.
2
14
7
e
-
0
0
2
2
.
3
56
7
e
-
0
0
2
3
.
3
28
8
e
-
0
0
2Table 4
Comparison of performance on ZDT6.
ZDT6
CMOABC
MOABC
NSGA-II
MOPSO
Convergence metric
Average
3
.
0
22
1
e
-
0
0
4
2
.
8
62
3
e
-
0
03
3
.
1
40
3
e
-
0
0
1
4
.
9
31
5
e
-
0
0
1
Median
3
.
6
05
9
e
-
0
0
4
9
.
2
56
7
e
-
0
0
4
2
.
4
21
4
e
-
0
0
1
8
.
6
06
1
e
-
0
0
4
Min
9
.
2
46
6
e
-
0
0
5
1
.
2
07
3
e
-
0
0
5
6
.
6
36
0
e
-
0
0
2
6
.
6
70
9
e
-
0
0
4
Max
2
.
7
19
9
e
-
0
0
3
1
.
9
54
2
e
-
0
0
2
8
.
7
22
0
e
-
0
0
1
6
.
2
1
1
4
e
-
0
0
2
Std.
9
.
0
66
2
e
-
0
0
4
7
.
8
52
6
e
-
0
0
3
2
.
8
45
2
e
-
0
0
1
6
.
6
57
1
e
-
0
0
3
Diversity metric
Average
6
.
8
40
7
e
-
0
0
2
5
.
9
75
2
e
-
0
0
1
7
.
5
90
5
e
-
0
0
1
6
.
5
90
2
e
-
0
0
1
Median
6
.
7
38
1
e
-
0
0
2
4
.
0
18
5
e
-
0
0
1
7
.
4
81
2
e
-
0
0
1
6
.
6
97
1
e
-
0
0
1
Min
6
.
2
05
9
e
-
0
0
2
3
.
9
85
1
e
-
0
0
1
5
.
7
13
5
e
-
0
0
1
5
.
8
89
5
e
-
0
0
1
Max
6
.
1
56
7
e
-
0
0
1
1
.
1
45
5
e
+
0
0
0
9
.
4
40
3
e
-
0
0
1
7
.
0
70
2
e
-
0
0
1
Std.
2
.
0
52
2
e
-
0
0
1
2
.
2
14
7
e
-
0
0
1
9
.
7
50
0
e
-
0
0
2
4
.
4
95
8
e
-
0
0
2Figure 4
Pareto fronts obtained by CMOABC, MOPSO, MOABC, and NSGA-II on ZDT1.Figure 5
Pareto fronts obtained by CMOABC, MOPSO, MOABC, and NSGA-II on ZDT2.Figure 6
Pareto fronts obtained by CMOABC, MOPSO, MOABC, and NSGA-II on ZDT3.Figure 7
Pareto fronts obtained by CMOABC, MOPSO, MOABC, and NSGA-II on ZDT6.On ZDT1 function, when given 10,000 function evaluations for four algorithms, Figure4 shows that CMOABC, MOPSO, NSGA-II, and MOABC have found a well distributed and diverse solution set in the entire Pareto-optimal region. However, NSGA-II and MOABC cannot achieve the true Pareto front for ZDT1. Table 1 shows that the performance of CMOABC in convergence metric is of three orders of magnitude better than NSGA-II, two orders better than MOABC, and one order better than MOPSO. For diversity metric, it can be observed that all algorithms have almost the same performance on ZDT1.On ZDT2 function, Figure5 shows that MOPSO and NSGA-II produce poor results on this test function and both of them cannot achieve the true Pareto front, and the results of the performance measures show that CMOABC and MOABC have better performance compared to the MOPSO and NSGA-II. From Table 3, it can be noticed that the performance of CMOABC and MOABC in convergence metric and diversity metric is better than that of MOPSO and NSGA-II. Moreover, CMOABC outperforms MOABC by one order of magnitude in terms of diversity metric.On ZDT3 function, from Figure6, it can be noticed that NSGA-II has difficulties in converging the true Pareto front, and other algorithms find a better converged set of nondominated solutions in ZDT3 compared to NSGA-II. Further, Table 3 also shows that NSGA-II is worse than CMOABC, MOPSO, and MOABC in convergence metric and diversity metric.On ZDT6 function, Figure7 and Table 4 show that the algorithms achieve similar performance as on ZDT3.


#### 3.4.5. Three Objective Functions

The optimization results of CMOABC, MOPSO, MOABC, and NSGA-II algorithms for three objective DTLZ series problems are shown in Tables5 and 6. The Pareto-optimal region and all nondominated solutions obtained by the four algorithms are shown in Figures 8–11.Table 5
Comparison of performance on DTLZ2.
DTLZ2
CMOABC
MOABC
NSGA-II
MOPSO
Convergence metric
Average
3
.
0
89
5
e
-
0
0
4
8
.
0
02
3
e
-
00
3
8
.
7
89
6
e
-
0
0
2
3
.
8
68
4
e
-
0
0
3
Median
3
.
6
05
9
e
-
0
0
4
6
.
5
53
7
e
-
0
03
8
.
0
43
8
e
-
0
0
2
3
.
8
43
1
e
-
0
0
3
Min
9
.
8
02
1
e
-
0
0
5
5
.
2
18
3
e
-
0
03
6
.
1
40
9
e
-
0
0
2
3
.
7
63
1
e
-
0
0
3
Max
1
.
9
76
5
e
-
0
0
3
1
.
5
69
0
e
-
0
02
1
.
0
62
5
e
-
0
0
1
3
.
0
77
2
e
-
0
0
3
Std.
9
.
7
69
9
e
-
0
0
4
2
.
5
38
0
e
-
0
03
1
.
6
79
8
e
-
0
0
2
8
.
6
92
4
e
-
0
0
5
Diversity metric
Average
4
.
4
90
6
e
-
0
0
2
3
.
9
89
6
e
-
0
0
1
6
.
1
30
2
e
-
0
0
1
4
.
2
63
7
e
-
0
0
1
Median
6
.
7
38
1
e
-
0
0
2
3
.
9
68
5
e
-
0
0
1
5
.
2
99
9
e
-
0
0
1
4
.
1
52
1
e
-
0
0
1
Min
6
.
2
05
9
e
-
0
0
2
2
.
9
57
2
e
-
0
0
1
3
.
6
61
7
e
-
0
0
1
3
.
8
79
0
e
-
0
0
1
Max
1
.
4
90
6
e
-
0
0
1
5
.
2
56
1
e
-
0
0
1
7
.
2
54
2
e
-
0
0
1
4
.
5
61
8
e
-
0
0
1
Std.
1
.
8
64
9
e
-
0
0
2
3
.
2
23
4
e
-
0
0
2
1
.
7
56
7
e
-
0
0
2
2
.
1
69
2
e
-
0
0
2Table 6
Comparison of performance on DTLZ6.
DTLZ6
CMOABC
MOABC
NSGA-II
MOPSO
Convergence metric
Average
1
.
8
36
5
e
-
0
0
3
1
.
6
50
6
e
-
0
0
2
6
.
0
47
0
e
-
0
0
1
1
.
7
01
2
e
-
0
0
2
Median
3
.
6
21
9
e
-
0
0
4
3
.
4
04
1
e
-
0
0
3
3
.
2
73
3
e
-
0
0
1
1
.
8
00
2
e
-
0
0
2
Min
1
.
2
43
6
e
-
0
0
4
5
.
2
73
1
e
-
0
0
3
1
.
7
49
8
e
-
0
0
1
8
.
2
98
5
e
-
0
0
2
Max
9
.
1
82
6
e
-
0
0
3
2
.
8
04
3
e
-
0
0
2
2
.
9
72
8
e
+
0
0
0
2
.
8
96
1
e
-
0
0
2
Std.
4
.
1
06
6
e
-
0
0
3
4
.
9
12
7
e
-
0
0
3
8
.
4
99
7
e
-
0
0
1
8
.
3
34
6
e
-
0
0
3
Diversity metric
Average
9
.
0
41
3
e
-
0
0
2
2
.
9
02
1
e
-
0
0
1
5
.
6
47
7
e
-
0
0
1
5
.
5
00
2
e
-
0
0
1
Median
9
.
0
32
8
e
-
0
0
2
2
.
8
92
6
e
-
0
0
1
5
.
5
86
6
e
-
0
0
1
5
.
4
97
6
e
-
0
0
1
Min
7
.
8
77
2
e
-
0
0
2
9
.
8
26
6
e
-
0
0
2
5
.
1
84
9
e
-
0
0
1
5
.
2
01
5
e
-
0
0
1
Max
1
.
1
52
5
e
-
0
0
1
6
.
3
23
1
e
-
0
0
1
6
.
3
24
5
e
-
0
0
1
5
.
7
30
2
e
-
0
0
1
Std.
2
.
0
64
5
e
-
0
0
2
2
.
8
02
5
e
-
0
0
2
3
.
7
46
9
e
-
0
0
2
1
.
7
49
8
e
-
0
0
2Figure 8
The true Pareto-optimal front on DTLZ2.Figure 9
Pareto fronts obtained by CMOABC, MOPSO, MOABC, and NSGA-II on DTLZ2.Figure 10
The true Pareto-optimal front on DTLZ6.Figure 11
Pareto fronts obtained by CMOABC, MOPSO, MOABC, and NSGA-II on DTLZ6.From Figure9, it is observed that all algorithms obtain similar performance on DTLZ2. However, Table 5 shows the performance of NSGA-II is worse than that of the other three algorithms in convergence, while the performance of CMOABC is better than that of the other three algorithms in diversity.From Figure11, it can be seen that CMOABC and MOABC have better performance on DTLZ6. The CMOABC especially has attained the nondominated set properly. The front obtained from the other two algorithms cannot be uniformly distributed on the true Pareto-optimal front on DTLZ6.Table6 also shows that MOABC has obtained better convergence and diversity than NSGA-II and MOPSO on DTLZ6, but it is found to be second to the CMOABC algorithm.On the whole, for two and three objective problems, CMOABC has a more competitive performance than MOPSO, MOABC, and NSGA-II on both diversity metric and convergence metric. This is due to the fact that, in the proposed CMOABC algorithm, the modified employed bees’ phase with a self-adaptive strength of information communication effectively enhances its abilities of exploration in the early stage and exploitation in the late stage. Moreover, benefiting from the swarm decomposition operation with a dynamicK-means clustering method, CMOABC could ensure more sufficient exploitation around the potential solutions in the less crowded regions. When the number of clusters is large, the probability of the individuals with similar performances allocated into the same subcomponent increases, which improves the ability of global search. When the number of clusters becomes small, the individuals trapped into local optimal may move toward the current best individual of their cluster, which effectively avoids the inefficiency of local search. In addition, the regeneration operation also enriches the diversity of external archive. Therefore, CMOABC can reasonably accommodate a considerable potential for solving more complex multiobjective optimization problems.


## 3.1. Canonical ABC Algorithm

The artificial bee colony (ABC) algorithm, proposed by Karaboga and Basturk [39] and further developed by Basturk and Akay [11, 40] for real-parameter optimization, which simulates the intelligent foraging behavior of a honeybee swarm, is one of the most recently introduced swarm-based optimization techniques.The entire bee colony contains three groups of bees: employed bees, onlookers, and scouts. Employed bees explore the specific food sources; meanwhile they pass the food information to onlooker bees. The number of employed bees is equal to that of food sources; in other words, each food source owns only one employed bee. Then onlooker bees choose good food sources based on the received information and then further exploit the food near their selected food source. The food source with higher quality would have a larger opportunity to be selected by onlookers. There is a control parameter called “limit” in the canonical ABC algorithm. If a food source is not improved anymore whenlimit is exceeded, it is assumed to be abandoned by its employed bee and the employed bee associated with that food source becomes a scout to search for a new food source randomly. The fundamental mathematic representations are listed as follows [41].Step 1 (initialization phase).
In initialization phase, a group of food sources are generated randomly in the search space using the following equation:(14)xi,j=xjmin⁡+rand⁡0,1xjmax⁡-xjmin⁡,where i=1,2,…,SN, j=1,2,…,D. SN is the number of food sources. D is the number of variables, that is, problem dimension. xjmin⁡ and xjmax⁡ are the lower and upper bounds of the jth variable, respectively.Step 2 (employed bees’ phase).
In the employed bees’ phase, the neighbor food source (candidate solution) can be generated from the old food source of each employed bee in its memory using the following expression:(15)vi,j=xi,j+δxi,j-xk,j,where k is a randomly selected food source and must be different from i; j is randomly chosen indexes; δ is a random number in range [-1,1].Step 3 (onlooker bees’ phase).
In the onlooker bees’ phase, an onlooker bee selects a food source depending on the probability value associated with that food source andPi can be calculated as follows:(16)Pi=fitnessi∑j=1SNfitnessj,where fitnessj is the fitness value of the jth solution.Step 4 (scout bees’ phase).
In scout bees’ phase, if a food source cannot be improved further through a predetermined cycle (called “limit” in ABC), the food source is supposed to be abandoned. The employed bee subsequently becomes a scout. A new food source will be produced randomly in the search space using (14).The employed, onlooker, and scout bees’ phase will recycle until the termination condition is met.


## 3.2. CMOABC for Optimization

This section presents the detailed description of the proposed multiobjective artificial bee colony algorithm based onK-means clustering (named CMOABC). The proposed algorithm makes improvements on two aspects: (1) modifying the method of information communication in employed bees’ phase and (2) partitioning the population using K-means clustering. The number of the clusters will be chosen from the predefined set G. After each specific iteration, when the defined criteria are met, parts of individuals in clusters will be removed and an equal number of new individuals are regenerated randomly. Then the number of clusters will change into a new value of the set G and the population will be reclustered based on K-means clustering.


### 3.2.1. Modified Employed Bees’ Phase of ABC

To increase the exploitation of ABC algorithm and fasten the convergence rate, the employed bees’ phase is modified. The new food source (candidate solution) is generated using the following way:(17)vij=xij+ϕxij-xkjrand≥Cr,vij=xij+φxbest,j-xijrand<Cr,(18)Cr=iteritermax⁡,where xbest,j is the position of the global best solution and the term xbest,j can drive the new candidate solution to the global best solution, iter is the current iteration, and itermax⁡ is the maximum iteration.The parameterCr in (18) plays an important role in balancing the exploration and exploitation of the candidate solution search. If Cr=0, (18) is identical to (15). With the increase of Cr, the probability of the candidate solution learning to the best solution increases correspondingly. In this way, in the beginning of optimization process, the proposed algorithm operates as the canonical artificial bee colony algorithm, which can well keep the population diversity; in the end of optimization process, the modified ABC algorithm has a considerable improvement on both convergence rate and local search.Furthermore, in (17), the different learning intensities, which are determined by ϕ and φ, are adopted. We adopt the beta probability distribution [42] to tune these two parameters, which is based on the similar strategy of tuning parameters in [43]. The beta distribution is flexible for modeling data that are measured in a continuous scale on a truncated interval in range [0,1], since its density is a versatile way to represent different shapes depending on the values of the two parameters that index the distribution. The probability density (Bd_pd) of the beta distribution in the range [0,1] is calculated as follows:(19)Bd_pd=xα-11-xβ-1Bα,β,where B is the beta function.The mean (E) and variance (Var) of the beta distribution are calculated as follows:(20)EX=αα+β,Var⁡X=αβα+β2α+β+1.The values ofϕ and φ are tuned by MATLAB∗ script (see (21)). An example of ϕ and φ values generated using beta distribution is shown in Figure 1. Consider(21)α=rand,β=0.7∗iteritermax⁡+0.2∗rand,ϕ=betarnda,b,11,α=rand,β=0.6∗iteritermax⁡+0.1∗rand,φ=betarnda,b,11,where itermax⁡ denotes the max iteration and iter is the current iteration.Figure 1
(a)F and (b) CR values for a run generated using beta distribution in the beta distribution.
(a)
(b)


### 3.2.2. Cluster Setting

TheK-means clustering method is employed to partition the population into subpopulations. The basic concepts of K-means clustering are presented firstly, and then we will give a detailed description of its application in our proposed algorithm.(I) Basic Parameters of
K
-Means. The cluster centers are substituted for center positions of food sources and the formula of computing the centers is shown in (22). If the ith cluster contains ni members and the members are denoted by x1i,x2i,…,xnii, then the center (clustericenter) is determined as(22)clustericenter=∑i=1nixini.The radius (R) of a cluster is defined as the mean distance (Euclidean) of the members from the center of the cluster. Thus, R can be written as (23)Ri=∑p=1nixp-clustericenterni.(II) Clusters in the Proposed Algorithm CMOABC. In the proposed algorithm, the stochastically generated population is partitioned into n subpopulations based on the widely adopted K-means cluster method (Algorithm 1) [44]. The number of clusters is determined by the predefined set G={g1,g2,…,gm}, where g1>g2>⋯>gm. Every cluster operates as the modified ABC introduced in the above section. During optimization, it may happen that two or more clusters come close to each other or get overlapped to a high degree. Then, they will practically search the same domain of the functional landscape. To avoid this scenario, the distances between each two clusters are calculated as follows:(24)Dis_cluster=clustericenter-Nei_clustericenter,where Dis_cluster is the distance between one cluster and its neighbor and Nei_clustericenter is the center of the ith cluster’s neighbor. clustericenter is the center of the ith cluster.Algorithm 1:Pseudocode of CMOABC.
//Step  1. Initialization
Generalize a population of NP individuals in the search region randomly;
Set the numbers of clustersG; Create the external archive EA; initial the current iteration iter = 1;
//Step  2. Loop
while stopping criteria are not satisfied do
Calculate TI according to (26)
Decide the numbergi of clusters according to G
Partition the whole population based onK-means clustering
For all cluster
ComputeDis_cluster according to  (24)
Decide whether some clusters need to be removed according toDISm
Selectxup individuals based on non-domination.
End for
For iter = 1 : TI
Update the individuals’ position in each cluster according to  (17)
Update the external archive and updateNum_Update
iter = iter + 1;
End for
Calculate the number of individuals in each cluster needing to be regenerated, according to (27)
Generate a certain number of the new individuals in every cluster;
End whileIf the distance between one cluster and its neighbors is smaller than the specific distanceDISm, one of the clusters will be removed and its nondomination solutions are stored:(25)DISm=0.2∗min⁡⁡Ri,Ri_neighbor,where Ri is the radius of clusteri and Ri_neighbor is the radius of the neighbors of clusteri.During optimization, information communication is not an inconsiderable aspect. For keeping the good ability of exploration, in the proposed algorithm, there is no information communication amonggi clusters in specific iterations. Therefore, in order to exchange information among individuals, the whole population is repartitioned into gi+1 clusters based on K-means clustering after each TI iteration, where gi and gi+1 are orderly chosen from the predefined set G={g1,g2,…,gm}. That is, the individuals in a cluster may be distributed into different new clusters when the number of the clusters is changing. Moreover, to balance the exploration and exploitation, the value of TI is not a constant. To strengthen exploration, the rate of sharing information should keep a low value in the first half of search process. To improve the ability of exploitation, in the second half, the interval iteration TI changes to a small value. The formula of TI is listed as follows:(26)TI=floor0.03∗itermax⁡if  iter≤0.5∗itermax⁡,floor0.06∗itermax⁡if  iter>0.5∗itermax⁡,where itermax⁡ is the maximum iterations; iter is the current iteration.After eachTI iteration, a certain number of individuals in each cluster should be regenerated according to this cluster’s contribution to the external archive that is used to store the nondomination solutions. For the jth cluster, the number of solutions updating to the external archive during each TI iteration is recorded in Num_Update(j). Then, according to its position in the sort of Num_Update, the number of individuals needed to regenerate in the jth cluster is calculated in (27). The individuals that will be removed in cluster j are determined by nondomination sort:(27)Num_regeneratej=Sort_Updatejgi∗Num_indj2,where Num_indj is the number of individuals in the jth cluster; Sort_Updatej indicates the jth cluster’s position in the sort of Num_Update; and gi is the current number of clusters.Additionally, to prevent the size of a cluster (Num_ind) from becoming too large, we set a limit Cmax⁡ as the maximum of cluster size. If the number of individuals in a cluster exceeds the limit, the best Cmax⁡ solutions are reserved based on nondomination and others are removed to other clusters stochastically. Setting formula of Cmax⁡ is shown in (28)Cmax⁡=1.5∗NPgi,where NP is the number of individuals in the population and gi is the number of clusters in current.


## 3.2.1. Modified Employed Bees’ Phase of ABC

To increase the exploitation of ABC algorithm and fasten the convergence rate, the employed bees’ phase is modified. The new food source (candidate solution) is generated using the following way:(17)vij=xij+ϕxij-xkjrand≥Cr,vij=xij+φxbest,j-xijrand<Cr,(18)Cr=iteritermax⁡,where xbest,j is the position of the global best solution and the term xbest,j can drive the new candidate solution to the global best solution, iter is the current iteration, and itermax⁡ is the maximum iteration.The parameterCr in (18) plays an important role in balancing the exploration and exploitation of the candidate solution search. If Cr=0, (18) is identical to (15). With the increase of Cr, the probability of the candidate solution learning to the best solution increases correspondingly. In this way, in the beginning of optimization process, the proposed algorithm operates as the canonical artificial bee colony algorithm, which can well keep the population diversity; in the end of optimization process, the modified ABC algorithm has a considerable improvement on both convergence rate and local search.Furthermore, in (17), the different learning intensities, which are determined by ϕ and φ, are adopted. We adopt the beta probability distribution [42] to tune these two parameters, which is based on the similar strategy of tuning parameters in [43]. The beta distribution is flexible for modeling data that are measured in a continuous scale on a truncated interval in range [0,1], since its density is a versatile way to represent different shapes depending on the values of the two parameters that index the distribution. The probability density (Bd_pd) of the beta distribution in the range [0,1] is calculated as follows:(19)Bd_pd=xα-11-xβ-1Bα,β,where B is the beta function.The mean (E) and variance (Var) of the beta distribution are calculated as follows:(20)EX=αα+β,Var⁡X=αβα+β2α+β+1.The values ofϕ and φ are tuned by MATLAB∗ script (see (21)). An example of ϕ and φ values generated using beta distribution is shown in Figure 1. Consider(21)α=rand,β=0.7∗iteritermax⁡+0.2∗rand,ϕ=betarnda,b,11,α=rand,β=0.6∗iteritermax⁡+0.1∗rand,φ=betarnda,b,11,where itermax⁡ denotes the max iteration and iter is the current iteration.Figure 1
(a)F and (b) CR values for a run generated using beta distribution in the beta distribution.
(a)
(b)


## 3.2.2. Cluster Setting

TheK-means clustering method is employed to partition the population into subpopulations. The basic concepts of K-means clustering are presented firstly, and then we will give a detailed description of its application in our proposed algorithm.(I) Basic Parameters of
K
-Means. The cluster centers are substituted for center positions of food sources and the formula of computing the centers is shown in (22). If the ith cluster contains ni members and the members are denoted by x1i,x2i,…,xnii, then the center (clustericenter) is determined as(22)clustericenter=∑i=1nixini.The radius (R) of a cluster is defined as the mean distance (Euclidean) of the members from the center of the cluster. Thus, R can be written as (23)Ri=∑p=1nixp-clustericenterni.(II) Clusters in the Proposed Algorithm CMOABC. In the proposed algorithm, the stochastically generated population is partitioned into n subpopulations based on the widely adopted K-means cluster method (Algorithm 1) [44]. The number of clusters is determined by the predefined set G={g1,g2,…,gm}, where g1>g2>⋯>gm. Every cluster operates as the modified ABC introduced in the above section. During optimization, it may happen that two or more clusters come close to each other or get overlapped to a high degree. Then, they will practically search the same domain of the functional landscape. To avoid this scenario, the distances between each two clusters are calculated as follows:(24)Dis_cluster=clustericenter-Nei_clustericenter,where Dis_cluster is the distance between one cluster and its neighbor and Nei_clustericenter is the center of the ith cluster’s neighbor. clustericenter is the center of the ith cluster.Algorithm 1:Pseudocode of CMOABC.
//Step  1. Initialization
Generalize a population of NP individuals in the search region randomly;
Set the numbers of clustersG; Create the external archive EA; initial the current iteration iter = 1;
//Step  2. Loop
while stopping criteria are not satisfied do
Calculate TI according to (26)
Decide the numbergi of clusters according to G
Partition the whole population based onK-means clustering
For all cluster
ComputeDis_cluster according to  (24)
Decide whether some clusters need to be removed according toDISm
Selectxup individuals based on non-domination.
End for
For iter = 1 : TI
Update the individuals’ position in each cluster according to  (17)
Update the external archive and updateNum_Update
iter = iter + 1;
End for
Calculate the number of individuals in each cluster needing to be regenerated, according to (27)
Generate a certain number of the new individuals in every cluster;
End whileIf the distance between one cluster and its neighbors is smaller than the specific distanceDISm, one of the clusters will be removed and its nondomination solutions are stored:(25)DISm=0.2∗min⁡⁡Ri,Ri_neighbor,where Ri is the radius of clusteri and Ri_neighbor is the radius of the neighbors of clusteri.During optimization, information communication is not an inconsiderable aspect. For keeping the good ability of exploration, in the proposed algorithm, there is no information communication amonggi clusters in specific iterations. Therefore, in order to exchange information among individuals, the whole population is repartitioned into gi+1 clusters based on K-means clustering after each TI iteration, where gi and gi+1 are orderly chosen from the predefined set G={g1,g2,…,gm}. That is, the individuals in a cluster may be distributed into different new clusters when the number of the clusters is changing. Moreover, to balance the exploration and exploitation, the value of TI is not a constant. To strengthen exploration, the rate of sharing information should keep a low value in the first half of search process. To improve the ability of exploitation, in the second half, the interval iteration TI changes to a small value. The formula of TI is listed as follows:(26)TI=floor0.03∗itermax⁡if  iter≤0.5∗itermax⁡,floor0.06∗itermax⁡if  iter>0.5∗itermax⁡,where itermax⁡ is the maximum iterations; iter is the current iteration.After eachTI iteration, a certain number of individuals in each cluster should be regenerated according to this cluster’s contribution to the external archive that is used to store the nondomination solutions. For the jth cluster, the number of solutions updating to the external archive during each TI iteration is recorded in Num_Update(j). Then, according to its position in the sort of Num_Update, the number of individuals needed to regenerate in the jth cluster is calculated in (27). The individuals that will be removed in cluster j are determined by nondomination sort:(27)Num_regeneratej=Sort_Updatejgi∗Num_indj2,where Num_indj is the number of individuals in the jth cluster; Sort_Updatej indicates the jth cluster’s position in the sort of Num_Update; and gi is the current number of clusters.Additionally, to prevent the size of a cluster (Num_ind) from becoming too large, we set a limit Cmax⁡ as the maximum of cluster size. If the number of individuals in a cluster exceeds the limit, the best Cmax⁡ solutions are reserved based on nondomination and others are removed to other clusters stochastically. Setting formula of Cmax⁡ is shown in (28)Cmax⁡=1.5∗NPgi,where NP is the number of individuals in the population and gi is the number of clusters in current.


## 3.3. External Archive and Crowding Distance

### 3.3.1. External Archive

In multiobjective optimization, in the absence of preference information, none of the solutions can be said to be better than the others. Therefore, in order to keep the best solutions, a fixed-sized external archive is used by the CMOABC algorithm. The external archive is used to keep a historical record of the nondominated vectors found along the search process [6, 45].


### 3.3.2. Crowding Distance

In order to prevent the number of nondominated solutions from exceeding the size of EA, the crowding distance [4] is used to remove the crowded solutions. But in the proposed algorithm, the rules of removing the crowded solutions are modified. In each iteration, the distances between each nondominant solution and adjacent nondominant solution are calculated first, whether the size of external archive reaches the defined size or not. The detailed improvement is described as follows.(1)
The size of external archive exceeds the defined size of EA: calculate the crowding distance and remove the crowded solutions [4].(2)
The size of external archive does not outnumber the defined size: the crowding distance will be compared with a metric DIS. DIS is defined in (29); those nondominant solutions in the EA whose distances are all less than DIS will be removed from the EA [44]. Consider(29)DIS=∑i=1n-1fi-fi-1n+1.


### 3.3.3. Constraints Handling Rules

To solve constraints, this paper uses a modified way to handle constraints which is introduced by Mohamed and Sabry [45]. Based on the initial constraints handling methods [46], three modified rules are as follows.(1)
The trial vector is feasible and the target vector is not.(2)
The trial vector and target vector are both feasible, but the trial vector has fitness value smaller than or equal to the corresponding target vector.(3)
The trial vector and target vector are both infeasible, but the trial vector has overall constraint violation smaller than or equal to the corresponding target vector.The above selection procedure allows the trial vector to be entered in the new population if it has the same amount of constraint violation or objective function value as the target vector. Therefore, this simple modification can help the algorithm to spread out and pass through the search space, so the algorithm can escape from stagnation.


## 3.3.1. External Archive

In multiobjective optimization, in the absence of preference information, none of the solutions can be said to be better than the others. Therefore, in order to keep the best solutions, a fixed-sized external archive is used by the CMOABC algorithm. The external archive is used to keep a historical record of the nondominated vectors found along the search process [6, 45].


## 3.3.2. Crowding Distance

In order to prevent the number of nondominated solutions from exceeding the size of EA, the crowding distance [4] is used to remove the crowded solutions. But in the proposed algorithm, the rules of removing the crowded solutions are modified. In each iteration, the distances between each nondominant solution and adjacent nondominant solution are calculated first, whether the size of external archive reaches the defined size or not. The detailed improvement is described as follows.(1)
The size of external archive exceeds the defined size of EA: calculate the crowding distance and remove the crowded solutions [4].(2)
The size of external archive does not outnumber the defined size: the crowding distance will be compared with a metric DIS. DIS is defined in (29); those nondominant solutions in the EA whose distances are all less than DIS will be removed from the EA [44]. Consider(29)DIS=∑i=1n-1fi-fi-1n+1.


## 3.3.3. Constraints Handling Rules

To solve constraints, this paper uses a modified way to handle constraints which is introduced by Mohamed and Sabry [45]. Based on the initial constraints handling methods [46], three modified rules are as follows.(1)
The trial vector is feasible and the target vector is not.(2)
The trial vector and target vector are both feasible, but the trial vector has fitness value smaller than or equal to the corresponding target vector.(3)
The trial vector and target vector are both infeasible, but the trial vector has overall constraint violation smaller than or equal to the corresponding target vector.The above selection procedure allows the trial vector to be entered in the new population if it has the same amount of constraint violation or objective function value as the target vector. Therefore, this simple modification can help the algorithm to spread out and pass through the search space, so the algorithm can escape from stagnation.


## 3.4. Benchmark Test

### 3.4.1. Test Function

To fully evaluate the performance of the CMOABC algorithm without a biased conclusion towards some chosen problems, in this paper, six commonly recognized benchmark functions have been used. Four of these test problems ZDT1, ZDT2, ZDT3, and ZDT6 [32] are of two objectives, while the other two, DLTZ2 and DLTZ6 [33], are of three objectives.


### 3.4.2. Performance Measures

In order to facilitate the quantitative assessment of the performance of a multiobjective optimization algorithm, two performance metrics are taken into consideration: (1) convergence metricλ and (2) diversity metric η [4].Convergence Metric. This metric measures the extent of convergence to a known set of Pareto-optimal solutions:(30)λ=∑i=1NdiN,where N is the number of nondominated solutions obtained with an algorithm and di is the Euclidean distance between each of the nondominated solutions and the nearest member of the true Pareto-optimal front. To calculate this metric, we find a set of H=500 uniformly spaced solutions from the true Pareto-optimal front in the objective space. For each solution obtained with an algorithm, we compute the minimum Euclidean distance of it from H chosen solutions on the Pareto-optimal front. The average of these distances is used as the convergence metric λ. Figure 2 shows the calculation procedure of this metric.Figure 2
Convergence metricλ.Diversity Metric. This metric measures the extent of spread achieved among the obtained solutions. Here, we are interested in getting a set of solutions that spans the entire Pareto-optimal region. This metric is defined as(31)η=df+dl+∑i=1N-1di-d¯df+dl+N-1d¯,where di is the Euclidean distance between consecutive solutions in the obtained nondominated set of solutions and N is the number of nondominated solutions obtained by an algorithm. d¯ is the average value of these distances. df and dl are the Euclidean distances between the extreme solutions and the boundary solutions of the obtained nondominated set, as depicted in Figure 3.Figure 3
Diversity metricη.


### 3.4.3. Experimental Setting

The effectiveness of CMOABC will be demonstrated on various standard test problems in this section. In addition, the proposed algorithm will be compared with the nondominated sorting genetic algorithm II (NSGA-II) [4], the multiobjective particle swarm optimization (MOPSO) [47], and the multiobjective artificial bee colony algorithm (MOABC) [13]. In order to compare the different algorithms with a fair time measure, the number of function evaluations (Fes = 10,000) is used for the termination criterion.For two test examples presented in this study, the CMOABC algorithm parameters are set as follows: the population size is 500 and the setG is defined as {50,20,10,5}.For NSGA-II, the population size, crossover, and mutation probabilities are selected as 500, 0.8, and 0.2, respectively, for the two test examples. For MOPSO, the population size, mutation rate, and divisions for the adaptive grid are selected as 500, 0.5, and 30. For MOABC, the population size is selected as 500, and for other parameter settings refer to [13].


### 3.4.4. Two Objective Functions’ Results

The results reported in terms of the best, worst, average, median, and standard deviation of the performance measures are listed in Tables1–4. Figures 4–7 show the optimal front obtained by four algorithms for each two objective functions, respectively. The continuous lines represent the Pareto-optimal front, while mark spots represent found nondominated solutions by the algorithms.Table 1
Comparison of performance on ZDT1.
ZDT1
CMOABC
MOABC
NSGA-II
MOPSO
Convergence metric
Average
8
.
4
93
2
e
-
0
0
4
1
.
3
87
9
e
-
0
0
2
2
.
3
01
2
e
-
00
1
1
.
6
89
2
e
-
00
3
Median
8
.
2
39
8
e
-
0
0
4
1
.
9
08
6
e
-
0
0
2
1
.
4
32
1
e
-
00
1
1
.
4
30
9
e
-
00
3
Min
7
.
5
65
6
e
-
0
0
4
9
.
8
86
7
e
-
0
0
3
7
.
2
36
2
e
-
00
2
8
.
4
79
8
e
-
00
4
Max
1
.
1
02
1
e
-
0
0
3
3
.
6
60
9
e
-
00
2
8
.
8
03
1
e
-
00
1
2
.
1
97
9
e
-
00
3
Std.
8
.
0
98
8
e
-
0
0
5
1
.
0
15
7
e
-
00
2
2
.
3
71
5
e
-
00
1
6
.
7
12
1
e
-
00
4
Diversity metric
Average
6
.
7
12
9
e
-
0
0
2
2
.
8
96
2
e
-
0
0
1
7
.
5
98
9
e
-
0
0
1
6
.
3
10
2
e
-
0
0
1
Median
6
.
2
48
3
e
-
0
0
1
3
.
3
66
3
e
-
0
0
1
5
.
1
98
6
e
-
0
0
1
6
.
8
96
3
e
-
0
0
1
Min
5
.
3
26
6
e
-
0
0
1
2
.
9
63
1
e
-
0
0
1
4
.
6
01
2
e
-
0
0
1
6
.
5
6
1
1
e
-
0
0
1
Max
6
.
2
30
9
e
-
0
0
1
3
.
7
35
4
e
-
0
0
1
9
.
6
12
3
e
-
0
0
1
7
.
2
02
2
e
-
0
0
1
Std.
4
.
4
15
6
e
-
0
0
2
2
.
4
21
2
e
-
0
0
2
1
.
5
30
1
e
-
0
0
1
2
.
4
26
2
e
-
0
0
2Table 2
Comparison of performance on ZDT2.
ZDT2
CMOABC
MOABC
NSGA-II
MOPSO
Convergence metric
Average
2
.
0
30
6
e
-
0
0
4
1
.
0
15
1
e
-
0
0
3
8
.
1
21
8
e
-
0
0
2
1
.
5
20
1
e
-
0
0
1
Median
3
.
3
72
7
e
_
00
4
9
.
4
90
1
e
-
0
0
4
9
.
5
23
1
e
-
0
0
2
7
.
9
93
1
e
-
0
0
2
Min
9
.
2
46
6
e
_
00
5
6
.
9
88
3
e
-
0
0
4
2
.
1
83
4
e
-
0
0
2
6
.
8
64
1
e
-
0
0
4
Max
1
.
8
27
6
e
-
0
0
3
1
.
0
92
4
e
-
0
0
3
1
.
1
63
2
e
-
0
0
1
6
.
1
07
3
e
-
0
0
1
Std.
6
.
0
91
9
e
-
0
0
4
2
.
5
06
1
e
-
0
0
4
3
.
3
87
6
e
-
0
0
2
6
.
5
66
3
e
-
0
0
1
Diversity metric
Average
6
.
4
83
2
e
-
0
0
2
2
.
8
95
7
e
-
0
0
1
6
.
0
09
7
e
-
0
0
1
6
.
2
90
7
e
-
0
0
1
Median
6
.
7
37
9
e
-
00
2
3
.
0
01
5
e
-
0
0
1
5
.
9
50
5
e
-
0
0
1
6
.
3
0
1
1
e
-
0
0
1
Min
6
.
1
05
8
e
-
00
2
2
.
7
63
1
e
-
0
0
1
5
.
4
03
8
e
-
0
0
1
5
.
1
87
6
e
-
0
0
1
Max
5
.
8
34
9
e
-
0
0
1
3
.
3
12
0
e
-
0
0
1
6
.
8
72
4
e
-
0
0
1
7
.
2
27
2
e
-
0
0
1
Std.
1
.
9
45
0
e
-
0
0
1
2
.
2
09
2
e
-
0
0
2
4
.
5
62
8
e
-
0
0
2
7
.
6
24
1
e
-
0
0
2Table 3
Comparison of performance on ZDT3.
ZDT3
CMOABC
MOABC
NSGA-II
MOPSO
Convergence metric
Average
4
.
3
25
6
e
-
0
0
4
1
.
0
01
2
e
-
0
0
3
2
.
9
47
6
e
+
0
0
0
1
.
6
42
1
e
-
0
0
2
Median
1
.
5
90
8
e
-
0
0
4
9
.
7
59
7
e
-
0
0
4
2
.
8
10
5
e
+
0
0
0
2
.
9
83
1
e
-
0
0
3
Min
9
.
3
66
4
e
-
0
0
5
8
.
9
18
7
e
-
0
0
4
5
.
0
66
8
e
-
0
0
1
1
.
8
27
7
e
-
0
0
3
Max
3
.
8
93
0
e
-
0
0
3
2
.
5
60
9
e
-
0
0
3
4
.
8
21
8
e
+
0
0
0
5
.
2
1
1
4
e
-
0
0
2
Std.
1
.
2
97
7
e
-
0
0
3
2
.
4
48
1
e
-
0
0
4
1
.
6
28
8
e
+
0
0
0
2
.
6
67
1
e
-
0
0
2
Diversity metric
Average
7
.
5
59
2
e
-
0
0
2
6
.
4
35
2
e
-
0
0
1
9
.
9
80
2
e
-
0
0
1
8
.
0
00
2
e
-
0
0
1
Median
6
.
9
82
5
e
-
0
0
2
6
.
4
18
5
e
-
0
0
1
8
.
9
97
9
e
-
0
0
1
6
.
9
87
1
e
-
0
0
1
Min
2
.
11
1
8
e
-
0
0
2
5
.
9
63
1
e
-
0
0
1
4
.
6
61
7
e
-
0
0
1
5
.
0
80
5
e
-
0
0
1
Max
6
.
8
03
2
e
-
0
0
1
6
.
8
11
1
e
-
0
0
1
1
.
0
54
2
e
+
0
0
0
8
.
9
39
2
e
-
0
0
1
Std.
2
.
2
67
7
e
-
0
0
1
1
.
2
14
7
e
-
0
0
2
2
.
3
56
7
e
-
0
0
2
3
.
3
28
8
e
-
0
0
2Table 4
Comparison of performance on ZDT6.
ZDT6
CMOABC
MOABC
NSGA-II
MOPSO
Convergence metric
Average
3
.
0
22
1
e
-
0
0
4
2
.
8
62
3
e
-
0
03
3
.
1
40
3
e
-
0
0
1
4
.
9
31
5
e
-
0
0
1
Median
3
.
6
05
9
e
-
0
0
4
9
.
2
56
7
e
-
0
0
4
2
.
4
21
4
e
-
0
0
1
8
.
6
06
1
e
-
0
0
4
Min
9
.
2
46
6
e
-
0
0
5
1
.
2
07
3
e
-
0
0
5
6
.
6
36
0
e
-
0
0
2
6
.
6
70
9
e
-
0
0
4
Max
2
.
7
19
9
e
-
0
0
3
1
.
9
54
2
e
-
0
0
2
8
.
7
22
0
e
-
0
0
1
6
.
2
1
1
4
e
-
0
0
2
Std.
9
.
0
66
2
e
-
0
0
4
7
.
8
52
6
e
-
0
0
3
2
.
8
45
2
e
-
0
0
1
6
.
6
57
1
e
-
0
0
3
Diversity metric
Average
6
.
8
40
7
e
-
0
0
2
5
.
9
75
2
e
-
0
0
1
7
.
5
90
5
e
-
0
0
1
6
.
5
90
2
e
-
0
0
1
Median
6
.
7
38
1
e
-
0
0
2
4
.
0
18
5
e
-
0
0
1
7
.
4
81
2
e
-
0
0
1
6
.
6
97
1
e
-
0
0
1
Min
6
.
2
05
9
e
-
0
0
2
3
.
9
85
1
e
-
0
0
1
5
.
7
13
5
e
-
0
0
1
5
.
8
89
5
e
-
0
0
1
Max
6
.
1
56
7
e
-
0
0
1
1
.
1
45
5
e
+
0
0
0
9
.
4
40
3
e
-
0
0
1
7
.
0
70
2
e
-
0
0
1
Std.
2
.
0
52
2
e
-
0
0
1
2
.
2
14
7
e
-
0
0
1
9
.
7
50
0
e
-
0
0
2
4
.
4
95
8
e
-
0
0
2Figure 4
Pareto fronts obtained by CMOABC, MOPSO, MOABC, and NSGA-II on ZDT1.Figure 5
Pareto fronts obtained by CMOABC, MOPSO, MOABC, and NSGA-II on ZDT2.Figure 6
Pareto fronts obtained by CMOABC, MOPSO, MOABC, and NSGA-II on ZDT3.Figure 7
Pareto fronts obtained by CMOABC, MOPSO, MOABC, and NSGA-II on ZDT6.On ZDT1 function, when given 10,000 function evaluations for four algorithms, Figure4 shows that CMOABC, MOPSO, NSGA-II, and MOABC have found a well distributed and diverse solution set in the entire Pareto-optimal region. However, NSGA-II and MOABC cannot achieve the true Pareto front for ZDT1. Table 1 shows that the performance of CMOABC in convergence metric is of three orders of magnitude better than NSGA-II, two orders better than MOABC, and one order better than MOPSO. For diversity metric, it can be observed that all algorithms have almost the same performance on ZDT1.On ZDT2 function, Figure5 shows that MOPSO and NSGA-II produce poor results on this test function and both of them cannot achieve the true Pareto front, and the results of the performance measures show that CMOABC and MOABC have better performance compared to the MOPSO and NSGA-II. From Table 3, it can be noticed that the performance of CMOABC and MOABC in convergence metric and diversity metric is better than that of MOPSO and NSGA-II. Moreover, CMOABC outperforms MOABC by one order of magnitude in terms of diversity metric.On ZDT3 function, from Figure6, it can be noticed that NSGA-II has difficulties in converging the true Pareto front, and other algorithms find a better converged set of nondominated solutions in ZDT3 compared to NSGA-II. Further, Table 3 also shows that NSGA-II is worse than CMOABC, MOPSO, and MOABC in convergence metric and diversity metric.On ZDT6 function, Figure7 and Table 4 show that the algorithms achieve similar performance as on ZDT3.


### 3.4.5. Three Objective Functions

The optimization results of CMOABC, MOPSO, MOABC, and NSGA-II algorithms for three objective DTLZ series problems are shown in Tables5 and 6. The Pareto-optimal region and all nondominated solutions obtained by the four algorithms are shown in Figures 8–11.Table 5
Comparison of performance on DTLZ2.
DTLZ2
CMOABC
MOABC
NSGA-II
MOPSO
Convergence metric
Average
3
.
0
89
5
e
-
0
0
4
8
.
0
02
3
e
-
00
3
8
.
7
89
6
e
-
0
0
2
3
.
8
68
4
e
-
0
0
3
Median
3
.
6
05
9
e
-
0
0
4
6
.
5
53
7
e
-
0
03
8
.
0
43
8
e
-
0
0
2
3
.
8
43
1
e
-
0
0
3
Min
9
.
8
02
1
e
-
0
0
5
5
.
2
18
3
e
-
0
03
6
.
1
40
9
e
-
0
0
2
3
.
7
63
1
e
-
0
0
3
Max
1
.
9
76
5
e
-
0
0
3
1
.
5
69
0
e
-
0
02
1
.
0
62
5
e
-
0
0
1
3
.
0
77
2
e
-
0
0
3
Std.
9
.
7
69
9
e
-
0
0
4
2
.
5
38
0
e
-
0
03
1
.
6
79
8
e
-
0
0
2
8
.
6
92
4
e
-
0
0
5
Diversity metric
Average
4
.
4
90
6
e
-
0
0
2
3
.
9
89
6
e
-
0
0
1
6
.
1
30
2
e
-
0
0
1
4
.
2
63
7
e
-
0
0
1
Median
6
.
7
38
1
e
-
0
0
2
3
.
9
68
5
e
-
0
0
1
5
.
2
99
9
e
-
0
0
1
4
.
1
52
1
e
-
0
0
1
Min
6
.
2
05
9
e
-
0
0
2
2
.
9
57
2
e
-
0
0
1
3
.
6
61
7
e
-
0
0
1
3
.
8
79
0
e
-
0
0
1
Max
1
.
4
90
6
e
-
0
0
1
5
.
2
56
1
e
-
0
0
1
7
.
2
54
2
e
-
0
0
1
4
.
5
61
8
e
-
0
0
1
Std.
1
.
8
64
9
e
-
0
0
2
3
.
2
23
4
e
-
0
0
2
1
.
7
56
7
e
-
0
0
2
2
.
1
69
2
e
-
0
0
2Table 6
Comparison of performance on DTLZ6.
DTLZ6
CMOABC
MOABC
NSGA-II
MOPSO
Convergence metric
Average
1
.
8
36
5
e
-
0
0
3
1
.
6
50
6
e
-
0
0
2
6
.
0
47
0
e
-
0
0
1
1
.
7
01
2
e
-
0
0
2
Median
3
.
6
21
9
e
-
0
0
4
3
.
4
04
1
e
-
0
0
3
3
.
2
73
3
e
-
0
0
1
1
.
8
00
2
e
-
0
0
2
Min
1
.
2
43
6
e
-
0
0
4
5
.
2
73
1
e
-
0
0
3
1
.
7
49
8
e
-
0
0
1
8
.
2
98
5
e
-
0
0
2
Max
9
.
1
82
6
e
-
0
0
3
2
.
8
04
3
e
-
0
0
2
2
.
9
72
8
e
+
0
0
0
2
.
8
96
1
e
-
0
0
2
Std.
4
.
1
06
6
e
-
0
0
3
4
.
9
12
7
e
-
0
0
3
8
.
4
99
7
e
-
0
0
1
8
.
3
34
6
e
-
0
0
3
Diversity metric
Average
9
.
0
41
3
e
-
0
0
2
2
.
9
02
1
e
-
0
0
1
5
.
6
47
7
e
-
0
0
1
5
.
5
00
2
e
-
0
0
1
Median
9
.
0
32
8
e
-
0
0
2
2
.
8
92
6
e
-
0
0
1
5
.
5
86
6
e
-
0
0
1
5
.
4
97
6
e
-
0
0
1
Min
7
.
8
77
2
e
-
0
0
2
9
.
8
26
6
e
-
0
0
2
5
.
1
84
9
e
-
0
0
1
5
.
2
01
5
e
-
0
0
1
Max
1
.
1
52
5
e
-
0
0
1
6
.
3
23
1
e
-
0
0
1
6
.
3
24
5
e
-
0
0
1
5
.
7
30
2
e
-
0
0
1
Std.
2
.
0
64
5
e
-
0
0
2
2
.
8
02
5
e
-
0
0
2
3
.
7
46
9
e
-
0
0
2
1
.
7
49
8
e
-
0
0
2Figure 8
The true Pareto-optimal front on DTLZ2.Figure 9
Pareto fronts obtained by CMOABC, MOPSO, MOABC, and NSGA-II on DTLZ2.Figure 10
The true Pareto-optimal front on DTLZ6.Figure 11
Pareto fronts obtained by CMOABC, MOPSO, MOABC, and NSGA-II on DTLZ6.From Figure9, it is observed that all algorithms obtain similar performance on DTLZ2. However, Table 5 shows the performance of NSGA-II is worse than that of the other three algorithms in convergence, while the performance of CMOABC is better than that of the other three algorithms in diversity.From Figure11, it can be seen that CMOABC and MOABC have better performance on DTLZ6. The CMOABC especially has attained the nondominated set properly. The front obtained from the other two algorithms cannot be uniformly distributed on the true Pareto-optimal front on DTLZ6.Table6 also shows that MOABC has obtained better convergence and diversity than NSGA-II and MOPSO on DTLZ6, but it is found to be second to the CMOABC algorithm.On the whole, for two and three objective problems, CMOABC has a more competitive performance than MOPSO, MOABC, and NSGA-II on both diversity metric and convergence metric. This is due to the fact that, in the proposed CMOABC algorithm, the modified employed bees’ phase with a self-adaptive strength of information communication effectively enhances its abilities of exploration in the early stage and exploitation in the late stage. Moreover, benefiting from the swarm decomposition operation with a dynamicK-means clustering method, CMOABC could ensure more sufficient exploitation around the potential solutions in the less crowded regions. When the number of clusters is large, the probability of the individuals with similar performances allocated into the same subcomponent increases, which improves the ability of global search. When the number of clusters becomes small, the individuals trapped into local optimal may move toward the current best individual of their cluster, which effectively avoids the inefficiency of local search. In addition, the regeneration operation also enriches the diversity of external archive. Therefore, CMOABC can reasonably accommodate a considerable potential for solving more complex multiobjective optimization problems.


## 3.4.1. Test Function

To fully evaluate the performance of the CMOABC algorithm without a biased conclusion towards some chosen problems, in this paper, six commonly recognized benchmark functions have been used. Four of these test problems ZDT1, ZDT2, ZDT3, and ZDT6 [32] are of two objectives, while the other two, DLTZ2 and DLTZ6 [33], are of three objectives.


## 3.4.2. Performance Measures

In order to facilitate the quantitative assessment of the performance of a multiobjective optimization algorithm, two performance metrics are taken into consideration: (1) convergence metricλ and (2) diversity metric η [4].Convergence Metric. This metric measures the extent of convergence to a known set of Pareto-optimal solutions:(30)λ=∑i=1NdiN,where N is the number of nondominated solutions obtained with an algorithm and di is the Euclidean distance between each of the nondominated solutions and the nearest member of the true Pareto-optimal front. To calculate this metric, we find a set of H=500 uniformly spaced solutions from the true Pareto-optimal front in the objective space. For each solution obtained with an algorithm, we compute the minimum Euclidean distance of it from H chosen solutions on the Pareto-optimal front. The average of these distances is used as the convergence metric λ. Figure 2 shows the calculation procedure of this metric.Figure 2
Convergence metricλ.Diversity Metric. This metric measures the extent of spread achieved among the obtained solutions. Here, we are interested in getting a set of solutions that spans the entire Pareto-optimal region. This metric is defined as(31)η=df+dl+∑i=1N-1di-d¯df+dl+N-1d¯,where di is the Euclidean distance between consecutive solutions in the obtained nondominated set of solutions and N is the number of nondominated solutions obtained by an algorithm. d¯ is the average value of these distances. df and dl are the Euclidean distances between the extreme solutions and the boundary solutions of the obtained nondominated set, as depicted in Figure 3.Figure 3
Diversity metricη.


## 3.4.3. Experimental Setting

The effectiveness of CMOABC will be demonstrated on various standard test problems in this section. In addition, the proposed algorithm will be compared with the nondominated sorting genetic algorithm II (NSGA-II) [4], the multiobjective particle swarm optimization (MOPSO) [47], and the multiobjective artificial bee colony algorithm (MOABC) [13]. In order to compare the different algorithms with a fair time measure, the number of function evaluations (Fes = 10,000) is used for the termination criterion.For two test examples presented in this study, the CMOABC algorithm parameters are set as follows: the population size is 500 and the setG is defined as {50,20,10,5}.For NSGA-II, the population size, crossover, and mutation probabilities are selected as 500, 0.8, and 0.2, respectively, for the two test examples. For MOPSO, the population size, mutation rate, and divisions for the adaptive grid are selected as 500, 0.5, and 30. For MOABC, the population size is selected as 500, and for other parameter settings refer to [13].


## 3.4.4. Two Objective Functions’ Results

The results reported in terms of the best, worst, average, median, and standard deviation of the performance measures are listed in Tables1–4. Figures 4–7 show the optimal front obtained by four algorithms for each two objective functions, respectively. The continuous lines represent the Pareto-optimal front, while mark spots represent found nondominated solutions by the algorithms.Table 1
Comparison of performance on ZDT1.
ZDT1
CMOABC
MOABC
NSGA-II
MOPSO
Convergence metric
Average
8
.
4
93
2
e
-
0
0
4
1
.
3
87
9
e
-
0
0
2
2
.
3
01
2
e
-
00
1
1
.
6
89
2
e
-
00
3
Median
8
.
2
39
8
e
-
0
0
4
1
.
9
08
6
e
-
0
0
2
1
.
4
32
1
e
-
00
1
1
.
4
30
9
e
-
00
3
Min
7
.
5
65
6
e
-
0
0
4
9
.
8
86
7
e
-
0
0
3
7
.
2
36
2
e
-
00
2
8
.
4
79
8
e
-
00
4
Max
1
.
1
02
1
e
-
0
0
3
3
.
6
60
9
e
-
00
2
8
.
8
03
1
e
-
00
1
2
.
1
97
9
e
-
00
3
Std.
8
.
0
98
8
e
-
0
0
5
1
.
0
15
7
e
-
00
2
2
.
3
71
5
e
-
00
1
6
.
7
12
1
e
-
00
4
Diversity metric
Average
6
.
7
12
9
e
-
0
0
2
2
.
8
96
2
e
-
0
0
1
7
.
5
98
9
e
-
0
0
1
6
.
3
10
2
e
-
0
0
1
Median
6
.
2
48
3
e
-
0
0
1
3
.
3
66
3
e
-
0
0
1
5
.
1
98
6
e
-
0
0
1
6
.
8
96
3
e
-
0
0
1
Min
5
.
3
26
6
e
-
0
0
1
2
.
9
63
1
e
-
0
0
1
4
.
6
01
2
e
-
0
0
1
6
.
5
6
1
1
e
-
0
0
1
Max
6
.
2
30
9
e
-
0
0
1
3
.
7
35
4
e
-
0
0
1
9
.
6
12
3
e
-
0
0
1
7
.
2
02
2
e
-
0
0
1
Std.
4
.
4
15
6
e
-
0
0
2
2
.
4
21
2
e
-
0
0
2
1
.
5
30
1
e
-
0
0
1
2
.
4
26
2
e
-
0
0
2Table 2
Comparison of performance on ZDT2.
ZDT2
CMOABC
MOABC
NSGA-II
MOPSO
Convergence metric
Average
2
.
0
30
6
e
-
0
0
4
1
.
0
15
1
e
-
0
0
3
8
.
1
21
8
e
-
0
0
2
1
.
5
20
1
e
-
0
0
1
Median
3
.
3
72
7
e
_
00
4
9
.
4
90
1
e
-
0
0
4
9
.
5
23
1
e
-
0
0
2
7
.
9
93
1
e
-
0
0
2
Min
9
.
2
46
6
e
_
00
5
6
.
9
88
3
e
-
0
0
4
2
.
1
83
4
e
-
0
0
2
6
.
8
64
1
e
-
0
0
4
Max
1
.
8
27
6
e
-
0
0
3
1
.
0
92
4
e
-
0
0
3
1
.
1
63
2
e
-
0
0
1
6
.
1
07
3
e
-
0
0
1
Std.
6
.
0
91
9
e
-
0
0
4
2
.
5
06
1
e
-
0
0
4
3
.
3
87
6
e
-
0
0
2
6
.
5
66
3
e
-
0
0
1
Diversity metric
Average
6
.
4
83
2
e
-
0
0
2
2
.
8
95
7
e
-
0
0
1
6
.
0
09
7
e
-
0
0
1
6
.
2
90
7
e
-
0
0
1
Median
6
.
7
37
9
e
-
00
2
3
.
0
01
5
e
-
0
0
1
5
.
9
50
5
e
-
0
0
1
6
.
3
0
1
1
e
-
0
0
1
Min
6
.
1
05
8
e
-
00
2
2
.
7
63
1
e
-
0
0
1
5
.
4
03
8
e
-
0
0
1
5
.
1
87
6
e
-
0
0
1
Max
5
.
8
34
9
e
-
0
0
1
3
.
3
12
0
e
-
0
0
1
6
.
8
72
4
e
-
0
0
1
7
.
2
27
2
e
-
0
0
1
Std.
1
.
9
45
0
e
-
0
0
1
2
.
2
09
2
e
-
0
0
2
4
.
5
62
8
e
-
0
0
2
7
.
6
24
1
e
-
0
0
2Table 3
Comparison of performance on ZDT3.
ZDT3
CMOABC
MOABC
NSGA-II
MOPSO
Convergence metric
Average
4
.
3
25
6
e
-
0
0
4
1
.
0
01
2
e
-
0
0
3
2
.
9
47
6
e
+
0
0
0
1
.
6
42
1
e
-
0
0
2
Median
1
.
5
90
8
e
-
0
0
4
9
.
7
59
7
e
-
0
0
4
2
.
8
10
5
e
+
0
0
0
2
.
9
83
1
e
-
0
0
3
Min
9
.
3
66
4
e
-
0
0
5
8
.
9
18
7
e
-
0
0
4
5
.
0
66
8
e
-
0
0
1
1
.
8
27
7
e
-
0
0
3
Max
3
.
8
93
0
e
-
0
0
3
2
.
5
60
9
e
-
0
0
3
4
.
8
21
8
e
+
0
0
0
5
.
2
1
1
4
e
-
0
0
2
Std.
1
.
2
97
7
e
-
0
0
3
2
.
4
48
1
e
-
0
0
4
1
.
6
28
8
e
+
0
0
0
2
.
6
67
1
e
-
0
0
2
Diversity metric
Average
7
.
5
59
2
e
-
0
0
2
6
.
4
35
2
e
-
0
0
1
9
.
9
80
2
e
-
0
0
1
8
.
0
00
2
e
-
0
0
1
Median
6
.
9
82
5
e
-
0
0
2
6
.
4
18
5
e
-
0
0
1
8
.
9
97
9
e
-
0
0
1
6
.
9
87
1
e
-
0
0
1
Min
2
.
11
1
8
e
-
0
0
2
5
.
9
63
1
e
-
0
0
1
4
.
6
61
7
e
-
0
0
1
5
.
0
80
5
e
-
0
0
1
Max
6
.
8
03
2
e
-
0
0
1
6
.
8
11
1
e
-
0
0
1
1
.
0
54
2
e
+
0
0
0
8
.
9
39
2
e
-
0
0
1
Std.
2
.
2
67
7
e
-
0
0
1
1
.
2
14
7
e
-
0
0
2
2
.
3
56
7
e
-
0
0
2
3
.
3
28
8
e
-
0
0
2Table 4
Comparison of performance on ZDT6.
ZDT6
CMOABC
MOABC
NSGA-II
MOPSO
Convergence metric
Average
3
.
0
22
1
e
-
0
0
4
2
.
8
62
3
e
-
0
03
3
.
1
40
3
e
-
0
0
1
4
.
9
31
5
e
-
0
0
1
Median
3
.
6
05
9
e
-
0
0
4
9
.
2
56
7
e
-
0
0
4
2
.
4
21
4
e
-
0
0
1
8
.
6
06
1
e
-
0
0
4
Min
9
.
2
46
6
e
-
0
0
5
1
.
2
07
3
e
-
0
0
5
6
.
6
36
0
e
-
0
0
2
6
.
6
70
9
e
-
0
0
4
Max
2
.
7
19
9
e
-
0
0
3
1
.
9
54
2
e
-
0
0
2
8
.
7
22
0
e
-
0
0
1
6
.
2
1
1
4
e
-
0
0
2
Std.
9
.
0
66
2
e
-
0
0
4
7
.
8
52
6
e
-
0
0
3
2
.
8
45
2
e
-
0
0
1
6
.
6
57
1
e
-
0
0
3
Diversity metric
Average
6
.
8
40
7
e
-
0
0
2
5
.
9
75
2
e
-
0
0
1
7
.
5
90
5
e
-
0
0
1
6
.
5
90
2
e
-
0
0
1
Median
6
.
7
38
1
e
-
0
0
2
4
.
0
18
5
e
-
0
0
1
7
.
4
81
2
e
-
0
0
1
6
.
6
97
1
e
-
0
0
1
Min
6
.
2
05
9
e
-
0
0
2
3
.
9
85
1
e
-
0
0
1
5
.
7
13
5
e
-
0
0
1
5
.
8
89
5
e
-
0
0
1
Max
6
.
1
56
7
e
-
0
0
1
1
.
1
45
5
e
+
0
0
0
9
.
4
40
3
e
-
0
0
1
7
.
0
70
2
e
-
0
0
1
Std.
2
.
0
52
2
e
-
0
0
1
2
.
2
14
7
e
-
0
0
1
9
.
7
50
0
e
-
0
0
2
4
.
4
95
8
e
-
0
0
2Figure 4
Pareto fronts obtained by CMOABC, MOPSO, MOABC, and NSGA-II on ZDT1.Figure 5
Pareto fronts obtained by CMOABC, MOPSO, MOABC, and NSGA-II on ZDT2.Figure 6
Pareto fronts obtained by CMOABC, MOPSO, MOABC, and NSGA-II on ZDT3.Figure 7
Pareto fronts obtained by CMOABC, MOPSO, MOABC, and NSGA-II on ZDT6.On ZDT1 function, when given 10,000 function evaluations for four algorithms, Figure4 shows that CMOABC, MOPSO, NSGA-II, and MOABC have found a well distributed and diverse solution set in the entire Pareto-optimal region. However, NSGA-II and MOABC cannot achieve the true Pareto front for ZDT1. Table 1 shows that the performance of CMOABC in convergence metric is of three orders of magnitude better than NSGA-II, two orders better than MOABC, and one order better than MOPSO. For diversity metric, it can be observed that all algorithms have almost the same performance on ZDT1.On ZDT2 function, Figure5 shows that MOPSO and NSGA-II produce poor results on this test function and both of them cannot achieve the true Pareto front, and the results of the performance measures show that CMOABC and MOABC have better performance compared to the MOPSO and NSGA-II. From Table 3, it can be noticed that the performance of CMOABC and MOABC in convergence metric and diversity metric is better than that of MOPSO and NSGA-II. Moreover, CMOABC outperforms MOABC by one order of magnitude in terms of diversity metric.On ZDT3 function, from Figure6, it can be noticed that NSGA-II has difficulties in converging the true Pareto front, and other algorithms find a better converged set of nondominated solutions in ZDT3 compared to NSGA-II. Further, Table 3 also shows that NSGA-II is worse than CMOABC, MOPSO, and MOABC in convergence metric and diversity metric.On ZDT6 function, Figure7 and Table 4 show that the algorithms achieve similar performance as on ZDT3.


## 3.4.5. Three Objective Functions

The optimization results of CMOABC, MOPSO, MOABC, and NSGA-II algorithms for three objective DTLZ series problems are shown in Tables5 and 6. The Pareto-optimal region and all nondominated solutions obtained by the four algorithms are shown in Figures 8–11.Table 5
Comparison of performance on DTLZ2.
DTLZ2
CMOABC
MOABC
NSGA-II
MOPSO
Convergence metric
Average
3
.
0
89
5
e
-
0
0
4
8
.
0
02
3
e
-
00
3
8
.
7
89
6
e
-
0
0
2
3
.
8
68
4
e
-
0
0
3
Median
3
.
6
05
9
e
-
0
0
4
6
.
5
53
7
e
-
0
03
8
.
0
43
8
e
-
0
0
2
3
.
8
43
1
e
-
0
0
3
Min
9
.
8
02
1
e
-
0
0
5
5
.
2
18
3
e
-
0
03
6
.
1
40
9
e
-
0
0
2
3
.
7
63
1
e
-
0
0
3
Max
1
.
9
76
5
e
-
0
0
3
1
.
5
69
0
e
-
0
02
1
.
0
62
5
e
-
0
0
1
3
.
0
77
2
e
-
0
0
3
Std.
9
.
7
69
9
e
-
0
0
4
2
.
5
38
0
e
-
0
03
1
.
6
79
8
e
-
0
0
2
8
.
6
92
4
e
-
0
0
5
Diversity metric
Average
4
.
4
90
6
e
-
0
0
2
3
.
9
89
6
e
-
0
0
1
6
.
1
30
2
e
-
0
0
1
4
.
2
63
7
e
-
0
0
1
Median
6
.
7
38
1
e
-
0
0
2
3
.
9
68
5
e
-
0
0
1
5
.
2
99
9
e
-
0
0
1
4
.
1
52
1
e
-
0
0
1
Min
6
.
2
05
9
e
-
0
0
2
2
.
9
57
2
e
-
0
0
1
3
.
6
61
7
e
-
0
0
1
3
.
8
79
0
e
-
0
0
1
Max
1
.
4
90
6
e
-
0
0
1
5
.
2
56
1
e
-
0
0
1
7
.
2
54
2
e
-
0
0
1
4
.
5
61
8
e
-
0
0
1
Std.
1
.
8
64
9
e
-
0
0
2
3
.
2
23
4
e
-
0
0
2
1
.
7
56
7
e
-
0
0
2
2
.
1
69
2
e
-
0
0
2Table 6
Comparison of performance on DTLZ6.
DTLZ6
CMOABC
MOABC
NSGA-II
MOPSO
Convergence metric
Average
1
.
8
36
5
e
-
0
0
3
1
.
6
50
6
e
-
0
0
2
6
.
0
47
0
e
-
0
0
1
1
.
7
01
2
e
-
0
0
2
Median
3
.
6
21
9
e
-
0
0
4
3
.
4
04
1
e
-
0
0
3
3
.
2
73
3
e
-
0
0
1
1
.
8
00
2
e
-
0
0
2
Min
1
.
2
43
6
e
-
0
0
4
5
.
2
73
1
e
-
0
0
3
1
.
7
49
8
e
-
0
0
1
8
.
2
98
5
e
-
0
0
2
Max
9
.
1
82
6
e
-
0
0
3
2
.
8
04
3
e
-
0
0
2
2
.
9
72
8
e
+
0
0
0
2
.
8
96
1
e
-
0
0
2
Std.
4
.
1
06
6
e
-
0
0
3
4
.
9
12
7
e
-
0
0
3
8
.
4
99
7
e
-
0
0
1
8
.
3
34
6
e
-
0
0
3
Diversity metric
Average
9
.
0
41
3
e
-
0
0
2
2
.
9
02
1
e
-
0
0
1
5
.
6
47
7
e
-
0
0
1
5
.
5
00
2
e
-
0
0
1
Median
9
.
0
32
8
e
-
0
0
2
2
.
8
92
6
e
-
0
0
1
5
.
5
86
6
e
-
0
0
1
5
.
4
97
6
e
-
0
0
1
Min
7
.
8
77
2
e
-
0
0
2
9
.
8
26
6
e
-
0
0
2
5
.
1
84
9
e
-
0
0
1
5
.
2
01
5
e
-
0
0
1
Max
1
.
1
52
5
e
-
0
0
1
6
.
3
23
1
e
-
0
0
1
6
.
3
24
5
e
-
0
0
1
5
.
7
30
2
e
-
0
0
1
Std.
2
.
0
64
5
e
-
0
0
2
2
.
8
02
5
e
-
0
0
2
3
.
7
46
9
e
-
0
0
2
1
.
7
49
8
e
-
0
0
2Figure 8
The true Pareto-optimal front on DTLZ2.Figure 9
Pareto fronts obtained by CMOABC, MOPSO, MOABC, and NSGA-II on DTLZ2.Figure 10
The true Pareto-optimal front on DTLZ6.Figure 11
Pareto fronts obtained by CMOABC, MOPSO, MOABC, and NSGA-II on DTLZ6.From Figure9, it is observed that all algorithms obtain similar performance on DTLZ2. However, Table 5 shows the performance of NSGA-II is worse than that of the other three algorithms in convergence, while the performance of CMOABC is better than that of the other three algorithms in diversity.From Figure11, it can be seen that CMOABC and MOABC have better performance on DTLZ6. The CMOABC especially has attained the nondominated set properly. The front obtained from the other two algorithms cannot be uniformly distributed on the true Pareto-optimal front on DTLZ6.Table6 also shows that MOABC has obtained better convergence and diversity than NSGA-II and MOPSO on DTLZ6, but it is found to be second to the CMOABC algorithm.On the whole, for two and three objective problems, CMOABC has a more competitive performance than MOPSO, MOABC, and NSGA-II on both diversity metric and convergence metric. This is due to the fact that, in the proposed CMOABC algorithm, the modified employed bees’ phase with a self-adaptive strength of information communication effectively enhances its abilities of exploration in the early stage and exploitation in the late stage. Moreover, benefiting from the swarm decomposition operation with a dynamicK-means clustering method, CMOABC could ensure more sufficient exploitation around the potential solutions in the less crowded regions. When the number of clusters is large, the probability of the individuals with similar performances allocated into the same subcomponent increases, which improves the ability of global search. When the number of clusters becomes small, the individuals trapped into local optimal may move toward the current best individual of their cluster, which effectively avoids the inefficiency of local search. In addition, the regeneration operation also enriches the diversity of external archive. Therefore, CMOABC can reasonably accommodate a considerable potential for solving more complex multiobjective optimization problems.


## 4. Multiobjective Optimal Power Flow Based on CMOABC

In order to prove the numerical correctness, efficiency, and validation of CMOABC algorithm for the model, in this section, the multiobjective OPF problem will be used to implement simulation experiment. The proposed algorithm will be applied to two objectives and three objectives of OPF problem and it is compared with NSGA-II, MOPSO, and MOABC to solve such problem.


### 4.1. Implementation of the CMOABC Algorithm for the OPF Problem

Step 1.
Input the system parameters and the minimum and maximum limits of control variables.Step 2.
Input parameters of the CMOABC algorithm and the lower and upper limits of each variable.Step 3.
Produce the initial population.Step 4.
Cluster the population byK-means and compute the fitness.Step 5.
Store nondominated solutions in the memory table.Step 6.
Optimize every cluster individually.Step 7.
Handle the constraints.Step 8.
Compute the fitness of the solutions and store nondominated solutions into EA.Step 9.
Update each EA of each hive according to greedy selecting strategy, sort the EA based on nondomination, and select the nondomination solutions to stay in EA. If the number of nondominated solutions exceeds the allocated size of EA, apply crowding distance to remove the crowded members.Step 10.
Judge the condition of reclustering. If it meets the condition, the number of clusters can be reset and go to Step4; if not, go to Step 6.


### 4.2. Best Compromise Solution Using Fuzzy Membership Approach

In this work, a fuzzy-based mechanism is employed to extract the best compromise solution over the trade-off curve and assist the decision maker to adjust the generation levels efficiently [48]. Due to imprecise nature of the decision maker’s judgment, each objective function of the ith solution is represented by a membership function μi defined as follows:(32)μi=1Fi≤min⁡⁡Fi,max⁡⁡Fi-Fimax⁡⁡Fi-min⁡⁡Fimin⁡⁡Fi≤Fi≤max⁡⁡Fi0Fi≥max⁡⁡Fi,where min⁡⁡(Fi) and max⁡⁡(Fi) are lower and upper bounds of the ith objective function. The higher the values of the membership function are, the greater the solution satisfaction is.For each nondominated solution, the normalized membership functionμk is calculated as(33)μk=∑i=1Nobjμik∑k=1M∑I=1Nobjμik,where M is the number of nondominated solutions and Nobj is the number of objects. The best compromise solution is the one having the maximum of μk.


### 4.3. Simulation Results

In order to validate the robustness of the proposed CMOABC method, a standard IEEE 30-bus system has been used as the test system [28]. The system represents a portion of the American Electric Power System (in the Midwestern USA) as of December 1961, which consists of six generators, 41 transmission lines, and 4 transformers with off-nominal tap ratio in lines 6–9, 6–10, 4–12, and 27-28. The system data with minimum and maximum limits of control variables are shown in Table 7 [49]. The limits of generator buses and load buses are between 0.95 and 1.1 p.u. and 0.9 and 1.05 p.u., respectively. The lower and upper limits of transformer taps are 0.9 p.u. and 1.05 p.u., respectively, and the step size is 0.01 p.u. The parameter settings for these four algorithms are the same as in “performance setting” in Section 3.4.3. In the following section, there are two cases of the test system considered.(1)
Two objectives: emission cost, emission loss, and loss cost are considered, respectively.(2)
Three objectives: they are emission, cost, and loss.Table 7
Characteristics of the generation units.
G1
G2
G3
G4
G5
G6
Generator limits
P
C
max
⁡ (MW)
150
150
150
150
150
150
P
C
min
⁡ (MW)
5
5
5
5
5
5
Cost coefficients
a
10
10
20
10
20
10
b
200
150
180
100
180
150
c
100
120
40
60
40
100
Emission coefficients
α
4
.
0
9
1
2
.
5
4
3
4
.
2
5
8
5
.
3
2
6
4
.
2
5
8
6
.
1
3
1
β
-
5
.
55
4
-
6
.
40
7
-
5
.
09
4
-
3
.
55
0
-
5
.
09
4
-
5
.
55
5
γ
6
.
4
9
0
5
.
6
3
8
4
.
5
8
6
3
.
3
8
0
4
.
5
8
6
5
.
1
5
1
δ
2
.
0
e
-
4
5
.
0
e
-
4
1
.
0
e
-
6
2
.
0
e
-
3
1
.
0
e
-
6
1
.
0
e
-
5
ε
2
.
8
5
7
3
.
3
3
3
8
.
0
0
0
2
.
0
0
0
8
.
0
0
0
6
.
6
6
7


#### 4.3.1. Case I: Two-Objective OPF Optimization

In this case, the results obtained from the proposed CMOABC have been compared with MOPSO, MOABC, and NSGA-II in three cases: minimum cost, minimum emission, and minimum loss. Figures12–14 show the convergence characteristics of the four algorithms on three two-objective OPF problems, respectively. For the case of emission cost, Table 8 gives the best solutions for minimum cost and minimum emission; the data of emission loss and loss cost are shown in Tables 10 and 12. Finally, the best compromise solutions for each objective in the two-dimensional Pareto front are shown in Tables 9, 11, and 13.Table 8
The best solutions for cost and emission from the Pareto front based on four multiobjective algorithms.
CMOABC
MOABC
NSGA-II
MOPSO
Bestf1
Bestf2
Bestf1
Bestf2
Bestf1
Bestf2
Bestf1
Bestf2
PG1
12.36
41.30
11.09
7.92
36.25
42.01
26.11
56.01
PG2
29.99
45.41
30.16
42.08
55.45
47.91
32.99
46.62
PG3
55.40
55.31
57.50
60.89
49.60
50.97
62.14
56.12
PG4
105.16
43.40
105.41
101.06
46.20
43.10
102.13
40.31
PG5
45.26
52.79
45.39
42.61
56.92
52.20
45.21
56.03
PG6
38.01
52.01
37.81
33.88
44.32
52.06
38.15
48.96
f
1
(cost)
606.52
639.09
607.18
645.01
610.80
639.89
607.12
643.90
f
2 (emission)
0.2255
0.1931
0.2259
0.2009
0.2317
0.1992
0.2590
0.1959Table 9
The best compromise solutions for emission and cost using different multiobjective algorithms.
CMOABC
MOABC
NSGA-II
MOPSO
PG1
24.4905
24.9981
36.2012
24.3132
PG2
38.9856
39.5037
55.4215
37.8013
PG3
56.6471
39.0035
50.1034
55.4107
PG4
72.0001
73.9851
45.9804
72.1045
PG5
51.0322
48.1023
56.7635
52.0511
PG6
45.0258
42.9861
44.2014
43.9648
f
1 (cost)
616.7902
617.4896
627.0113
617.5842
f
2 (emission)
0.2015
0.2033
0.1968
0.2019Table 10
The best solutions for cost and loss from the Pareto front based on four multiobjective algorithms.
CMOABC
MOABC
NSGA-II
MOPSO
Bestf1
Bestf3
Bestf1
Bestf3
Bestf1
Bestf3
Bestf1
Bestf3
PG1
15.13
3.69
9.62
4.02
19.06
3.85
66.98
9.54
PG2
29.92
27.41
32.16
11.24
20.65
23.66
37.99
10.02
PG3
55.89
105.93
57.99
109.92
39.62
72.49
68.16
109.34
PG4
112.66
57.01
95.41
67.53
118.69
99.04
102.03
69.36
PG5
31.23
5.23
44.38
11.64
42.02
5.27
38.21
71.81
PG6
38.95
87.56
39.85
80.63
47.14
82.31
28.15
33.65
f
1 (cost)
607.26
626.09
609.99
626.12
611.01
643.26
609.42
669.91
f
3 (loss)
2.4514
1.5196
2.4643
1.5489
3.8167
1.5841
2.5389
1.5826Table 11
The best compromise solutions for loss and cost using different multiobjective algorithms.
CMOABC
MOABC
NSGA-II
MOPSO
PG1
4.4112
8.6014
3.6012
11.6012
PG2
25.0134
29.1642
12.4035
20.7035
PG3
79.4701
74.8653
87.6634
83.0107
PG4
106.2073
105.2366
91.9804
96.1342
PG5
8.0312
17.5038
28.0135
12.0511
PG6
62.0252
50.0031
61.8014
70.9652
f
1 (cost)
622.6002
613.3213
623.2813
630.5862
f
3 (loss)
1.6691
1.7097
1.8703
1.6258Table 12
The best solutions for emission and loss from the Pareto front based on four multiobjective algorithms.
CMOABC
MOABC
NSGA-II
MOPSO
Bestf2
Bestf3
Bestf2
Bestf3
Bestf2
Bestf3
Bestf2
Bestf3
PG1
45.13
9.09
33.90
12.06
44.01
63.30
33.01
10.97
PG2
44.92
5.02
60.03
6.48
41.86
144.51
44.69
77.01
PG3
55.96
91.01
49.89
70.01
57.29
108.58
61.98
23.11
PG4
41.02
70.98
34.65
95.32
40.66
62.03
43.03
63.35
PG5
46.63
5.13
63.61
6.39
42.61
5.12
49.96
82.49
PG6
54.02
103.86
46.96
95.22
34.59
34.59
52.59
58.60
f
2 (emission)
0.1958
0.2489
0.1969
0.2519
0.2501
0.2334
0.1962
0.2599
f
3 (loss)
4.1213
1.6072
4.7821
1.6501
4.2996
2.1239
3.8300
1.6640Table 13
The best compromise solutions for loss and emission using different multiobjective algorithms.
CMOABC
MOABC
NSGA-II
MOPSO
PG1
21.1132
68.9014
30.6011
22.6312
PG2
30.9814
23.0042
32.1035
25.1034
PG3
79.4721
30.8973
77.6934
103.0112
PG4
66.2076
65.9956
65.9894
64.0343
PG5
18.0382
77.5468
34.0243
6.4802
PG6
72.0258
19.0032
48.0014
85.9639
f
2 (emission)
2.0065
2.0143
2.6470
1.5959
f
3 (loss)
0.2035
0.2132
0.2055
0.2439Figure 12
Pareto fronts obtained by CMOABC, MOPSO, MOABC, and NSGA-II for emission and cost.Figure 13
Pareto fronts obtained by CMOABC, MOPSO, MOABC, and NSGA-II for cost and loss.Figure 14
Pareto fronts obtained by CMOABC, MOPSO, MOABC, and NSGA-II for emission and loss.First, two competing objectives, fuel cost and emission, are considered. From the Pareto-optimal solution shown in Figure12, it is clear that the proposed CMOABC method is able to give well distributed solutions compared to the other three algorithms. From Table 8, it is clear that the minimum fuel cost obtained by CMOABC is 606.52 $/h, which is better than the results optimized by the other three algorithms. And when the best fuel cost effect case is considered, CMOABC method can provide smaller loss. The best compromise solution is 616.7902 $/h and 0.2015 ton/h of CMOABC which can be observed in Table 9. It can once again prove that the proposed method has obtained better solutions.The cost loss objective functions are shown in Figure13 and Table 10. From the figure, we can see that the proposed algorithm obtains better results than the other three algorithms, though the obtained Pareto fronts are not evenly distributed. From Table 10, it is easily observed that CMOABC method can provide lower economic fuel cost when the minimum loss case is considered and smaller loss when the best fuel cost effect case is considered.For emission loss objective functions, we can observe that the algorithms achieve similar performance ranking as for fuel cost emission objective functions.


#### 4.3.2. Case II: Three-Objective OPF Optimization

In this case, the three objectives are optimized simultaneously by the four algorithms and the corresponding best solutions are given in Table14. Figure 15 also shows the result values of three competing objectives.Table 14
The best solutions for cost, emission, and loss based on four multiobjective algorithms.
CMOABC
MOABC
MOPSO
NSGA-II
Best
Best
Best
Best
Best
Best
Best
Best
Best
Best
Best
Best
f
1
f
2
f
3
f
1
f
2
f
3
f
1
f
2
f
3
f
1
f
2
f
3
PG1
20.02
44.96
6.01
35.66
3.21
18.95
22.01
36.97
9.01
25.86
36.63
18.65
PG2
21.03
45.01
15.26
49.25
16.74
30.01
25.63
45.76
25.75
28.61
49.87
27.83
PG3
70.25
59.79
106.92
62.73
81.51
61.36
64.59
61.02
83.83
58.02
50.76
95.62
PG4
101.05
38.12
68.95
39.03
69.01
102.65
101.81
45.53
64.92
48.72
38.52
66.75
PG5
46.12
52.68
5.00
45.67
8.96
36.57
40.92
49.89
68.66
43.15
51.34
7.97
PG6
35.19
46.11
87.65
33.52
76.72
39.96
32.56
51.64
55.65
42.84
49.42
68.74
f
1 (cost)
612.01
649.86
639.96
651.02
626.84
613.02
616.98
645.09
663.72
618.02
651.64
645.52
f
2 (emission)
0.2498
0.1956
0.2398
0.2069
0.2497
0.2501
0.3805
0.2101
0.3964
0.2426
0.2193
0.2467
f
3 (loss)
1.9968
4.0117
1.6046
3.9964
1.7013
1.9904
2.0972
2.4657
1.5206
2.3994
2.3652
2.0985Figure 15
Pareto fronts obtained by CMOABC, MOPSO, MOABC, and NSGA-II for fuel cost, emission, and loss.
(a)
CMOABC
(b)
MOABC
(c)
MOPSO
(d)
NSGA-IIAs seen in Figure15, compared to the other three algorithms, it can be observed that the Pareto-optimal solutions obtained by the proposed CMOABC are better distributed on the front with good diversity. Among the other three algorithms, the Pareto-optimal solutions obtained by the standard MOABC are also well distributed and the diversification of them is not as good as the ones obtained by the proposed CMOABC.Furthermore, from Table14, it can be also concluded that the CMOABC is able to discover a well distributed and diverse solution set for this three-objective problem. However, the other three algorithms cannot archive the true Pareto front for the three-objective OPF problem.The compromise solutions are shown in Table15. From this table we can see the proposed algorithm obtains the best emission, cost, though the system loss obtained by CMOABC is 2.1601 MW, which is a little more than 2.1598 MW obtained by MOABC.Table 15
The best compromise solutions for cost, emission, and loss using different multiobjective algorithms.
CMOABC
MOABC
NSGA-II
MOPSO
PG1
19.0732
18.9234
35.6214
22.0152
PG2
32.9875
27.9542
53.9065
15.1023
PG3
68.0549
70.8965
47.6936
90.0136
PG4
82.0132
85.9831
45.4762
84.2253
PG5
29.0385
27.0467
54.9945
7.2104
PG6
52.1248
53.0102
46.0154
65.3609
f
1 fuel cost
612.0513
614.0154
622.5149
631.4003
f
2 (emission)
0.2120
0.2171
0.2225
0.2341
f
3 (loss)
2.1601
2.1598
3.0301
2.8998


## 4.1. Implementation of the CMOABC Algorithm for the OPF Problem

Step 1.
Input the system parameters and the minimum and maximum limits of control variables.Step 2.
Input parameters of the CMOABC algorithm and the lower and upper limits of each variable.Step 3.
Produce the initial population.Step 4.
Cluster the population byK-means and compute the fitness.Step 5.
Store nondominated solutions in the memory table.Step 6.
Optimize every cluster individually.Step 7.
Handle the constraints.Step 8.
Compute the fitness of the solutions and store nondominated solutions into EA.Step 9.
Update each EA of each hive according to greedy selecting strategy, sort the EA based on nondomination, and select the nondomination solutions to stay in EA. If the number of nondominated solutions exceeds the allocated size of EA, apply crowding distance to remove the crowded members.Step 10.
Judge the condition of reclustering. If it meets the condition, the number of clusters can be reset and go to Step4; if not, go to Step 6.


## 4.2. Best Compromise Solution Using Fuzzy Membership Approach

In this work, a fuzzy-based mechanism is employed to extract the best compromise solution over the trade-off curve and assist the decision maker to adjust the generation levels efficiently [48]. Due to imprecise nature of the decision maker’s judgment, each objective function of the ith solution is represented by a membership function μi defined as follows:(32)μi=1Fi≤min⁡⁡Fi,max⁡⁡Fi-Fimax⁡⁡Fi-min⁡⁡Fimin⁡⁡Fi≤Fi≤max⁡⁡Fi0Fi≥max⁡⁡Fi,where min⁡⁡(Fi) and max⁡⁡(Fi) are lower and upper bounds of the ith objective function. The higher the values of the membership function are, the greater the solution satisfaction is.For each nondominated solution, the normalized membership functionμk is calculated as(33)μk=∑i=1Nobjμik∑k=1M∑I=1Nobjμik,where M is the number of nondominated solutions and Nobj is the number of objects. The best compromise solution is the one having the maximum of μk.


## 4.3. Simulation Results

In order to validate the robustness of the proposed CMOABC method, a standard IEEE 30-bus system has been used as the test system [28]. The system represents a portion of the American Electric Power System (in the Midwestern USA) as of December 1961, which consists of six generators, 41 transmission lines, and 4 transformers with off-nominal tap ratio in lines 6–9, 6–10, 4–12, and 27-28. The system data with minimum and maximum limits of control variables are shown in Table 7 [49]. The limits of generator buses and load buses are between 0.95 and 1.1 p.u. and 0.9 and 1.05 p.u., respectively. The lower and upper limits of transformer taps are 0.9 p.u. and 1.05 p.u., respectively, and the step size is 0.01 p.u. The parameter settings for these four algorithms are the same as in “performance setting” in Section 3.4.3. In the following section, there are two cases of the test system considered.(1)
Two objectives: emission cost, emission loss, and loss cost are considered, respectively.(2)
Three objectives: they are emission, cost, and loss.Table 7
Characteristics of the generation units.
G1
G2
G3
G4
G5
G6
Generator limits
P
C
max
⁡ (MW)
150
150
150
150
150
150
P
C
min
⁡ (MW)
5
5
5
5
5
5
Cost coefficients
a
10
10
20
10
20
10
b
200
150
180
100
180
150
c
100
120
40
60
40
100
Emission coefficients
α
4
.
0
9
1
2
.
5
4
3
4
.
2
5
8
5
.
3
2
6
4
.
2
5
8
6
.
1
3
1
β
-
5
.
55
4
-
6
.
40
7
-
5
.
09
4
-
3
.
55
0
-
5
.
09
4
-
5
.
55
5
γ
6
.
4
9
0
5
.
6
3
8
4
.
5
8
6
3
.
3
8
0
4
.
5
8
6
5
.
1
5
1
δ
2
.
0
e
-
4
5
.
0
e
-
4
1
.
0
e
-
6
2
.
0
e
-
3
1
.
0
e
-
6
1
.
0
e
-
5
ε
2
.
8
5
7
3
.
3
3
3
8
.
0
0
0
2
.
0
0
0
8
.
0
0
0
6
.
6
6
7


### 4.3.1. Case I: Two-Objective OPF Optimization

In this case, the results obtained from the proposed CMOABC have been compared with MOPSO, MOABC, and NSGA-II in three cases: minimum cost, minimum emission, and minimum loss. Figures12–14 show the convergence characteristics of the four algorithms on three two-objective OPF problems, respectively. For the case of emission cost, Table 8 gives the best solutions for minimum cost and minimum emission; the data of emission loss and loss cost are shown in Tables 10 and 12. Finally, the best compromise solutions for each objective in the two-dimensional Pareto front are shown in Tables 9, 11, and 13.Table 8
The best solutions for cost and emission from the Pareto front based on four multiobjective algorithms.
CMOABC
MOABC
NSGA-II
MOPSO
Bestf1
Bestf2
Bestf1
Bestf2
Bestf1
Bestf2
Bestf1
Bestf2
PG1
12.36
41.30
11.09
7.92
36.25
42.01
26.11
56.01
PG2
29.99
45.41
30.16
42.08
55.45
47.91
32.99
46.62
PG3
55.40
55.31
57.50
60.89
49.60
50.97
62.14
56.12
PG4
105.16
43.40
105.41
101.06
46.20
43.10
102.13
40.31
PG5
45.26
52.79
45.39
42.61
56.92
52.20
45.21
56.03
PG6
38.01
52.01
37.81
33.88
44.32
52.06
38.15
48.96
f
1
(cost)
606.52
639.09
607.18
645.01
610.80
639.89
607.12
643.90
f
2 (emission)
0.2255
0.1931
0.2259
0.2009
0.2317
0.1992
0.2590
0.1959Table 9
The best compromise solutions for emission and cost using different multiobjective algorithms.
CMOABC
MOABC
NSGA-II
MOPSO
PG1
24.4905
24.9981
36.2012
24.3132
PG2
38.9856
39.5037
55.4215
37.8013
PG3
56.6471
39.0035
50.1034
55.4107
PG4
72.0001
73.9851
45.9804
72.1045
PG5
51.0322
48.1023
56.7635
52.0511
PG6
45.0258
42.9861
44.2014
43.9648
f
1 (cost)
616.7902
617.4896
627.0113
617.5842
f
2 (emission)
0.2015
0.2033
0.1968
0.2019Table 10
The best solutions for cost and loss from the Pareto front based on four multiobjective algorithms.
CMOABC
MOABC
NSGA-II
MOPSO
Bestf1
Bestf3
Bestf1
Bestf3
Bestf1
Bestf3
Bestf1
Bestf3
PG1
15.13
3.69
9.62
4.02
19.06
3.85
66.98
9.54
PG2
29.92
27.41
32.16
11.24
20.65
23.66
37.99
10.02
PG3
55.89
105.93
57.99
109.92
39.62
72.49
68.16
109.34
PG4
112.66
57.01
95.41
67.53
118.69
99.04
102.03
69.36
PG5
31.23
5.23
44.38
11.64
42.02
5.27
38.21
71.81
PG6
38.95
87.56
39.85
80.63
47.14
82.31
28.15
33.65
f
1 (cost)
607.26
626.09
609.99
626.12
611.01
643.26
609.42
669.91
f
3 (loss)
2.4514
1.5196
2.4643
1.5489
3.8167
1.5841
2.5389
1.5826Table 11
The best compromise solutions for loss and cost using different multiobjective algorithms.
CMOABC
MOABC
NSGA-II
MOPSO
PG1
4.4112
8.6014
3.6012
11.6012
PG2
25.0134
29.1642
12.4035
20.7035
PG3
79.4701
74.8653
87.6634
83.0107
PG4
106.2073
105.2366
91.9804
96.1342
PG5
8.0312
17.5038
28.0135
12.0511
PG6
62.0252
50.0031
61.8014
70.9652
f
1 (cost)
622.6002
613.3213
623.2813
630.5862
f
3 (loss)
1.6691
1.7097
1.8703
1.6258Table 12
The best solutions for emission and loss from the Pareto front based on four multiobjective algorithms.
CMOABC
MOABC
NSGA-II
MOPSO
Bestf2
Bestf3
Bestf2
Bestf3
Bestf2
Bestf3
Bestf2
Bestf3
PG1
45.13
9.09
33.90
12.06
44.01
63.30
33.01
10.97
PG2
44.92
5.02
60.03
6.48
41.86
144.51
44.69
77.01
PG3
55.96
91.01
49.89
70.01
57.29
108.58
61.98
23.11
PG4
41.02
70.98
34.65
95.32
40.66
62.03
43.03
63.35
PG5
46.63
5.13
63.61
6.39
42.61
5.12
49.96
82.49
PG6
54.02
103.86
46.96
95.22
34.59
34.59
52.59
58.60
f
2 (emission)
0.1958
0.2489
0.1969
0.2519
0.2501
0.2334
0.1962
0.2599
f
3 (loss)
4.1213
1.6072
4.7821
1.6501
4.2996
2.1239
3.8300
1.6640Table 13
The best compromise solutions for loss and emission using different multiobjective algorithms.
CMOABC
MOABC
NSGA-II
MOPSO
PG1
21.1132
68.9014
30.6011
22.6312
PG2
30.9814
23.0042
32.1035
25.1034
PG3
79.4721
30.8973
77.6934
103.0112
PG4
66.2076
65.9956
65.9894
64.0343
PG5
18.0382
77.5468
34.0243
6.4802
PG6
72.0258
19.0032
48.0014
85.9639
f
2 (emission)
2.0065
2.0143
2.6470
1.5959
f
3 (loss)
0.2035
0.2132
0.2055
0.2439Figure 12
Pareto fronts obtained by CMOABC, MOPSO, MOABC, and NSGA-II for emission and cost.Figure 13
Pareto fronts obtained by CMOABC, MOPSO, MOABC, and NSGA-II for cost and loss.Figure 14
Pareto fronts obtained by CMOABC, MOPSO, MOABC, and NSGA-II for emission and loss.First, two competing objectives, fuel cost and emission, are considered. From the Pareto-optimal solution shown in Figure12, it is clear that the proposed CMOABC method is able to give well distributed solutions compared to the other three algorithms. From Table 8, it is clear that the minimum fuel cost obtained by CMOABC is 606.52 $/h, which is better than the results optimized by the other three algorithms. And when the best fuel cost effect case is considered, CMOABC method can provide smaller loss. The best compromise solution is 616.7902 $/h and 0.2015 ton/h of CMOABC which can be observed in Table 9. It can once again prove that the proposed method has obtained better solutions.The cost loss objective functions are shown in Figure13 and Table 10. From the figure, we can see that the proposed algorithm obtains better results than the other three algorithms, though the obtained Pareto fronts are not evenly distributed. From Table 10, it is easily observed that CMOABC method can provide lower economic fuel cost when the minimum loss case is considered and smaller loss when the best fuel cost effect case is considered.For emission loss objective functions, we can observe that the algorithms achieve similar performance ranking as for fuel cost emission objective functions.


### 4.3.2. Case II: Three-Objective OPF Optimization

In this case, the three objectives are optimized simultaneously by the four algorithms and the corresponding best solutions are given in Table14. Figure 15 also shows the result values of three competing objectives.Table 14
The best solutions for cost, emission, and loss based on four multiobjective algorithms.
CMOABC
MOABC
MOPSO
NSGA-II
Best
Best
Best
Best
Best
Best
Best
Best
Best
Best
Best
Best
f
1
f
2
f
3
f
1
f
2
f
3
f
1
f
2
f
3
f
1
f
2
f
3
PG1
20.02
44.96
6.01
35.66
3.21
18.95
22.01
36.97
9.01
25.86
36.63
18.65
PG2
21.03
45.01
15.26
49.25
16.74
30.01
25.63
45.76
25.75
28.61
49.87
27.83
PG3
70.25
59.79
106.92
62.73
81.51
61.36
64.59
61.02
83.83
58.02
50.76
95.62
PG4
101.05
38.12
68.95
39.03
69.01
102.65
101.81
45.53
64.92
48.72
38.52
66.75
PG5
46.12
52.68
5.00
45.67
8.96
36.57
40.92
49.89
68.66
43.15
51.34
7.97
PG6
35.19
46.11
87.65
33.52
76.72
39.96
32.56
51.64
55.65
42.84
49.42
68.74
f
1 (cost)
612.01
649.86
639.96
651.02
626.84
613.02
616.98
645.09
663.72
618.02
651.64
645.52
f
2 (emission)
0.2498
0.1956
0.2398
0.2069
0.2497
0.2501
0.3805
0.2101
0.3964
0.2426
0.2193
0.2467
f
3 (loss)
1.9968
4.0117
1.6046
3.9964
1.7013
1.9904
2.0972
2.4657
1.5206
2.3994
2.3652
2.0985Figure 15
Pareto fronts obtained by CMOABC, MOPSO, MOABC, and NSGA-II for fuel cost, emission, and loss.
(a)
CMOABC
(b)
MOABC
(c)
MOPSO
(d)
NSGA-IIAs seen in Figure15, compared to the other three algorithms, it can be observed that the Pareto-optimal solutions obtained by the proposed CMOABC are better distributed on the front with good diversity. Among the other three algorithms, the Pareto-optimal solutions obtained by the standard MOABC are also well distributed and the diversification of them is not as good as the ones obtained by the proposed CMOABC.Furthermore, from Table14, it can be also concluded that the CMOABC is able to discover a well distributed and diverse solution set for this three-objective problem. However, the other three algorithms cannot archive the true Pareto front for the three-objective OPF problem.The compromise solutions are shown in Table15. From this table we can see the proposed algorithm obtains the best emission, cost, though the system loss obtained by CMOABC is 2.1601 MW, which is a little more than 2.1598 MW obtained by MOABC.Table 15
The best compromise solutions for cost, emission, and loss using different multiobjective algorithms.
CMOABC
MOABC
NSGA-II
MOPSO
PG1
19.0732
18.9234
35.6214
22.0152
PG2
32.9875
27.9542
53.9065
15.1023
PG3
68.0549
70.8965
47.6936
90.0136
PG4
82.0132
85.9831
45.4762
84.2253
PG5
29.0385
27.0467
54.9945
7.2104
PG6
52.1248
53.0102
46.0154
65.3609
f
1 fuel cost
612.0513
614.0154
622.5149
631.4003
f
2 (emission)
0.2120
0.2171
0.2225
0.2341
f
3 (loss)
2.1601
2.1598
3.0301
2.8998


## 4.3.1. Case I: Two-Objective OPF Optimization

In this case, the results obtained from the proposed CMOABC have been compared with MOPSO, MOABC, and NSGA-II in three cases: minimum cost, minimum emission, and minimum loss. Figures12–14 show the convergence characteristics of the four algorithms on three two-objective OPF problems, respectively. For the case of emission cost, Table 8 gives the best solutions for minimum cost and minimum emission; the data of emission loss and loss cost are shown in Tables 10 and 12. Finally, the best compromise solutions for each objective in the two-dimensional Pareto front are shown in Tables 9, 11, and 13.Table 8
The best solutions for cost and emission from the Pareto front based on four multiobjective algorithms.
CMOABC
MOABC
NSGA-II
MOPSO
Bestf1
Bestf2
Bestf1
Bestf2
Bestf1
Bestf2
Bestf1
Bestf2
PG1
12.36
41.30
11.09
7.92
36.25
42.01
26.11
56.01
PG2
29.99
45.41
30.16
42.08
55.45
47.91
32.99
46.62
PG3
55.40
55.31
57.50
60.89
49.60
50.97
62.14
56.12
PG4
105.16
43.40
105.41
101.06
46.20
43.10
102.13
40.31
PG5
45.26
52.79
45.39
42.61
56.92
52.20
45.21
56.03
PG6
38.01
52.01
37.81
33.88
44.32
52.06
38.15
48.96
f
1
(cost)
606.52
639.09
607.18
645.01
610.80
639.89
607.12
643.90
f
2 (emission)
0.2255
0.1931
0.2259
0.2009
0.2317
0.1992
0.2590
0.1959Table 9
The best compromise solutions for emission and cost using different multiobjective algorithms.
CMOABC
MOABC
NSGA-II
MOPSO
PG1
24.4905
24.9981
36.2012
24.3132
PG2
38.9856
39.5037
55.4215
37.8013
PG3
56.6471
39.0035
50.1034
55.4107
PG4
72.0001
73.9851
45.9804
72.1045
PG5
51.0322
48.1023
56.7635
52.0511
PG6
45.0258
42.9861
44.2014
43.9648
f
1 (cost)
616.7902
617.4896
627.0113
617.5842
f
2 (emission)
0.2015
0.2033
0.1968
0.2019Table 10
The best solutions for cost and loss from the Pareto front based on four multiobjective algorithms.
CMOABC
MOABC
NSGA-II
MOPSO
Bestf1
Bestf3
Bestf1
Bestf3
Bestf1
Bestf3
Bestf1
Bestf3
PG1
15.13
3.69
9.62
4.02
19.06
3.85
66.98
9.54
PG2
29.92
27.41
32.16
11.24
20.65
23.66
37.99
10.02
PG3
55.89
105.93
57.99
109.92
39.62
72.49
68.16
109.34
PG4
112.66
57.01
95.41
67.53
118.69
99.04
102.03
69.36
PG5
31.23
5.23
44.38
11.64
42.02
5.27
38.21
71.81
PG6
38.95
87.56
39.85
80.63
47.14
82.31
28.15
33.65
f
1 (cost)
607.26
626.09
609.99
626.12
611.01
643.26
609.42
669.91
f
3 (loss)
2.4514
1.5196
2.4643
1.5489
3.8167
1.5841
2.5389
1.5826Table 11
The best compromise solutions for loss and cost using different multiobjective algorithms.
CMOABC
MOABC
NSGA-II
MOPSO
PG1
4.4112
8.6014
3.6012
11.6012
PG2
25.0134
29.1642
12.4035
20.7035
PG3
79.4701
74.8653
87.6634
83.0107
PG4
106.2073
105.2366
91.9804
96.1342
PG5
8.0312
17.5038
28.0135
12.0511
PG6
62.0252
50.0031
61.8014
70.9652
f
1 (cost)
622.6002
613.3213
623.2813
630.5862
f
3 (loss)
1.6691
1.7097
1.8703
1.6258Table 12
The best solutions for emission and loss from the Pareto front based on four multiobjective algorithms.
CMOABC
MOABC
NSGA-II
MOPSO
Bestf2
Bestf3
Bestf2
Bestf3
Bestf2
Bestf3
Bestf2
Bestf3
PG1
45.13
9.09
33.90
12.06
44.01
63.30
33.01
10.97
PG2
44.92
5.02
60.03
6.48
41.86
144.51
44.69
77.01
PG3
55.96
91.01
49.89
70.01
57.29
108.58
61.98
23.11
PG4
41.02
70.98
34.65
95.32
40.66
62.03
43.03
63.35
PG5
46.63
5.13
63.61
6.39
42.61
5.12
49.96
82.49
PG6
54.02
103.86
46.96
95.22
34.59
34.59
52.59
58.60
f
2 (emission)
0.1958
0.2489
0.1969
0.2519
0.2501
0.2334
0.1962
0.2599
f
3 (loss)
4.1213
1.6072
4.7821
1.6501
4.2996
2.1239
3.8300
1.6640Table 13
The best compromise solutions for loss and emission using different multiobjective algorithms.
CMOABC
MOABC
NSGA-II
MOPSO
PG1
21.1132
68.9014
30.6011
22.6312
PG2
30.9814
23.0042
32.1035
25.1034
PG3
79.4721
30.8973
77.6934
103.0112
PG4
66.2076
65.9956
65.9894
64.0343
PG5
18.0382
77.5468
34.0243
6.4802
PG6
72.0258
19.0032
48.0014
85.9639
f
2 (emission)
2.0065
2.0143
2.6470
1.5959
f
3 (loss)
0.2035
0.2132
0.2055
0.2439Figure 12
Pareto fronts obtained by CMOABC, MOPSO, MOABC, and NSGA-II for emission and cost.Figure 13
Pareto fronts obtained by CMOABC, MOPSO, MOABC, and NSGA-II for cost and loss.Figure 14
Pareto fronts obtained by CMOABC, MOPSO, MOABC, and NSGA-II for emission and loss.First, two competing objectives, fuel cost and emission, are considered. From the Pareto-optimal solution shown in Figure12, it is clear that the proposed CMOABC method is able to give well distributed solutions compared to the other three algorithms. From Table 8, it is clear that the minimum fuel cost obtained by CMOABC is 606.52 $/h, which is better than the results optimized by the other three algorithms. And when the best fuel cost effect case is considered, CMOABC method can provide smaller loss. The best compromise solution is 616.7902 $/h and 0.2015 ton/h of CMOABC which can be observed in Table 9. It can once again prove that the proposed method has obtained better solutions.The cost loss objective functions are shown in Figure13 and Table 10. From the figure, we can see that the proposed algorithm obtains better results than the other three algorithms, though the obtained Pareto fronts are not evenly distributed. From Table 10, it is easily observed that CMOABC method can provide lower economic fuel cost when the minimum loss case is considered and smaller loss when the best fuel cost effect case is considered.For emission loss objective functions, we can observe that the algorithms achieve similar performance ranking as for fuel cost emission objective functions.


## 4.3.2. Case II: Three-Objective OPF Optimization

In this case, the three objectives are optimized simultaneously by the four algorithms and the corresponding best solutions are given in Table14. Figure 15 also shows the result values of three competing objectives.Table 14
The best solutions for cost, emission, and loss based on four multiobjective algorithms.
CMOABC
MOABC
MOPSO
NSGA-II
Best
Best
Best
Best
Best
Best
Best
Best
Best
Best
Best
Best
f
1
f
2
f
3
f
1
f
2
f
3
f
1
f
2
f
3
f
1
f
2
f
3
PG1
20.02
44.96
6.01
35.66
3.21
18.95
22.01
36.97
9.01
25.86
36.63
18.65
PG2
21.03
45.01
15.26
49.25
16.74
30.01
25.63
45.76
25.75
28.61
49.87
27.83
PG3
70.25
59.79
106.92
62.73
81.51
61.36
64.59
61.02
83.83
58.02
50.76
95.62
PG4
101.05
38.12
68.95
39.03
69.01
102.65
101.81
45.53
64.92
48.72
38.52
66.75
PG5
46.12
52.68
5.00
45.67
8.96
36.57
40.92
49.89
68.66
43.15
51.34
7.97
PG6
35.19
46.11
87.65
33.52
76.72
39.96
32.56
51.64
55.65
42.84
49.42
68.74
f
1 (cost)
612.01
649.86
639.96
651.02
626.84
613.02
616.98
645.09
663.72
618.02
651.64
645.52
f
2 (emission)
0.2498
0.1956
0.2398
0.2069
0.2497
0.2501
0.3805
0.2101
0.3964
0.2426
0.2193
0.2467
f
3 (loss)
1.9968
4.0117
1.6046
3.9964
1.7013
1.9904
2.0972
2.4657
1.5206
2.3994
2.3652
2.0985Figure 15
Pareto fronts obtained by CMOABC, MOPSO, MOABC, and NSGA-II for fuel cost, emission, and loss.
(a)
CMOABC
(b)
MOABC
(c)
MOPSO
(d)
NSGA-IIAs seen in Figure15, compared to the other three algorithms, it can be observed that the Pareto-optimal solutions obtained by the proposed CMOABC are better distributed on the front with good diversity. Among the other three algorithms, the Pareto-optimal solutions obtained by the standard MOABC are also well distributed and the diversification of them is not as good as the ones obtained by the proposed CMOABC.Furthermore, from Table14, it can be also concluded that the CMOABC is able to discover a well distributed and diverse solution set for this three-objective problem. However, the other three algorithms cannot archive the true Pareto front for the three-objective OPF problem.The compromise solutions are shown in Table15. From this table we can see the proposed algorithm obtains the best emission, cost, though the system loss obtained by CMOABC is 2.1601 MW, which is a little more than 2.1598 MW obtained by MOABC.Table 15
The best compromise solutions for cost, emission, and loss using different multiobjective algorithms.
CMOABC
MOABC
NSGA-II
MOPSO
PG1
19.0732
18.9234
35.6214
22.0152
PG2
32.9875
27.9542
53.9065
15.1023
PG3
68.0549
70.8965
47.6936
90.0136
PG4
82.0132
85.9831
45.4762
84.2253
PG5
29.0385
27.0467
54.9945
7.2104
PG6
52.1248
53.0102
46.0154
65.3609
f
1 fuel cost
612.0513
614.0154
622.5149
631.4003
f
2 (emission)
0.2120
0.2171
0.2225
0.2341
f
3 (loss)
2.1601
2.1598
3.0301
2.8998


## 5. Conclusions

An improved multiobjective ABC algorithm based onK-means cluster, called CMOABC, is proposed in this paper. In the proposed algorithm, we first modify the employed bees’ phase of the original ABC to improve its convergence rate. The proposed algorithm adopts K-means clustering method to partition the population into many clusters and the number of the clusters is changing to implement information exchange among the different clusters. With testing against a set of 6 mathematical benchmark functions (including both two- and three-objective cases), CMOABC has a better performance than the other MOPSO, MOABC, and NSGA-II.Then, the proposed algorithm is used to handle multiobjective OPF problem, and 30-bus IEEE test system is adopted to test the proposed algorithm. By comparing the simulation results of CMOABC, MOABC, MOPSO, and NSGA-II, the proposed method is able to give well distributed Pareto-optimal solutions compared to the other three methods for OPF problem with different objectives.



---

*Source: 762853-2015-05-04.xml*