# An MOEA/D-ACO with PBI for Many-Objective Optimization

Authors: Tianbai Ling; Chen Wang
Journal: Mathematical Problems in Engineering
 (2018)

---

## Abstract

Evolutionary algorithms (EAs) are an important instrument for solving the multiobjective optimization problems (MOPs). It has been observed that the combined ant colony (MOEA/D-ACO) based on decomposition is very promising for MOPs. However, as the number of optimization objectives increases, the selection pressure will be released, leading to a significant reduction in the performance of the algorithm. It is a significant problem and challenge in the MOEA/D-ACO to maintain the balance between convergence and diversity in many-objective optimization problems (MaOPs). In the proposed algorithm, an MOEA/D-ACO with the penalty based boundary intersection distance (PBI) method (MOEA/D-ACO-PBI) is intended to solve the MaOPs. PBI decomposes the problems with many single-objective problems, a weighted vector adjustment method based on clustering, and uses different pheromone matrices to solve different single objectives proposed. Then the solutions are constructed and pheromone was updated. Experimental results on both CF1-CF4 and suits of C-DTLZ benchmarks problems demonstrate the superiority of the proposed algorithm in comparison with three state-of-the-art algorithms in terms of both convergence and diversity.


---

## Body

## 1. Introduction

Over the past few decades, EAs have been successfully implemented to various real-world optimization problems with two or three objectives. However, EAs face many challenges when dealing with MOPs, which have more than three objectives and are commonly referred to as MaOPs [1]. The foremost reason is the case that the proportion of nondominated solutions in the population rises sharply with an increase in the number of objectives [2]. The Pareto dominate relation cannot provide sufficient selection pressure to ensure that the population evolve [3]. In recent years, more and more multiobjective evolutionary optimization algorithms (MOEAs) have been improved for solving MaOPs [4].In recent years, as system functions become more and more complex. MaOPs have become a hot topic in the field of evolutionary multiobjective optimization community [5]. It has been found that EAs have succeeded in solving various optimization problems. However, as the number of optimization objectives increases, sufficient selection pressure will be released, resulting in a significant reduction in the performance of the algorithm. It is a challenging major problem and challenge in the EAs to maintain the balance between convergence and diversity in MaOPs.At present, more and more MOEAs have made extraordinary improvement in solving MaOPs. The most popular method is the Pareto-dominance and modified Pareto-dominance that relaxes the form of Pareto-dominance to increase the selection pressure towards the Pareto front, for example, H Jain and Deb based on the NSGAIII framework to solve generic MaOPs [6] and fuzzy Pareto-dominance [7]. The disadvantage of the above methods is introduced in one or more parameters and needs to be adjusted heuristically. The indicator-based algorithm is promising method for MaOPs. It uses a single indicator value to guide the evolution of the population, such as R2 indicator [8] and hypervolume indicator [9]. The last method is based on decomposition. Decomposition strategy and neighbourhood concept are introduced in this method. As one of the popular algorithms, MOPs based on decomposition (MOEA/D) were proposed by Zhang and Li [10]. The MOEA/D uses aggregation function to compare the solution with the uniformly distributed weight vector to retain solution convergence and diversity [4, 11]. MOEA/D gets a real good convergence and diversity and low computational complexity and thus has an effective method.ACO originally proposed dealing with single-objective combination optimization problems [12]. The successful application of ACO in single-objective optimization has led some researchers to develop ACO algorithms for multiobjective optimization problems. Recently, a multiobjective evolutionary algorithm for multiobjective combinatorial optimization problem based on decomposition and ant colony algorithm (MOEA/D-ACO) was proposed. Ke et al. combined ACO and decomposition-based multiobjective evolutionary algorithm (MOEAD/ACO) to solve biobjective TSP [13]. Despite the advance in adapting MOEAs for leading with MaOPs, the MOEA/D-ACO performed much better than the BicriterionAnt in all the test instances. The MOEA/D-ACO algorithm successfully solves the biobjective TSP problem. However, its performance in cases with three or more objectives requires further study. Murilo Zangari Souza et al. proposed a parallel implementation of MOEA/D-ACO on a graphics processing unit (GPU) using NVIDIA CUDA to improve efficiency and achieve high-quality results within a reasonable execution time [14].At present, the MOEA/D-ACO optimization algorithm has been introduced to solve MaOPs in numerous studies. However, there are few studies on the further use of MOEA/D-ACO to solve the MaOPs. In this paper, we introduce the PBI method for MOEA/D-ACO to solve the MaOPs. And many-objective optimization techniques can be incorporated into our method to maintain convergence and diversity of algorithm.From the discussion above, we propose a new algorithm with MOEA/D-ACO-PBI designed to solve the MaOPs. This paper mainly focuses on the algorithm that promotes the performance of convergence and diversity in many-objective optimization.This remainder of this paper is organized as follows. In Section2, the related work and background were discussed. The MOEA/D-ACO-PBI algorithm is described in Section 3. In Section 5, the experimental studies are provided to demonstrate the efficiency of the proposed method as well as some discussions on this paper. Finally, Section 6 makes a conclusion and points out some future research directions.


## 2. Related Background

### 2.1. Basic Definitions

