(2010) Evolutionary multi-objective optimization algorithms with probabilistic representation based on pheromone trails. In: Proceedings of the 2010 IEEE Congress on Evolutionary Computation (CEC 2010), 18-23

—Recently, the research on quantum-inspired evolutionary algorithms (QEA) has attracted some attention in the area of evolutionary computation. QEA use a probabilistic representation, called Q-bit, to encode individuals in population. Unlike standard evolutionary algorithms, each Q-bit individual is a probability model, which can represent multiple solutions. Since probability models store global statistical information of good solutions found previously in the search, QEA have good potential to deal with hard optimization problems with many local optimal solutions. So far, not much work has been done on evolutionary multi-objective (EMO) algorithms with probabilistic representation. In this paper, we investigate the performance of two state-of-the-art EMO algorithms - MOEA/D and NSGA-II, with probabilistic representation based on pheromone trails, on the multi-objective travelling salesman problem. Our experimental results show that MOEA/D and NSGA-II with probabilistic presentation are very promising in sampling high-quality offspring solutions and in diversifying the search along the Pareto fronts.


I. INTRODUCTION
Evolutionary algorithms are a class of very popular optimization search techniques for tackling hard optimization problems from real-world applications. To design an effective evolutionary algorithm, the choice of proper solution representation and the development of appropriate search operators must be considered. This is because representation of solutions could affect some properties of the search landscap, such as redundancy, neighborhood structure, and ruggedness [1]. Some well-known examples of solution representations include: binary, real-value vector, messy encodings, permutation, and tree structures.
In the past two decades, evolutionary multi-objective (EMO) algorithms have attracted a lot of interest from researchers [2] [3]. Fitness assignment and diversity maintenance are two major research issues in evolutionary multiobjective optimization. However, little work has been done on designing multi-objective oriented representation of solutions. In [4], a hybrid representation was studied for multi-objective optimization. That scheme uses both binary and real-valued representations to encode individuals in the population. In [5], the influence of binary and order-based representations on the performance of EMO algorithms was investigated for multi-objective knapsack problems. Hui Li and Dario Landa-Silva are with the Automated Scheduling, Optimisation and Planning (ASAP) Research Group, in the School of Computer Science, University of Nottingham, United Kingdom (email: hzl@cs.nott.ac.uk and dario.landasilva@nottingham.ac.uk).
Xavier Gandibleux is with the Laboratoire d'Informatique de Nantes Atlantique (LINA), in the Department of Computer Science, University of Nantes, France (email: xavier.gandibleux@univ-nantes.fr).
Currently, the mainstream state-of-the-art EMO algorithms are based on Pareto dominance (e.g. NSGA-II [6] and SPEA-II [7]) or on aggregating functions (e.g. MOGLS [8] and MOEA/D [9]). The research literature shows that Pareto-based EMO algorithms are very effective to tackle continuous multi-objective optimization problems. In contrast, for tackling multi-objective combinatorial optimization problems, EMO algorithms based on aggregating functions seem to be more suitable. The main reason for this is that the later approaches directly use local search to intensify the exploration of promising regions in the search space. Moreover, these algorithms also have advantages in dealing with many-objective optimization problems.
In [10], a quantum-inspired evolutionary algorithm (QEA) was developed to solve combinatorial optimization problems. Unlike other evolutionary algorithms, a QEA uses Q-bit representation to encode individuals. Since each Q-bit individual encodes multiple solutions, a QEA has the ability to provide good diversity in the population. Essentially, each Q-bit individual can be regarded as a simple probabilistic model. Therefore, a QEA is also a kind of multi-model estimation of distribution algorithm (EDA) [11], which samples offspring solutions from probabilistic models. Therefore, there are no crossover and mutation operators used in QEA. In [12], a QEA has been extended to solve the multi-objective knapsack problem.
In order to find a set of diverse non-dominated solutions, many EMO algorithms incorporate mechanisms to encourage diversity of the population in the objective space. However, maintaining diversity in the decision space is also crucial to the performance of EMO algorithms. Based on this, we develop variants of NSGA-II and MOEA/D with probabilistic solution representation. Each probabilistic individual corresponds to one pheromone matrix, which is a probability model commonly used in ant colony optimization (ACO) algorithm [13]. We also compare the performance of NSGA-II and MOEA/D with probabilistic representation to their original versions on the multi-objective travelling salesman problem.
The rest of this paper is organized as follows. Section II introduces some basic definitions in multi-objective optimization. Section III briefly overviews ant colony optimization. Section IV presents the new versions of NSGA-II and MOEA/D with probabilistic presentation. Section V describes the pheromone trails and heuristics in the multiobjective travelling salesman problem. Experimental results are reported and discussed in Section VI while the final section concludes this paper.