General MaOPs can be define as follows:(1)minimizeFx=f1x,f2x,…,fmxTsubject  tox∈Ω⊆Rnwhere x =x1,…,xnT is the n-dimensional decision variable vector, Ω is the decision space, and F:Ω→Rm consists of m objective functions, fi(x) (i = 1, 2, …, m), m⩾ 4.Rm is called the objective space.Given two decision vectors x,y∈Ω, x is said to Pareto dominate y (x≺y), if(2)∀i∈1,2,…,m:fix≤fiy∃i∈1,2,…,m:fix<fiyA solution x∗∈Ω is Pareto optimal if there is no x≺x∗. f(x∗) is called a Pareto optimal objective vector. The set of all Pareto optimal objective vectors is called the Pareto front (PF) and the set of all Pareto optimal solutions is called the Pareto optimal set (PS).


### 2.2. MOEA/D

MOEA/D: it is the most popular multiobjective evolutionary algorithm in recent years. It first initializes N uniform and widely distributed vectors for each individual evolving in the direction of each vector, which is equivalent to dividing multiobjective optimization problems. The solution to the N scalar optimization subproblems is through the constant. Evolution can simultaneously obtain the optimal solution of these N subproblems. The algorithm proposes 3 weighted methods that can be used: weighted sum aggregation, Tchebycheff approach, and PBI approach.The traditional weighted sum aggregation is the most direct method to deal with multiobjective optimization problems, and it is also a method used earlier. This method combines or aggregates all the optimized subobjectives into a single-objective, and the multiobjective optimization problem is converted into a single-objective optimization problem.The expression of weighted sum aggregation is(3)mingwsx∣λ=∑i=1mλifixThe weighted sum aggregation method generates a set of different Pareto optimal solutions by generating different weight vectors in the above scalar optimization problem. However, in the case that the shape of the optimal Pareto surface is nonconvex, this method cannot obtain every Pareto optimal vector.The expression of Tchebycheff approach is(4)mingtex∣λ,z∗=max1≤i≤mλi∣fix-zi∗For each Pareto optimality x there is always a weight vector y such that the solution of (4) is a Pareto optimal solution, which corresponds to the Pareto optimal solution of the original multiobjective problem. And by modifying the weight vector, different solutions on the Pareto optimal surface can be obtained. One disadvantage of this approach is that the aggregate method curves are not smooth when dealing with continuous multiobjective problems.The PBI method aims to find the intersection of the uppermost boundary and a set of lines. In a sense, if the lines are evenly distributed, the resulting intersections can be seen as providing a good approximation to the overall Pareto optimal boundary. These methods can deal with the Pareto optimal boundary nonconvex problem. The PBI method will be described in detail later.


### 2.3. MOEA/D-ACO

The basic framework of MOEA/D-ACO first decomposes the MOPs intoNs single-objective subproblems by N weight vectors λ1,…,λN. Then, MOEA/D-ACO utilizes N ants to solve these single-objective subproblems. Ant ia is responsible for subproblems is. Finally, MOEA/D-ACO uses the concept of neighbourhood and group for decomposing the MaOPs. Figure 1 illustrates the method of above.Figure 1
MOEA/D-ACO.Firstly, MOEA/D-ACO decomposes the MOP intoNs single-objective subproblems by choosing Nw weight vectors (λ1…λNw). Then the subproblem is was associated with weight vector λis, and its objective function is denoted as g(x∣λis). MOEA/D-ACO utilizes Na ants for solving these single-objective subproblems, while ant ia is responsible for subproblem is. Both the weighted sum and Tchebycheff approaches are used in MOEA/D-ACO and are recognized as MOEA/D-ACOW and MOEA/D-ACOT, respectively. However, recently, Deb et al. proposed that PBI decomposition method is more suitable for solving MOPs.The MOEA/D-ACO is shown in Algorithm1. In Algorithm 1, for each generation, each ant constructs a solution in Step (3). In Step (4), the newly constructed solution is updated in the external archive, and an ant updates the pheromone matrix of its group in Step (5), if it has found a new nondominated solution. In Step (6), each ant updates its current solution if there is a solution that (1) is better than its current solutions; and (2) has not been used for updating any other old solution. The algorithm is terminated in Step (7).Algorithm 1:MOEA/D-ACO.
( 1 ) Initialization: generate Nini Solutions by weight vectors λ ( 2 ) Decomposition by some aggregation approach ( 3 ) Solution construction ( 4 ) Update of EP ( 5 ) Update of pheromone ( 6 ) Check the solutions in the neighbourhood and update the solutions ( 7 ) Termination: judgment and output EP (a set of high-quality solutions: XS)


### 2.4. PBI Method

In our paper, we use the PBI approach, because of its promising performance for many-objective optimization problems. The optimization problem of the PBI approach is defined as(5)minimizegpbix∣w,z∗=d1+θd2subjecttox∈Ωwhere(6)d1=Fx-z∗TwWd2=Fx-z∗+d1WWZ ∗= z1…z∗mT is the ideal objective vector with zi∗<minfi(x), x∈Ω, i∈{1…m}, a penalty parameter θ predefined by the user. Figure 2 presents an example to illustrate d1 and d2 of a solution xi. In the PBI approach, d1 is used to measure the convergence of xi towards the EF; d2 is a kind of measure for the population diversity. gpbi(x∣w,z∗) is a composite measure of x for both diversity and convergence. The distance between d1 and d2 is controlled by the parameter θ. In our paper, we set θ = 5.0 for empirical studies.Figure 2
Illustration of PBI approach.


### 2.5. Cluster-Based Weight Vector Adjustment

For the MaOPs problem, this paper proposes a weighted vector adjustment method based on clustering and uses different pheromone matrices to solve different single objectives. This scalable approach can help the ant colony algorithm to better solve the MaOPs problem.After the multiobjective problem is decomposed into N single-target subproblems by the PBI method, N ants are selected to solve these single-target subproblems. Ant i solves the subproblem i. The subproblem i is associated withλi. The corresponding objective function is g(x∣λi). N ants are divided into K groups according to their corresponding weight distances. Each group approaches a small area of the Pareto frontier. The same group of ants share a pheromone matrix that accommodates the location information of the learned Pareto frontier subregions. Each ant looks for a special point on the Pareto frontier for a single-target subproblem. Each ant has its own heuristic information matrix to store a priori knowledge of its subproblems. The ant i's neighbour B (i) contains T ants whose weights are closest to λi in N weights, and i∈ B (i), and ant i is its own neighbour.In each iteration of the algorithm, if the solution set in the dominating set remains unchanged for the first time, it means that the current weight vector can no longer meet the need of guiding population evolution, and weight adjustment is needed. In this case, nondominated populations are classified based on the clustering method, assuming that there are 7 individuals in the objective space (x1, x2… x7). They had three clusters, x3, x5, x2, x6, and x1, x7, x4. Among them, x3, x6, and x7 are cluster centres, as shown in Figure 3.Figure 3
Cluster-based weight vector adjustment.


## 2.1. Basic Definitions

General MaOPs can be define as follows:(1)minimizeFx=f1x,f2x,…,fmxTsubject  tox∈Ω⊆Rnwhere x =x1,…,xnT is the n-dimensional decision variable vector, Ω is the decision space, and F:Ω→Rm consists of m objective functions, fi(x) (i = 1, 2, …, m), m⩾ 4.Rm is called the objective space.Given two decision vectors x,y∈Ω, x is said to Pareto dominate y (x≺y), if(2)∀i∈1,2,…,m:fix≤fiy∃i∈1,2,…,m:fix<fiyA solution x∗∈Ω is Pareto optimal if there is no x≺x∗. f(x∗) is called a Pareto optimal objective vector. The set of all Pareto optimal objective vectors is called the Pareto front (PF) and the set of all Pareto optimal solutions is called the Pareto optimal set (PS).


## 2.2. MOEA/D

MOEA/D: it is the most popular multiobjective evolutionary algorithm in recent years. It first initializes N uniform and widely distributed vectors for each individual evolving in the direction of each vector, which is equivalent to dividing multiobjective optimization problems. The solution to the N scalar optimization subproblems is through the constant. Evolution can simultaneously obtain the optimal solution of these N subproblems. The algorithm proposes 3 weighted methods that can be used: weighted sum aggregation, Tchebycheff approach, and PBI approach.The traditional weighted sum aggregation is the most direct method to deal with multiobjective optimization problems, and it is also a method used earlier. This method combines or aggregates all the optimized subobjectives into a single-objective, and the multiobjective optimization problem is converted into a single-objective optimization problem.The expression of weighted sum aggregation is(3)mingwsx∣λ=∑i=1mλifixThe weighted sum aggregation method generates a set of different Pareto optimal solutions by generating different weight vectors in the above scalar optimization problem. However, in the case that the shape of the optimal Pareto surface is nonconvex, this method cannot obtain every Pareto optimal vector.The expression of Tchebycheff approach is(4)mingtex∣λ,z∗=max1≤i≤mλi∣fix-zi∗For each Pareto optimality x there is always a weight vector y such that the solution of (4) is a Pareto optimal solution, which corresponds to the Pareto optimal solution of the original multiobjective problem. And by modifying the weight vector, different solutions on the Pareto optimal surface can be obtained. One disadvantage of this approach is that the aggregate method curves are not smooth when dealing with continuous multiobjective problems.The PBI method aims to find the intersection of the uppermost boundary and a set of lines. In a sense, if the lines are evenly distributed, the resulting intersections can be seen as providing a good approximation to the overall Pareto optimal boundary. These methods can deal with the Pareto optimal boundary nonconvex problem. The PBI method will be described in detail later.


## 2.3. MOEA/D-ACO

The basic framework of MOEA/D-ACO first decomposes the MOPs intoNs single-objective subproblems by N weight vectors λ1,…,λN. Then, MOEA/D-ACO utilizes N ants to solve these single-objective subproblems. Ant ia is responsible for subproblems is. Finally, MOEA/D-ACO uses the concept of neighbourhood and group for decomposing the MaOPs. Figure 1 illustrates the method of above.Figure 1
MOEA/D-ACO.Firstly, MOEA/D-ACO decomposes the MOP intoNs single-objective subproblems by choosing Nw weight vectors (λ1…λNw). Then the subproblem is was associated with weight vector λis, and its objective function is denoted as g(x∣λis). MOEA/D-ACO utilizes Na ants for solving these single-objective subproblems, while ant ia is responsible for subproblem is. Both the weighted sum and Tchebycheff approaches are used in MOEA/D-ACO and are recognized as MOEA/D-ACOW and MOEA/D-ACOT, respectively. However, recently, Deb et al. proposed that PBI decomposition method is more suitable for solving MOPs.The MOEA/D-ACO is shown in Algorithm1. In Algorithm 1, for each generation, each ant constructs a solution in Step (3). In Step (4), the newly constructed solution is updated in the external archive, and an ant updates the pheromone matrix of its group in Step (5), if it has found a new nondominated solution. In Step (6), each ant updates its current solution if there is a solution that (1) is better than its current solutions; and (2) has not been used for updating any other old solution. The algorithm is terminated in Step (7).Algorithm 1:MOEA/D-ACO.
( 1 ) Initialization: generate Nini Solutions by weight vectors λ ( 2 ) Decomposition by some aggregation approach ( 3 ) Solution construction ( 4 ) Update of EP ( 5 ) Update of pheromone ( 6 ) Check the solutions in the neighbourhood and update the solutions ( 7 ) Termination: judgment and output EP (a set of high-quality solutions: XS)