II. MULTI-OBJECTIVE OPTIMIZATION
Mathematically, a multi-objective optimization problem (MOP) can be stated as: where x is the vector of decision variables, Ω is the feasible area in the decision search space and F (x) ∈ R m is a vector of objective functions.
When objective functions are conflicting, solutions in the objective space are not completely ordered. For this reason, the optimal solutions of a MOP are trade-offs among objectives, known as Pareto-optimal solutions. For two objective vectors u, v ∈ R m , u is said to dominate v if and only if u i ≤ v i for all i ∈ {1, . . . , m}, and there exists at least one j ∈ {1, . . . , m} satisfying u j < v j . A solution x * ∈ Ω is said to be Pareto-optimal if its objective vector is not dominated by the objective vector of any other solution in Ω. The set of objective vectors of all Pareto-optimal solutions is called Pareto-optimal front.
Traditional methods in mathematical programming often convert a MOP into one single objective optimization problem [14]. Among these traditional methods, weighted sum approach is one of the most commonly used. It minimizes the following weighted scalarizing function: where λ = (λ 1 , . . . , λ m ) is the weight vector with λ i ≥ 0 for all i ∈ {1, . . . , m} and m i=1 λ i = 1. The optimal solution of g is called supported. Under some mild conditions, it is also Pareto-optimal to the MOP in (1).

III. ANT COLONY OPTIMIZATION
Ant colony optimization (ACO) [13] is a constructive meta-heuristic that mimics the behavior of ants when seeking a path between the ant colony and a source of food. In ACO, the priori information (i.e. problem-specific heuristic) on the structure of good solutions and the posteriori information (i.e. artificial pheromone trail) on the probabilistic distribution of previously obtained good solutions are used to sample offspring solutions in promising areas of the search space. Over the years, ACO algorithms have been applied to solve many combinatorial optimization problems [15] [16], such as travelling salesman problem and quadratic assignment problem.
Consider a combinatorial optimization problem (S, f, Ω) with the following characteristics: • C = {c 1 , . . . , c Nc } is a finite set of solution components; • each state x of the problem can be represented by a sequence of solution components c i , c j , . . . , c k , . . . . The set of all states is denoted by S, the set of candidate solutions.
• Ω ⊂ S is the feasible set of the search space satisfying certain constraints and f (s) is the cost function for each candidate solution s ∈ S.
In ACO, artificial ants build solutions by incrementally adding solution components to a partial solution, which is initially empty. This procedure is stopped until a complete solution is generated. The decision rule for accepting one component as new member of the partial solution is determined by both heuristic information and pheromone trail. The ACO process for constructing solutions can also be understood as moving on a graph G = (C, L), where the vertices are the components in C and L = {l ij } is the connection between components. Connections l ij are associated with both pheromone trails τ ij , storing global search information obtained from previous search, and heuristic values η ij , representing problem-specific greediness. Assume that x r =< x r−1 , c i > is the partial solution constructed so far and c i is the current node, the next node c j is chosen from the feasible neighborhood N i of node c i by applying some probabilistic decision rules. Then, the new partial solution x r+1 is set to < x r , j >. The decision rule is determined by pheromone trails and heuristic values.
In the past twenty years, many variants of ACO algorithms have been proposed. The main differences among them lie in the design of pheromone representation (i.e. how to define solution components) and probabilistic decision rules. For example, in the first implementation of ACO algorithm for the travelling salesman problem, called ant system (AS) [15], the probability of moving from node c i to node c j was defined as: where • η ij = 1/d ij is the priori heuristic value and d ij is the distance between cities i and j; • α and β are two parameters representing the relative importance of heuristic information and pheromone information. • N i is the feasible neighborhood of city i not visited yet. When α = 0, the closest city is more likely to be selected as the new component of the current partial solution.
Once a complete tour x is generated, the pheromone trails need to be updated in two steps. First, the pheromone trail for each edge is modified by where ρ ∈ (0, 1] is the evaporation rate. Second, for each edge (i, j) belonging to the tour x, the related pheromone trail is updated as follows: where ∆τ ij (t) = 1/L(t)(L(t) is the length of x).
In [17], one of the successors of AS, called Ant Colony System (ACS), used a different probability decision rule. At each step, the node j with the maximal greediness was chosen with probability q 0 ∈ [0, 1]. That is, The ACS approach uses the same rules as the original AS with probability set to 1 − q 0 . Like evolutionary algorithms, the application of ACO algorithms for multi-objective optimization has also attracted interest for the research community. In [18], an empirical analysis of multi-objective ACO algorithms was conducted on bi-objective TSP problems.

IV. MOEA/D AND NSGA-II WITH PROBABILISTIC REPRESENTATION
In this section, we first describe a quantum-inspired genetic algorithm (QEA), in which the probabilistic representation was first used. Then, the modified versions of MOEA/D and NSGA-II with probabilistic representation are presented.

A. QEA
QEA is an evolutionary algorithm paradigm based on the concepts of qubits and superposition of states in quantum mechanics. In this algorithm, solutions are encoded using a Q-bit representation system. Each Q-bit corresponds to two states '0' or '1'. The probabilities of finding these two states are denoted by α 2 and β 2 respectively. Normalization of states needs to guarantee α 2 + β 2 = 1. An individual in a n Q-bit system can be defined as: where n is the number of Q-bits, and (α i , β i ), i = 1, . . . , n is a pair of complex numbers specifying the probability amplitudes of states '0' and '1' of the i-th Q-bit, respectively. Essentially, each Q-bit individual is a simple probabilistic model, which can represent multiple states of Q-bits. In QEA, a population of probabilistic Q-bit individuals is maintained. It can also be viewed as a multi-modal estimation of distribution algorithm (EDA). Like classical EDA algorithms, QEA produces offspring solutions by sampling from a probabilistic model, i.e., observing Q-bit probabilistic individuals. Since each Q-bit individual stores global statistical information of good solutions obtained previously in the search, evolutionary algorithms with Q-bit representation have advantage in exploring promising areas of the search space. In [12], a multi-objective version of quantum-inspired evolutionary algorithm, called QMEA, was also proposed. That algorithm uses NSGA-II as its baseline algorithm and samples new solutions by observing Q-bit individuals.