## 2.4. PBI Method

In our paper, we use the PBI approach, because of its promising performance for many-objective optimization problems. The optimization problem of the PBI approach is defined as(5)minimizegpbix∣w,z∗=d1+θd2subjecttox∈Ωwhere(6)d1=Fx-z∗TwWd2=Fx-z∗+d1WWZ ∗= z1…z∗mT is the ideal objective vector with zi∗<minfi(x), x∈Ω, i∈{1…m}, a penalty parameter θ predefined by the user. Figure 2 presents an example to illustrate d1 and d2 of a solution xi. In the PBI approach, d1 is used to measure the convergence of xi towards the EF; d2 is a kind of measure for the population diversity. gpbi(x∣w,z∗) is a composite measure of x for both diversity and convergence. The distance between d1 and d2 is controlled by the parameter θ. In our paper, we set θ = 5.0 for empirical studies.Figure 2
Illustration of PBI approach.


## 2.5. Cluster-Based Weight Vector Adjustment

For the MaOPs problem, this paper proposes a weighted vector adjustment method based on clustering and uses different pheromone matrices to solve different single objectives. This scalable approach can help the ant colony algorithm to better solve the MaOPs problem.After the multiobjective problem is decomposed into N single-target subproblems by the PBI method, N ants are selected to solve these single-target subproblems. Ant i solves the subproblem i. The subproblem i is associated withλi. The corresponding objective function is g(x∣λi). N ants are divided into K groups according to their corresponding weight distances. Each group approaches a small area of the Pareto frontier. The same group of ants share a pheromone matrix that accommodates the location information of the learned Pareto frontier subregions. Each ant looks for a special point on the Pareto frontier for a single-target subproblem. Each ant has its own heuristic information matrix to store a priori knowledge of its subproblems. The ant i's neighbour B (i) contains T ants whose weights are closest to λi in N weights, and i∈ B (i), and ant i is its own neighbour.In each iteration of the algorithm, if the solution set in the dominating set remains unchanged for the first time, it means that the current weight vector can no longer meet the need of guiding population evolution, and weight adjustment is needed. In this case, nondominated populations are classified based on the clustering method, assuming that there are 7 individuals in the objective space (x1, x2… x7). They had three clusters, x3, x5, x2, x6, and x1, x7, x4. Among them, x3, x6, and x7 are cluster centres, as shown in Figure 3.Figure 3
Cluster-based weight vector adjustment.


## 3. Algorithm Description

### 3.1. MOEA/D-ACO-PBI

In this paper, we propose an algorithm called MOEA/D-ACO-PBI. The original MOEA/D-ACO used a weighted sum and Tchebycheff method to solve the problem. However, this is difficult; because the correlation between the algebraic and geometrical interpretation of the weight sets is not valid [15]. Recent studies conducted by [16] have shown that MOEA/D-PBI is the better method for many-objective optimization problems. Therefore, we introduce the PBI approach to MOEA/D-ACO, known as MOEA/D-ACO-PBI. The MOEA/D-ACO-PBI algorithm is presented in Algorithm 2.Algorithm 2:MOEA/D-ACO-PBI.
( 1 ) Initialization ( 2 ) Decomposition by PBI approach ( 3 ) Solution construction ( 4 ) Update of EP ( 5 ) Update of pheromone ( 6 ) Update the solutions ( 7 ) Termination


### 3.2. Initialization Procedure

Each antia has a heuristic information value ηk,lia, each group j has pheromone trail τia,jk, and each individual solution xia is a tour represented by permutation.The initialization procedure of MOEA/D-ACO-PBI contains six main steps.( 1 ) Setting of Nw and λNw: the sets of weight vectors λNw = {λ1…λNw} are generated by a systematic approach, developed from Das and Dennis's method [17]. Here, the size of λ is Nw.(7)Nw=CH+m-1m-1where H > 0 is the quantity of divisions considered along each objective coordinate. Therefore, a subproblemis has a weight vector λis that satisfies ∑l=1mλlis=1 and λkis≥0, where m is the number of objectives.( 2 ) Initialize the solution xi to the subproblem is and set Fi= F (xi), for i = 1…n.( 3 ) Set group: the number of groups is J, which is a value parameter that requires to be selected. The number of subproblems of each group is equally distributed by clustering the weights in terms of the Euclidean distance.( 4 ) Initialize the pheromone information matrix τi,jk=1, for all j = 1…k and k = 1…n. The motivation is to encourage the search to focus on exploration in its early stage.( 5 ) Initialize ηk,li=1/∑j=1mλjick,lj, which is the heuristic information matrix to the subproblem is, for j = 1…k.( 6 ) Initialize EP as the set of the entire high-quality solutions in {F1…FN}; EP is initialized to be empty.


### 3.3. Solution Construction

Assume that antia is in group j and its current solution xi = (x1i…xni). Ant ia constructs its new solution through the following steps.( 1 ) The probability of choice of the link is set. For k, l = 1…n, set(8)Φk,l=Γk,lia+Δ∗lnxia,k,laηk,liaβΦ k , l is representing the attractiveness of the link between the k and l to ant ia. The indicator function in (xia, (k. l)) is equal to 1 if the link (k, l) is in tour xia or 0 if otherwise. α, β, and Δ are the three control parameters.( 2 ) First, ant ia chooses an ant randomly as its starting point and builds its tour. Suppose that its current position is l and it has not completed its tour. From S (not visited so far), it selects the ant k to visit next, according to the following probability by the roulette wheel selection:(9)Φk,l∑s∈CΦs,l( 3 ) If the ant has visited all the points, it returns its tour.


### 3.4. Update of Pheromone

The pheromone trail value of the link (k, l) for group j is updated as follows:(10)τk,lia≔ρτk,lj+∑x∈Π1gx∣λiaρ is the persistence rate of the pheromone trail. Π is the set of all the new solutions; it was constructed by the ants in group j in the current iteration and was just added to EP. It contains the link between the ants’ k and l. τmin and τmax are used to limit the range of the pheromone. τmin is given by(11)τmin=ετmaxwhere0<ε<1 is a control parameter. Then, if τk,lia< τmin set τk,lia = τmin and if τk,lia > τmax set τk,lia = τmax.τmax is set(12)τmax=B+11-ρgminwhere B is the number of nondominated solutions found in the current iteration andgmin is the smallest value obtained for the objective functions of all the Ns subproblems.


## 3.1. MOEA/D-ACO-PBI

In this paper, we propose an algorithm called MOEA/D-ACO-PBI. The original MOEA/D-ACO used a weighted sum and Tchebycheff method to solve the problem. However, this is difficult; because the correlation between the algebraic and geometrical interpretation of the weight sets is not valid [15]. Recent studies conducted by [16] have shown that MOEA/D-PBI is the better method for many-objective optimization problems. Therefore, we introduce the PBI approach to MOEA/D-ACO, known as MOEA/D-ACO-PBI. The MOEA/D-ACO-PBI algorithm is presented in Algorithm 2.Algorithm 2:MOEA/D-ACO-PBI.
( 1 ) Initialization ( 2 ) Decomposition by PBI approach ( 3 ) Solution construction ( 4 ) Update of EP ( 5 ) Update of pheromone ( 6 ) Update the solutions ( 7 ) Termination


## 3.2. Initialization Procedure

Each antia has a heuristic information value ηk,lia, each group j has pheromone trail τia,jk, and each individual solution xia is a tour represented by permutation.The initialization procedure of MOEA/D-ACO-PBI contains six main steps.( 1 ) Setting of Nw and λNw: the sets of weight vectors λNw = {λ1…λNw} are generated by a systematic approach, developed from Das and Dennis's method [17]. Here, the size of λ is Nw.(7)Nw=CH+m-1m-1where H > 0 is the quantity of divisions considered along each objective coordinate. Therefore, a subproblemis has a weight vector λis that satisfies ∑l=1mλlis=1 and λkis≥0, where m is the number of objectives.( 2 ) Initialize the solution xi to the subproblem is and set Fi= F (xi), for i = 1…n.( 3 ) Set group: the number of groups is J, which is a value parameter that requires to be selected. The number of subproblems of each group is equally distributed by clustering the weights in terms of the Euclidean distance.( 4 ) Initialize the pheromone information matrix τi,jk=1, for all j = 1…k and k = 1…n. The motivation is to encourage the search to focus on exploration in its early stage.( 5 ) Initialize ηk,li=1/∑j=1mλjick,lj, which is the heuristic information matrix to the subproblem is, for j = 1…k.( 6 ) Initialize EP as the set of the entire high-quality solutions in {F1…FN}; EP is initialized to be empty.


## 3.3. Solution Construction

Assume that antia is in group j and its current solution xi = (x1i…xni). Ant ia constructs its new solution through the following steps.( 1 ) The probability of choice of the link is set. For k, l = 1…n, set(8)Φk,l=Γk,lia+Δ∗lnxia,k,laηk,liaβΦ k , l is representing the attractiveness of the link between the k and l to ant ia. The indicator function in (xia, (k. l)) is equal to 1 if the link (k, l) is in tour xia or 0 if otherwise. α, β, and Δ are the three control parameters.( 2 ) First, ant ia chooses an ant randomly as its starting point and builds its tour. Suppose that its current position is l and it has not completed its tour. From S (not visited so far), it selects the ant k to visit next, according to the following probability by the roulette wheel selection:(9)Φk,l∑s∈CΦs,l( 3 ) If the ant has visited all the points, it returns its tour.


## 3.4. Update of Pheromone

The pheromone trail value of the link (k, l) for group j is updated as follows:(10)τk,lia≔ρτk,lj+∑x∈Π1gx∣λiaρ is the persistence rate of the pheromone trail. Π is the set of all the new solutions; it was constructed by the ants in group j in the current iteration and was just added to EP. It contains the link between the ants’ k and l. τmin and τmax are used to limit the range of the pheromone. τmin is given by(11)τmin=ετmaxwhere0<ε<1 is a control parameter. Then, if τk,lia< τmin set τk,lia = τmin and if τk,lia > τmax set τk,lia = τmax.τmax is set(12)τmax=B+11-ρgminwhere B is the number of nondominated solutions found in the current iteration andgmin is the smallest value obtained for the objective functions of all the Ns subproblems.


## 4. Simulation Results