B. MOEA/D
In [9], Zhang and Li proposed a decomposition-based evolutionary multi-objective optimization algorithm, called MOEA/D. It decomposes a multi-objective optimization problem into NP single objective subproblems with objective functions g(x, λ (i) ), i ∈ {1, . . . , NP}. The main goal of MOEA/D is to find the optimal solutions of NP subproblems simultaneously. The information from the optimization of one subproblem is used to assist the optimization of similar subproblems. MOEA/D provides a general framework, which Step 1: Initialization Step 1.1: Initialize NP weight vectors λ (s) , s = 1, . . . , NP and set t := 0. For each weight vector λ (s) , calculate its neighborhood B(s), which includes the indexes of K closest weight vectors to λ (s) ; Step 1.2: Initialize pheromone trails τ  Step 2.2: Improve y by applying local search; Step 2.3: Update EP with y; Step 2.4: Compete y with all neighbors.
Step 4: Termination If the stopping conditions are satisfied, then output EP and stop; otherwise set t := t + 1 and go to Step 2.
In the previous implementation of MOEA/D, each subproblem is only associated with one individual, which is its current best solution. As discussed in [19], the competition between neighboring solutions might cause loss of diversity. To avoid this weakness, more complex memory should be used to store the information about the history of the search for each subproblem. This memory could be a small population or a probability model. So far, no work has yet been done on MOEA/D with probabilistic representation. In this paper, we investigate a new version of MOEA/D, in which each subproblem is associated with one probabilistic individual represented by pheromone trails. At each generation t, MOEA/D maintains the following data structures: x (s) (t), s = 1, . . . , NP is the current best solution to g(x, λ (s) ); • Pheromone population: {τ  ij (t) is the pheromone trail for the subproblem associated with g(x, λ (s) ).
• Subproblem heuristics: {η is the heuristic for the subproblem associated with g(x, λ (s) ). • External population EP which is used to store all nondominated solutions found during the search. The algorithmic framework of MOEA/D with probabilistic representation is described in Procedure 1. Step

C. NSGA-II
Step 1: Initialization Step 1.1: Initialize NP weight vectors randomly and set t := 0; Step 1.2: Initialize {τ  Step 3.2: Copy the set Q of the best NP solutions into the population into the next generation; Step 4: Pheromone Update for s = 1 to NP Step 4.1 Evaporate pheromone trails τ  Compared to MOEA/D, NSGA-II has no bias in searching any particular part of the Pareto front. All non-dominated solutions in the current population have equal chance to be selected for reproduction. However, this might not be efficient when sampling offspring solutions due to the following two reasons. First, the non-dominated solutions might have very different structures in the decision space. Therefore, the possibility of generating high-quality offspring solutions by recombining these solutions is low. Second, the design of recombination operators is often problem-dependent. Efficient recombination operators for some combinatorial optimization problems are not always readily available.
Here, we also investigate the version of NSGA-II with probabilistic representation, in which multiple probability models (pheromone trails) are maintained. Each probability model stores long-term search information. All offspring solutions are sampled from these probability models. Similar to MOEA/D in Procedure 1, NSGA-II also needs to maintain a set of weight vectors, current population, pheromone population, subproblem heuristics and external population. The details of NSGA-II with pheromone representation are described in Procedure 2.
Steps 1.1 and 1.2 initialize pheromone population, current population, and subproblem heuristics. These are similar to the related steps in MOEA/D. Step 2.1 samples offspring solutions from a probabilistic individual. Unlike the original NSGA-II, tournament selection and crossover are not needed here. In Step 2.2, a local search is applied to improve offspring solutions. The acceptance function used in local search is a scalar function g(x, λ) with normalized weight vector generated randomly. The external population is updated in Step 2.3, which use the same rules as in MOEA/D. Step 3 ranks all solutions in the union of current and offspring populations by applying non-dominated sorting. The best NP solutions are copied into the population Q in the next generation. In Step 4.1, the pheromone trails of each probabilistic individual are evaporated. Then, they are updated by the offspring solution generated by this probabilistic individual as well as those members of Q that dominate the offspring solution.

V. PHEROMONE TRAILS AND HEURISTICS IN MULTI-OBJECTIVE TRAVELLING SALESMAN PROBLEM
In this paper, we tested MOEA/D and NSGA-II with probabilistic representation on the multi-objective travelling salesman problem. A travelling salesman problem with m objectives can be formulated as: where • π = (π 1 , . . . , π n ) is a permutation of {1, . . . , n} where n is the number of cities.
πiπj is the cost between city π i and city π j regarding criterion k. Like the single objective ACO algorithms for travelling salesman problem, solution components correspond to cities (nodes) in a graph while links between components correspond to edges between cities. In MOEA/D, the heuristic values of each subproblem g(x, λ (s) ) are defined as The definition of heuristic values in NSGA-II is almost the same as that in MOEA/D. The difference is that NSGA-II uses weight vectors generated randomly and changed at each construction step. In this paper, only bi-objective TSP instances are tested. MOEA/D uses a set of NP weight vectors with uniform distribution, which are generated by λ The total amount of pheromone values changed is: where S is the set of solutions selected for pheromone update.

VI. COMPUTATIONAL EXPERIMENTS
In this section, we compare the original MOEA/D and NSGA-II to their variants with probabilistic representation based on pheromone trails, which are denoted by MOEA/D-ACO and NSGA-II-ACO respectively.

A. Experimental Settings
In our experiments, we used two 50-city TSP instances (KROAB50 and KRCD50) as well as two 100-city TSP instances (KROAB100 and KROCD100) [22]. All four instances are bi-objective. They are constructed by combining two benchmark single objective TSP instances. The ACS   TABLE I  AVERAGE NUMBER OF FUNCTION EVALUATIONS USED BY MOEA/D   AND NSGA-II WITHOUT LOCAL SEARCH   Instance  MOEA/D  NSGA-II  MOEA/D-ACO  NSGA-II-ACO  KROAB50  1487523  825694  95985  67450  KROCD50  1489063  862066  95948  68460  KROAB100  2395545  1449649  49068  43800  KROCD100  2379413  1437141  49218  43610 probability decision rules are used in both algorithms. We also compared them with the original versions of MOEA/D and NSGA-II, in which cycle crossover is used for reproduction. The population size is set to 200 in all algorithms for all instances. Both α and β used in probabilistic decision rules are set to 1. The probability q 0 in ACS decision rule is 0.95 while the pheromone evaporation rate is 0.1. The neighborhood size in MOEA/D is 20. We also investigated the performance of MOEA/D and NSGA-II with local search. The 2-opt local move, which randomly exchanges two edges in the tour, is used to generate neighboring solutions. Each local search procedure is stopped after examining 100 neighbors. The total number of runs of each algorithm for each instance was set to 20. In each run, each algorithm is stopped after 50 seconds for the instances with 50 cities and 100 seconds for the instances with 100 cities. All algorithms were implemented in C++ and run in PC computer (Intel (R) Core (TM)2 CPU, 1.86 GHz, 2GB RAM) running Windows XP.
To measure the quality of the non-dominated solutions found by MOEA/D and NSGA-II, we use the inverted generation distance (IGD) metric in our experiments [9], which can be formulated as: IGD(P, P * ) =  Table I. Table II shows the mean and standard deviation IGD-metric values of non-dominated solutions found by the four algorithms. The average differences between the best and the worse values of the non-dominated solutions obtained by the four algorithms on each objective are summarized in Table III. From these results, we can make the following observations.
• Table I    NSGA-II. These results also indicate that the problemspecific constructive heuristics are more efficient than the standard crossover and mutation operators when sampling offspring solutions. • It can be seen from Fig. 1 and 2 Table II show the MOEA/D-ACO and NSGA-II-ACO found better IGDvalues than MOEA/D and NSGA-II. These results show that the non-dominated solutions found by two EMO algorithms with probabilistic representation are closer to the Pareto fronts than those found by the original versions. We can also observe that NSGA-II-ACO performed better than MOEA/D-ACO on all 4 instances in terms of the IGD-metric. • The results in Table III show that MOEA/D-ACO performed better than NSGA-II-ACO in diversity. The ranges of non-dominated fronts found by MOEA/D-ACO are wider than those of non-dominated solutions obtained by NSGA-II-ACO. The reason is that MOEA/D-ACO specifically spends equal computational effort in searching different parts of the Pareto fronts. In contrast, NSGA-II-ACO does not have any bias to search any part of the Pareto fronts.
2) MOEA/D and NSGA-II with Local Search: Fig. 3 and 4 plot the non-dominated solutions found by the four algorithms using local search in 20 runs. The average numbers of solutions examined by each algorithm are reported in Table  IV. The mean and standard deviation IGD-metric values of non-dominated solutions found by the four algorithms are shown in Table V. The average differences between the best and the worse values of the non-dominated solutions obtained by the four algorithms for each objective are summarized in Table VI. From these results, the following observations can be made.
• The results in Fig. 3 and 4 show that all four algorithms with local search are able to find non-dominated solutions close to the Pareto fronts. Again, NSGA-II-ACO performed better than MOEA/D-ACO in approximating the middle parts of the Pareto fronts. From the subfigures, we can see that NSGA-II performed better than MOEA/D in convergence towards the middle parts of the Pareto fronts no matter if probabilistic representation is used or not. However, it is also evident that MOEA/D has better performance in diversity. • The mean and standard deviation of IGD-metric values in Table V show that MOEA/D outperformed NSGA-II on KROAB50, KROCD50, and KROAB100, and performed slightly worse on KROCD100 in minimizing the closeness to the Pareto fronts. This is due to the advantage of MOEA/D in diversifying the search along the Pareto fronts. We can also note that MOEA/D-ACO is inferior to MOEA/D on all instances in terms of the IGD-metric. This suggests that MOEA/D with local search can have excellent performance in converging towards the Pareto fronts. However, it needs to consume more number of function evaluations (see Table IV). • Again, the differences between the best and worse values of the non-dominated solutions show the MOEA/D-ACO found wider ranges of non-dominated fronts than MOEA/D. This indicates that the use of probabilistic representation in MOEA/D is beneficial in diversifying the search along the Pareto fronts.

VII. CONCLUSION
In this paper, we investigated the performance of two stateof-the-art evolutionary multi-objective (EMO) algorithms, namely MOEA/D and NSGA-II, when using probabilistic representation based on pheromone trails. In both algorithms, an individual is encoded by a probability model, i.e. pheromone trails. Offspring solutions are sampled from the probabilistic individuals and problem-specific heuristics. We compared the four algorithms without local search -MOEA/D, NSGA-II, MOEA/D-ACO and NSGA-II-ACO on the multi-objective travelling salesman problem. Our experimental results showed that MOEA/D-ACO and NSGA-II-ACO performed much better than the original MOEA/D and NSGA-II without probabilistic representation. MOEA/D-ACO performed worse than NSGA-II in finding solutions in the middle parts of Pareto fronts but it found more diverse Pareto fronts. We also compared the four algorithms when using local search. Our results showed that MOEA/D found better IGD-metric values than MOEA/D-ACO. However, it needs more number of function evaluations. We conclude that probabilistic representation can improve the efficiency of EMO algorithms in sampling high-quality solutions by constructive heuristics and in diversifying the search along Pareto fronts by pheromone trails. In the future, we intend to apply EMO algorithms with probabilistic representation to other challenging multi-objective combinatorial optimization problems.