In this section, we compare the MOEA/D-ACO-PBI algorithm with three state-of-the-art algorithms including MOEA/D-ACOW[13], MOEA/D [10], and Multiple Ant Colony System (MACS) [18]. The test problems are introduced in Section 4.1. Section 4.2 describes the quality indicators. Three state-of-the-art algorithms used for comparison and the corresponding parameter settings are briefly introduced in Section 4.3.


### 4.1. Test Problems

The well-known MaOPs test function definitions of CF1–CF4 [19] are used to test the proposed algorithm. Table 1 shows the CF1–CF4 problems with the CF test suites. CF1-CF4 are test functions with biobjective problem. CF test suites are used to test the low dimensional objective function. To verify the ability of algorithms to deal with more objectives, we select C1-DTLZ1, C2-DTLZ2, and C3-DTLZ4 [20] used to test the high dimensional objective function. CF suites are popular test problem for MaOPs.Table 1
Test functions utilized in this experiment.
function Formulation CF1 f1(x)=x1+2J1∑j∈J1xj-x10.51.0+3j-2/(n-2)2 f 2 ( x ) =  1 - x 1 + 2 J 2 ∑ j ∈ J 2 x j - x 1 0.5 1.0 + 3 j - 2 / ( n - 2 ) 2 CF2 f1(x)=x1+2J1∑j∈J1xj-sin6πx1+jπn2 f 2 ( x ) =  1 - x 1 + 2 J 2 ∑ j ∈ J 2 x j - c o s 6 π x 1 + j π n 2 CF3 f1(x)=x1+2J14∑j∈J1yj2-2∏j∈J1cos20yjπj+2 CF4 f 2 ( x ) =  1 - x 1 2 + 2 J 2 4 ∑ j ∈ J 2 y j 2 - 2 ∏ j ∈ J 2 c o s 20 y j π j + 2 f 1 =  x 1 + ∑ j ∈ J 1 h j ( y j ) f 1 =  1 - x 1 + ∑ j ∈ J 2 h j ( y j )C1-DTLZ1 is a MaOPs problem; the problem is formulated as(13)Cx=1-fmx0.6-∑i=1m-1fix0.5≥0C2-DTLZ2 is a MaOPs problem; the problem is formulated as(14)Cx=-minmini=1m⁡fix-12+∑j=1,j≠imfxj2-r2,×∑i=1mfix-1m2-r2≥0C3-DTLZ4 is a MaOPs problem; the problem is formulated as(15)Cix=fi2x4+∑j=1,j≠imfi2x-1≥0


### 4.2. Quality Indicators

In our empirical study, the following two widely used quality indicators were considered. The first one reflects the convergence of an algorithm, while the second one reflects the convergence and diversity of the solutions simultaneously.


#### 4.2.1. Generational Distance (GD) Indicator

For any kind of algorithm, let P be the set of final nondominated points obtained from the objective space. LetP∗ be a set of points uniformly spread over the true PF. The GD can indicate only the convergence of an algorithm [21], and a smaller value indicates better quality. The GD is computed as(16)GDP,P∗=∑u∈Pdu,P∗P


#### 4.2.2. Inverted Generational Distance (IGD) Indicator

For any kind of algorithm, let P be the set of final nondominated points obtained from the objective space andP∗ be a set of points uniformly spread over the true PF. The IGD can indicate both the convergence and diversity [21], and a smaller value indicates better quality. The IGD is computed as(17)IGDP,P∗=∑v∈P∗dv,PP∗


### 4.3. Parameter Settings

The parameters for the four MOEAs considered in this paper are listed as follows.( 1 ) Some parameter settings in MOEA/D-ACO-PBI: set T=8, α=1, β=8, ρ=0.95, r=0.8, and θ=5.0.( 2 ) Some parameter settings in MOEA/D-ACOW and MOEA/D: T=8, α=1, β=8, ρ=0.95, and r=0.8. The other parameters can be obtained from [15].( 3 ) Some parameter settings in MACS: α=1, β=2, ρ=0.9, maxFcctor=5, and AntNum=10.( 4 ) Number of runs and termination condition: all algorithms were programmed in PlatEMO [22], and each algorithm was independently run 30 times for each test instance and stopped after 300 generations. All the algorithms were executed on a desktop PC with a 2.1-GHz CPU, 8-GB RAM, and Windows 10.


## 4.1. Test Problems

The well-known MaOPs test function definitions of CF1–CF4 [19] are used to test the proposed algorithm. Table 1 shows the CF1–CF4 problems with the CF test suites. CF1-CF4 are test functions with biobjective problem. CF test suites are used to test the low dimensional objective function. To verify the ability of algorithms to deal with more objectives, we select C1-DTLZ1, C2-DTLZ2, and C3-DTLZ4 [20] used to test the high dimensional objective function. CF suites are popular test problem for MaOPs.Table 1
Test functions utilized in this experiment.
function Formulation CF1 f1(x)=x1+2J1∑j∈J1xj-x10.51.0+3j-2/(n-2)2 f 2 ( x ) =  1 - x 1 + 2 J 2 ∑ j ∈ J 2 x j - x 1 0.5 1.0 + 3 j - 2 / ( n - 2 ) 2 CF2 f1(x)=x1+2J1∑j∈J1xj-sin6πx1+jπn2 f 2 ( x ) =  1 - x 1 + 2 J 2 ∑ j ∈ J 2 x j - c o s 6 π x 1 + j π n 2 CF3 f1(x)=x1+2J14∑j∈J1yj2-2∏j∈J1cos20yjπj+2 CF4 f 2 ( x ) =  1 - x 1 2 + 2 J 2 4 ∑ j ∈ J 2 y j 2 - 2 ∏ j ∈ J 2 c o s 20 y j π j + 2 f 1 =  x 1 + ∑ j ∈ J 1 h j ( y j ) f 1 =  1 - x 1 + ∑ j ∈ J 2 h j ( y j )C1-DTLZ1 is a MaOPs problem; the problem is formulated as(13)Cx=1-fmx0.6-∑i=1m-1fix0.5≥0C2-DTLZ2 is a MaOPs problem; the problem is formulated as(14)Cx=-minmini=1m⁡fix-12+∑j=1,j≠imfxj2-r2,×∑i=1mfix-1m2-r2≥0C3-DTLZ4 is a MaOPs problem; the problem is formulated as(15)Cix=fi2x4+∑j=1,j≠imfi2x-1≥0


## 4.2. Quality Indicators

In our empirical study, the following two widely used quality indicators were considered. The first one reflects the convergence of an algorithm, while the second one reflects the convergence and diversity of the solutions simultaneously.


### 4.2.1. Generational Distance (GD) Indicator

For any kind of algorithm, let P be the set of final nondominated points obtained from the objective space. LetP∗ be a set of points uniformly spread over the true PF. The GD can indicate only the convergence of an algorithm [21], and a smaller value indicates better quality. The GD is computed as(16)GDP,P∗=∑u∈Pdu,P∗P


### 4.2.2. Inverted Generational Distance (IGD) Indicator

For any kind of algorithm, let P be the set of final nondominated points obtained from the objective space andP∗ be a set of points uniformly spread over the true PF. The IGD can indicate both the convergence and diversity [21], and a smaller value indicates better quality. The IGD is computed as(17)IGDP,P∗=∑v∈P∗dv,PP∗


## 4.2.1. Generational Distance (GD) Indicator

For any kind of algorithm, let P be the set of final nondominated points obtained from the objective space. LetP∗ be a set of points uniformly spread over the true PF. The GD can indicate only the convergence of an algorithm [21], and a smaller value indicates better quality. The GD is computed as(16)GDP,P∗=∑u∈Pdu,P∗P


## 4.2.2. Inverted Generational Distance (IGD) Indicator

For any kind of algorithm, let P be the set of final nondominated points obtained from the objective space andP∗ be a set of points uniformly spread over the true PF. The IGD can indicate both the convergence and diversity [21], and a smaller value indicates better quality. The IGD is computed as(17)IGDP,P∗=∑v∈P∗dv,PP∗


## 4.3. Parameter Settings

The parameters for the four MOEAs considered in this paper are listed as follows.( 1 ) Some parameter settings in MOEA/D-ACO-PBI: set T=8, α=1, β=8, ρ=0.95, r=0.8, and θ=5.0.( 2 ) Some parameter settings in MOEA/D-ACOW and MOEA/D: T=8, α=1, β=8, ρ=0.95, and r=0.8. The other parameters can be obtained from [15].( 3 ) Some parameter settings in MACS: α=1, β=2, ρ=0.9, maxFcctor=5, and AntNum=10.( 4 ) Number of runs and termination condition: all algorithms were programmed in PlatEMO [22], and each algorithm was independently run 30 times for each test instance and stopped after 300 generations. All the algorithms were executed on a desktop PC with a 2.1-GHz CPU, 8-GB RAM, and Windows 10.


## 5. Results and Discussion

The GD indicator was utilized to compare the convergence capacity between the proposed MOEA/D-ACO-PBI algorithm and three other algorithms (Figure4). Table 2 shows the average and standard deviation of the GD values obtained by the four algorithms for test suits with different number of objectives. Table 2 shows that the MOEA/D-ACO-PBI algorithm could perform well on all the test functions, especially on CF1 and CF3 problems. When the number of objective test functions is higher, the performance of MOEA/D-ACO-PBI algorithm is better. The results of comparison of the MOEA/D-ACO-PBI algorithm with the other three MOEAs in terms of IGD values, using the test suite, are provided in Table 3. It shows the average IGD values for the four MOEAs over 30 independent runs and the best average and standard deviation values. From the experimental results of the test problem, the MOEA/D-ACO-PBI algorithm performs on CF2, CF4, C1-DTLZ1, C2-DTLZ2, and C3-DTLZ4 test problems. It is shown that the MOEA/D-ACO-PBI algorithm has better performance of convergence and diversity. In addition, performance of MOEA/D-ACOW algorithm is also better, and MACS performance is the worst.Table 2
The average and standard deviation of GD values obtained for CF1–CF4, C1_DTLZ1, C2_DTLZ2, and C3_DTLZ4 with four algorithms and different number of objectives. Best performance is in bold face.
M MOEA/D-ACO-PBI MOEA/D-ACOW MOEA/D MACS CF1 2 1.0013e-4 (1.14e-4) 1.1945e-4 (2.08e-4) 1.2524e-4 (2.12e-2) 3.4214e-2 (2.92e-2) CF2 2 2.3325e-2 (1.05e-2) 2.0382e-2 (2.59e-2) 2.9194e-2 (1.24e-2) 4.5986e-2 (4.28e-2) CF3 2 2.2667e-1 (1.68e-1) 3.7510e-1 (2.41e-1) 3.3477e-1 (3.19e-1) 4.2804e-1 (1.32e-1) CF4 2 4.4754e-2 (1.29e-1) 5.4257e-2 (1.38e-1) 4.2232e-2 (7.19e-2) 2.1410e-1 (2.78e-1) 3 3.0152e-3 (2.13e-4) 3.2859e-3 (2.36e-4) 3.4831e-3 (3.32e-4) 4.4339e-3 (1.33e-4) C1_DTLZ1 4 3.0221-3 (2.23e-4) 2.5532e-2 (2.23e-3) 2.1234e-2 (1.34e-3) 2.4452e-2 (2.12e-2) 5 4.4264e-3 (3.25e-4) 3.3354e-2 (2.23e-2) 5.4261e-2 (2.34e-3) 5.3634e-2 (2.84e-2) C2_DTLZ2 3 7.4664e-4 (5.31e-5) 2.2772e-3 (4.29e-4) 9.5259e-3 (1.67e-3) 7.9976e-3 (5.29e-3) 4 2.7056e-3 (1.15e-4) 1.1709e-2 (4.81e-3) 1.1867e-2 (1.28e-3) 2.7389e-2 (1.08e-2) 5 5.2944e-3 (2.99e-4) 2.5550e-2 (1.10e-2) 2.7631e-2 (6.95e-3) 5.9090e-2(2.84e-4) C3_DTLZ4 3 4.0236e-4 (2.23e-4) 6.9914e-4 (3.28e-4) 4.5365e-4 (2.65e-4) 3.514e-3 (3.26e-3) 4 3.0885e-3 (2.32e-3) 2.2035e-2 (3.82e-2) 1.3564e-2 (1.36e-2) 2.5658e-2 (1.25e-4) 5 3.0624e-3 (3.25e-4) 3.6412e-2 (2.21e-2) 3.7225e-2 (5.62e-3) 5.3361e-2 (3.26e-4)Table 3
The average and standard deviation of IGD values obtained for CF1–CF4, C1_DTLZ1, C2_DTLZ2, and C3_DTLZ4 with four algorithms and different number of objectives. Best performance is in bold face.
M MOEA/D-ACO-PBI MOEA/D-ACOW MOEA/D MACS CF1 2 6.4371e-1 (1.36e-1) 5.1532e-1 (1.27e-1) 6.5734e-1 (1.08e-1) 7.8365e-1 (1.89e-1) CF2 2 4.5894e-2 (1.28e-2) 4.8905e-2 (2.52e-2) 6.2428e2 (6.10e-2) 2.4116e-1 (3.61e-1) CF3 2 3.7740e-1 (7.78e-2) 3.1772e-1 (9.32e-2) 3.7877e-1 (1.77e-1) 3.7581e-1 (2.54e-1) CF4 2 1.6256e-1 (6.34e-2) 2.4364e-1 (1.18e-1) 1.9728e-1 (7.17e-2) 2.4654e-1 (1.32e-1) C1_DTLZ1 4 2.2653e-1 (3.25e-4) 3.4629e-1 (2.36e-2) 2.5698e-1 3.26e-3) 3.7241e-1 (3.21e-4) 5 3.1651e-1 (4.13e-4) 3.6234e-1 (1.12e-1) 3.7511e-1 (1.69e-1) 5.9283e-1 (3.87e-4) C2_DTLZ2 3 5.1902e-2 (5.09e-4) 6.5049e-2 (1.46e-3) 5.9751e-2 (4.85e-3) 9.3947e-2 (4.89e-4) 4 1.5255e-1 (8.49e-4) 2.3852e-1 (2.01e-2) 1.9833e-1 (5.52e-3) 4.7290e-1 (8.94e-4) 5 2.3283e-1 (4.29e-4) 2.9947e-1 (1.28e-1) 2.9391e-1 (1.04e-1) 5.5282e-1 (4.94e-4) C3_DTLZ4 3 5.1141e-2 (4.32e-4) 7.3349e-2 (2.32e-3) 5.3365e-2 (3.32e-3) 6.5254e-2 (2.35e-4) 4 1.1347e-1 (1.32e-4) 2.5983e-1 (3.25e-2) 1.4651e-1 (2.22e-3) 2.4014e-1 (7.52e-4) 5 2.1698e-1 (3.24e-4) 2.7321e-1 (2.36e-1) 2.6651e-1 (1.52e-1) 4.4652e-1 (5.32e-4)Figure 4
Parallel coordinates figure by four algorithms on the four-objective CF1, CF2, and CF3.(1) MOEA/D-ACO-PBI, (2) MOEA/D-ACOW, and (3) MOEA/D (4) MACS.The conclusions are given below.( 1 ) MOEA/D-ACO-PBI has the best convergence performance, as its approximate fronts cover most of those returned by the other algorithms. MOEA/D-ACOW can obtain well-spread approximately fronts; however, some of its solutions were dominated by MOEA/D-ACO-PBI.( 2 ) It can be observed in Figure 3 that MOEA/D-ACO-PBI dominates in most solutions and that the advantage of MOEA/D-ACO-PBI is more pronounced as the number of objectives increases.


## 6. Conclusions

The major contribution of this paper is the proposal of a MOEA/D-ACO-PBI algorithm for MaOPs. Due to the fact that the number of model optimization objectives increases, the PBI aggregate function was selected to decompose the objectives. MOEA/D-ACO-PBI algorithm can better ensure the convergence of algorithms and the balance of diversity, so as to obtain a better Pareto solution based on the simulation experiment; it could be seen that the results of the MOEA/D-ACO-PBI algorithm are significantly better than those of the other algorithms. Our future work aims at analysing further how to cope with constraint and dynamic constraint processing, and to optimize the number of objectives more effectively.



---

*Source: 5414869-2018-08-30.xml*