# An Enhanced Differential Evolution Based Algorithm with Simulated Annealing for Solving Multiobjective Optimization Problems

Authors: Bili Chen; Wenhua Zeng; Yangbin Lin; Qi Zhong
Journal: Journal of Applied Mathematics
 (2014)

---

## Abstract

An enhanced differential evolution based algorithm, named multi-objective differential evolution with simulated annealing algorithm (MODESA), is presented for solving multiobjective optimization problems (MOPs). The proposed algorithm utilizes the advantage of simulated annealing for guiding the algorithm to explore more regions of the search space for a better convergence to the true Pareto-optimal front. In the proposed simulated annealing approach, a new acceptance probability computation function based on domination is proposed and some potential solutions are assigned a life cycle to have a priority to be selected entering the next generation. Moreover, it incorporates an efficient diversity maintenance approach, which is used to prune the obtained nondominated solutions for a good distributed Pareto front. The feasibility of the proposed algorithm is investigated on a set of five biobjective and two triobjective optimization problems and the results are compared with three other algorithms. The experimental results illustrate the effectiveness of the proposed algorithm.


---

## Body

## 1. Introduction

Many real world problems are MOPs and it has prompted a wide research boom about the MOPs over the past few decades. A lot of multiobjective evolutionary algorithms (MOEAs) have been suggested such as these famous algorithms in [1–4]. For comprehensive overviews, the reader can refer to [5–7].In the present study, we emphasize on the differential evolution (DE) for solving MOPs. DE, one of the most popular evolutionary algorithms, was initially presented by Storn and Price [8, 9]. DE has a simple principle and is easy to be implemented. DE was originally used for solving single objective optimization problems in [10]. In 2001, Abbass et al. [11] presented a Pareto-frontier differential evolution algorithm (PDE) for solving MOPs. In 2002, Abbass [12] further improved the PDE with self-adaptive crossover and mutation operator. In the same year, the nondominated sorting and the concept of ranking used in NSGA-II are incorporated in the DE to solve MOPs in [13]. Xue et al. [14] proposed an algorithm called Pareto based multiobjective differential evolution (MODE), which incorporates the concept of Pareto optimal into the mutation operator. Parsopoulos et al. [15] presented a vector evaluated differential evolution (VEDE) for solving MOPs. A domination selection operator is adopted to improve the performance of the algorithm. Kukkonen and Lampinen [16, 17] proposed a generalized differential evolution (GDE) and an improved version of GDE for solving MOPs with constraints. In 2006, Hernández-Díaz et al. [18] combined the DE and the rough set theory for MOPs. Wang et al. [19] proposed a self-adapted differential evolution algorithm, which adopts an external elitist and a crowding entropy diversity measure tactic. Huang et al. [20, 21] improved their self-adaptive differential evolution (SaDE) to solve MOPs by a multiobjective adaptive differential evolution (MOSaDE). Moreover, they further extended MOSaDE by objective-wise learning strategies and called it as WO-MOSaDE. Ali et al. [22] proposed a modified differential evolution (MDE) for solving MOPs. In 2012, they [23] further continued their research by giving an enhanced version of MDE, named multiobjective differential evolution algorithm (MODEA), which incorporates the opposition-based learning [24] and the concept of random localization in mutation.Besides the DE method, in this paper, the concept of simulated annealing is utilized for controlling the acceptance of the candidate solutions. Simulated annealing, proposed by Kirkpatrick et al. [25] is a popular probabilistic metaheuristic for solving the optimization problems. Simulated annealing is inspired by the annealing in metallurgy, a technique of controlling the temperature of a material from a high temperature to low one in a probability for increasing the volume of its crystals and reducing their defects. The simulated annealing algorithm works for finding minimal cost solutions by minimizing the associated energy function. In the process, the temperature is first in a high level (T
max
⁡) and then decreases gradually, if the cooling is sufficiently slow, the global minimum will be reached. The viewpoint has been proven by S. Geman and D. Geman in [26]. There are some attempts of adopting simulated annealing to solve MOPs. Especially, there have been a few studies that utilize the concept of Pareto-dominance in the multiobjective simulated annealing, the reader can refer to [27] for a review about this research. In [28], an archived multiobjective simulated annealing algorithm named AMOSA is proposed. The algorithm incorporates the concept of archive and adopts a domination status to determine the acceptance of a new solution.The rest of the paper is organized as follows. In Section2, the problem formulation is introduced. The proposed algorithm is presented in Section 3. In Section 4, some numerical examples are given to show the performance of the proposed algorithm with the other algorithms, while the conclusion is reached in Section 5.


## 2. Problem Formulation

We assume that the MOPs given in this paper are minimization problems, while each maximization problem can be transformed into a minimization problem. A MOP can be formally described as follows:(1)
min
⁡
F
(
x
)
=
(
f
1
(
x
)
,
f
2
(
x
)
,
…
,
f
m
(
x
)
)
(2)
s
.
t
.
g
i
(
x
)
≤
0
,
i
=
1
,
…
,
p
(3)
h
i
(
x
)
=
0
,
i
=
1
,
…
,
q
,
where m is the number of objective functions and f
i
(
x
) is the ith objective function to be optimized. g
i
(
x
) and h
i
(
x
) are the set of inequality and equality constraints, respectively. p is the number of inequality constraints and q is the number of equality constraints. x
(
x
1
,
x
2
,
…
,
x
n
) is the vector of variables which is called a solution for a MOP, where n is the number of variables. That is to say, a MOP is to find all x
(
x
1
,
x
2
,
…
,
x
n
) which satisfies (2) and (3).Suppose that there are two solutionsx
(
x
1
,
x
2
,
…
,
x
n
) and y
(
y
1
,
y
2
,
…
,
y
n
), x is said to dominate y if ∀
i
∈
{
1,2
,
…
,
m
}, f
i
(
x
)
≤
f
i
(
y
) and∃
i
∈
{
1,2
,
…
,
m
}, f
i
(
x
)
<
f
i
(
y
). x dominates y which can be denoted as x
≺
y. A solution x is called a Pareto-optimal solution if there is not another solution x
′ satisfying x
′
≺
x. A set containing all the Pareto-optimal solutions is called the Pareto-optimal set (PS). A set containing all the Pareto-optimal objective vectors is called the Pareto front (PF) such that PF
=
{
F
(
x
)
∣
x
∈
PS
} [29].


## 3. The Proposed Algorithm

The proposed algorithm is an extension of MODEA algorithm proposed by Ali et al. [23]. In the proposed MODESA, we focus on improving the convergence of the obtained solution towards the Pareto-optimal front, and the diversity and distribution of the solutions. Besides using the simulated annealing process to control the acceptance of the candidate solutions, an efficient diversity maintenance mechanism [30] is also adopted to improve the distribution of these obtained solutions.


### 3.1. Population Initialization

In MODESA, initial population is generated in the same way as that of MODEA [23]. Because of the lack of a prior information about the solution, its corresponding opposite estimate is considered for getting a better solution. Here, the concept of opposite solution is given. Suppose that the current solution is x
(
x
1
,
x
2
,
…
,
x
i
,
…
,
x
n
) and the lower and upper bounds of x
i are v
i
,
min
⁡ and v
i
,
max
⁡, respectively. Then the opposite solution x
′
(
x
1
′
,
x
2
′
,
,
…
,
x
i
′
,
…
,
x
n
′
,
) of x is calculated as x
i
′
=
v
i
,
min
⁡
+
v
i
,
max
⁡
-
x
i [23].The advantage of doing this is that we can choose the better one between the solution and the opposite one as the initial solution. The procedure of initialization in the paper is given as follows. A population withN solutions is generated using uniform random distribution and a population containing their opposite solutions is constructed too. Next, a combined population with the population and the opposite population is constructed and the best N solutions of the combined population are selected by using the nondominated sorting in NSGA-II [2] and the diversity maintenance mechanism based on vicinity distance [30].


### 3.2. Differential Evolution (DE)

Differential Evolution (DE) presented by Storn and Price [8, 9] is an easy and efficient evolutionary algorithm. DE maintains a population of N solutions and contains three main operators: mutation, crossover, and selection.The mutation operation is defined as follows:(4)
z
=
x
r
1
+
f
*
(
x
r
2
-
x
r
3
)
,
where x
r
1, x
r
2, and x
r
3 are randomly selected from the current population and r
1
≠
r
2
≠
r
3
≠
i. f is a control parameter which is set in (
0,2
] in Storn and Price [8, 9]. For the sake of convenient description, the solution which is selected to plus a perturbed value is called as the base vector. x
r
1 (the base vector) is perturbed by adding to it the product of the control parameter and the difference of x
r
2 between x
r
3. z
(
z
1
,
z
2
,
…
,
z
n
) is the generated perturbed solution. It is important to note that if z
i is out of bounds, it must be repaired. As mentioned above, the lower and upper bounds of variable z
i are v
i
,
min
⁡ and v
i
,
max
⁡, respectively. In this study, if z
i is less than v
i
,
min
⁡, z
i is set to v
i
,
min
⁡; if z
i exceeds v
i
,
max
⁡, z
i is set to v
i
,
max
⁡. Ali et al. [23] have proved that using the one having the best fitness value in these three randomly selected individuals as the base vector has two advantages. The first advantage is that it is neither merely random nor merely greedy. The other one is that it can offer a localized effect to help in finding more different regions of the search space around the underlying solutions [23]. Therefore, in MODESA, after choosing the three individuals, we select the best one of them to be the base vector.After the perturbed solution is generated, a trial solutionw is produced by the crossover operation which can be defined as follows. For each variable j,
(5)
w
j
=
{
z
j
,
rand
j
≤
cr
∨
j
=
j
j
x
j
,
otherwise
,
where j
j
∈
{
1
,
…
,
n
} is a random number and cr
∈
[
0,1
] is the crossover rate. In general DE, the selection operator is utilized to decide whether the trial solution w replaces the target solution x
.The selection operation can be defined according to the following:
(6)
x
=
{
w
,
w
≺
x
x
,
otherwise
.In the proposed MODESA, the target individualx is not simply replaced by the trial solution w. A simulated annealing procedure introduced in Section 3.3 is employed to determine the acceptance of the trial solution.


### 3.3. Simulated Annealing Procedure

After the mutation phase, the crossover operator is performed, and the algorithm enters the selection phase. In MODESA, we incorporate the principle of simulated annealing into the selection phase. Before introducing the proposed selection phase, we first give a concept of domination in simulated annealing.


#### 3.3.1. The Concept of Domination

To calculate the acceptance probability of a new solution, the amount of domination used in AMOSA [28] is adopted here. Given two solutions x and y, the amount of domination is calculated by the following equation:
(7)
Δ
dom
⁡
a
,
b
=
∏
i
=
1
,
f
i
(
x
)
≠
f
i
(
y
)
m
(
|
f
i
(
x
)
-
f
i
(
y
)
|
R
i
)
,
where m is the number of objectives and R
i is the range of the ith objective. When R
i is unknown, the solutions in the current population are used for calculating it.


#### 3.3.2. Simulated Annealing Used in the Selection Operator of the Proposed MODESA

As (6), if the trial solution dominates the current solution, the general DE simply replaces the current solution by the trial solution. However, doing this may make the algorithm over converge and easy to trap into local optimal. As is known, there are three possibilities when comparing two solutions, the current solution x and the new solution y, as follows:(1)
xdominates y;
(2)
x is dominated by y;
(3)
x and y are incomparable.Ifx dominates y, it is obvious that y is not better than x, so that we can simply replace y by x. If x is dominated by y, we accept the new solution y distinctly. If x and y are incomparable, the solution y may have the potential to guide the solutions towards an unreached area of the search space. Therefore, in the DE process, we use the simulated annealing to accept y under a certain probability. Meanwhile, we need to control the acceptance probability, and the reason is that if we do not control the acceptance probability, all these potential solutions which are not absolutely good will enter to the next generation which will confuse the algorithm and make it unable to have a right direction towards the Pareto-optimal front. That is to say, we must decrease the acceptance probability. In simulated annealing, with the decreasing of the temperature, the acceptance probability is decreased too, which will avoid the algorithm producing too much potential solutions in the end. It is worth mentioning that if we allow some potential solutions existing during the procedure, the next nondominated sorting procedure may again remove these potential solutions, which will make the purpose of allowing some potential solutions to search more areas of the searching space fail. Therefore, in the proposed MODESA, each individual is assigned a prior life cycle. The number of prior life cycles means how many times the individual can be selected preferentially in the nondominated sorting. For the sake of convenience, the life cycle value of solution x
i is denoted as l
i
f
e
i. Once the individual is selected preferentially in the nondominated sorting, the prior life cycle of the individual is minus one when it is greater than zero. Note that, in the population initialization, the prior life cycle of each individual is set to zero.In this paper, given a temperaturetemp, the new solution is selected with a probability prob as follows:
(8)
prob
=
exp
⁡
⁡
(
-
Δ
dom
⁡
avg
temp
)
,
(9)
Δ
dom
⁡
avg
=
∑
i
=
1
k
Δ
dom
⁡
i
,
new-solution
k
,
where k is the number of solutions that dominates the new solution. It is worth mentioning that Δ
dom
⁡
avg represents the average amount of domination of the new solution and the k solutions that dominate it. Suppose that Δ
dom
⁡
avg is constant, the acceptance probability decreases when the temperature decreases. The simulated annealing used in the selection operator of MODESA is as shown in Algorithm 1.Algorithm 1:The proposed simulate annealing algorithm.
input: x
i, w
i, t
e
m
p
e
r
a
t
u
r
e, T
min
⁡, T
max
⁡, P
t, P
temp, m
a
x
_
l
i
f
e
_
c
i
c
l
e
If (t
e
m
p
e
r
a
t
u
r
e
>
T
min
⁡) {
Compute the value of Δ
dom
⁡
avg by (9).
Calculate the acceptance probability prob:
prob
=
ex
p
⁡
(
-
Δ
dom
⁡
avg
/
t
e
m
p
e
r
a
t
u
r
e
)
With a probability prob Do
Put the current individual x
i to the temporary population: P
temp
=
P
temp
∪
x
i.
Replace the current individual x
i by the trial individual w
i: x
i
=
w
i.
Assign the life value m
a
x
_
l
i
f
e
_
c
y
c
l
e to the new individual x
i:  l
i
f
e
i
=
m
a
x
_
l
i
f
e
_
c
y
c
l
e,
wherel
i
f
e
i represents the life cycle value of solution x
i.
otherwise
Put the current individual x
i to the temporary population: P
temp
=
P
temp
∪
x
i.
}// End If


### 3.4. Diversity Maintenance Mechanism

The diversity maintenance mechanism adopted in the paper is presented by Kukkonen and Deb [30] for pruning of nondominated solutions. Instead of using the crowding distance in Deb et al. [2], the method proposed a new crowding estimation called vicinity distance using nearest neighbors. In this paper, we call this method as the vicinity distance based pruning method for convenience. The basic process of the method is to remove the most crowded solutions in the nondominated solution one by one. Moreover, after each removal, the vicinity distance values of the remaining solutions must be updated again. The vicinity distance is used to get an estimate of the density of solutions surrounding a particular solution, it calculates the product of k nearest neighbors of the particular solution. The solutions with the smallest value of vicinity distance are considered to be the most crowded ones. In [30], the value of k is assigned to the number of objectives; that is, k
=
m. The reader can refer to [30] to get more information about this method.


### 3.5. The Framework of the Proposed MODESA Algorithm

#### 3.5.1. The Framework of the Proposed MODESA Algorithm

As mentioned before, the proposed algorithm first uses the population initialization phase to generate a population withN solutions. After the nondominated sorting, each solution has a ranking value. Next, the main loop begins. We introduce the tth generation of MODESA in the following.At thetth generation, an empty temporary population is created first. The temporary population is going to store the new good solutions generated by the DE in the loop. Then for each individual in the current population, the DE operators are used to generate a trial individual. Next, there are three possibilities when comparing the current individual and the trail individual. If the trial individual dominates the current individual, we simply replace the current solution by the trial solution. If the trial individual and the current solution are incomparable, the proposed simulated annealing process is adopted to determine whether the trial individual is accepted. If the trial individual is dominated by the current solution, we also put the trial individual into the temporary population as MODEA does. After all the individuals are handled, the current population and the temporary population are combined as a union population. Note that, if the current generation is the last generation, the members that are dominated by any other individuals must be removed in the union population. The reason of doing this is that, because we allow some potential solutions which are not absolutely good and have a positive life cycle to exist in the population; therefore, at the end of the generation, these solutions must be removed or they will influence the effectiveness of the results. Finally, the nondominated sorting and the vicinity distance based pruning approach are used for getting the best N solutions.In this paper, the parameters that need to be set a priori are listed in Table1 and the Pseudocode of MODESA is shown in Pseudocode 1.Table 1
The parameters used in the proposed algorithm.
Parameter's name
Description
Value
N
The population size
100
m
Number of objectives
As the test problem
n
Number of variables in a solution
As the test problem
T
Maximum number of generations
250
cr
The crossover rate used in DE
0.3
f
The control parameter used in DE
0.5
T
max
⁡
Maximum temperature in simulated annealing
100
T
min
⁡
Minimal temperature in simulated annealing
1
e
-
7
γ
The cooling rate in simulated annealing
0.6
m
a
x
_
l
i
f
e
_
c
i
c
l
e
Maximum prior life cycle value of the individual
1Pseudocode 1:Pseudocode of MODESA.
Step  1. Initialization
Initialize the values of the parameters in Table1.
Generate the initial populationP
0 using uniform distribution.
GenerateN opposite solutions denoted as P
obl by using OBL operation on P
0.
SelectN best individuals as the current population by using non-dominated sorting and the
vicinity distance based pruning method onP
0
∪
P
obl.
For (i
=
0; i
<
N; ++i) {
l
i
f
e
i
=
0.
}
/
/ End For
t
=
0
Step  2. Main loop
While (t
≤
T) {
t
e
m
p
e
r
a
t
u
r
e
=
T
max
⁡.
Construct a temporary population which is denoted asP
temp and set P
temp
=
∅.
For (i
=
0; i
<
N; ++i) {
Randomly select three distinct individuals:x
1, x
2 and x
3, which are different
from the current individualx
i from the current population.
Select the best individualx
best from the three distinct individuals as the base vector
for the next mutation operation.
Produce a trial individualw
i:
generate a random numberj
j
∈
{
1
,
…
,
n
}
For each variablej of w
i
{
With a probability cr or j is equal to j
j
do
w
j
,
i
=
x
j
,
1
+
f
*
(
x
j
,
2
-
x
j
,
3
).
otherwise
w
j
,
i
=
x
j
,
best.
If (w
j
,
i
>
v
j
,
max
⁡) w
j
,
i
=
v
j
,
max
⁡.
If (w
j
,
i
<
v
j
,
min
⁡) w
j
,
i
=
v
j
,
min
⁡.
}// End For
Check the domination status of the trial individual w
i and the current individual x
i.
If (w
i dominates x
i) {
Replace the current individual x
i by the trial individual w
i: x
i
=
w
i.
}
Else if (w
i and x
i are incomparable) {
Use the proposed simulated annealing to allocate the two individuals: w
i and x
i.
SimulateAnnealing (x
i, w
i, t
e
m
p
e
r
a
t
u
r
e, T
min
⁡, T
max
⁡, P
t, P
temp,  m
a
x
_
l
i
f
e
_
c
i
c
l
e).
t
e
m
p
e
r
a
t
u
r
e
=
γ
*
t
e
m
p
e
r
a
t
u
r
e.
}
Else if (w
i is dominated by x
i) {
Put the trial individual w
i to the temporary population.
}// End If
}// End For
Union the two population:P
t
∪
P
adv.
If (t
=
=
T) {
Use the non-dominated sorting to the union population to remove the individuals that
are dominated by any other solutions.
}// End If
SelectN individuals as the next population from the union population by using the
non-dominated sorting.
For (i
=
0; i
<
N; ++i) {
If (l
i
f
e
i
>
0)
l
i
f
e
i
-
-;
}// End For
Apply the vicinity distance based pruning method.
}// End While


#### 3.5.2. The Computational Complexity Analysis

The proposed MODESA approach is a simple algorithm to be implemented. In the following, the overall computational complexity of MODESA is analyzed. Basic operations and their worst case complexities are given as follows.(1)
The initialization phase is selectingN solutions out of 2
N solutions by using nondominated sorting and the vicinity distance base pruning method: O
(
m
×
4
N
2
)
+
O
(
m
2
×
2
N
log
⁡
⁡
(
2
N
)
).
(2)
The procedure of mutation isO
(
m
×
2
N
).
(3)
The procedure of checking the domination status of the trial solution and the current solution isO
(
m
×
N
), and the procedure of simulated annealing is O
(
(
m
+
m
×
N
)
×
N
).
(4)
The procedure of selectingN solutions out of 2
N solutions by using nondominated sorting and the vicinity distance base pruning method: O
(
m
×
4
N
2
)
+
O
(
m
2
×
2
N
log
⁡
⁡
(
2
N
)
).Here,m and N represent the number of objectives and the population size, respectively. We can see that the overall worst case complexity of the MODESA is O
(
m
2
N
log
⁡
⁡
N
+
m
N
2
).


#### 3.5.3. The Difference between MODESA and MODEA

Considering that the proposed algorithm MODESA is an improvement of MODEA, it is necessary to provide the difference between MODESA and MODEA. There are two main differences between MODESA and MODEA as follows.(1)
When comparing the new solution and the current solution, there exist two processing ways according to two different results of comparison. If the new solution dominates the current solution, the MODEA replaces the current solution with the new solution and puts the current solution into the temporary population; otherwise, the new solution is taken into the temporary population. After each solution in the population is handled, the two populations are combined. It indicates that the new solutions are all retained. However, MODESA works in three kinds of ways according to different results. If the new solution dominates the current solution, the current solution is replaced with the new solution and the current solution is removed directly, which means that if the current solution is worse than the new solution, it would not be retained. If the new solution and the current solution are incomparable, the proposed simulated annealing algorithm is adopted to determine the acceptance of the current solution; moreover, the usage of a prior life cycle enables some potential solutions preferentially to be selected in the next generation. Otherwise, the new solution is added to the temporary population and the current solution is retained in the current population.
(2)
The diversity maintenance mechanism is different. In MODEA, the diversity preserving strategy based on crowding distance in NSGA-II is utilized, while in MODESA, the vicinity distance based pruning method in [30] is adopted.


## 3.1. Population Initialization

In MODESA, initial population is generated in the same way as that of MODEA [23]. Because of the lack of a prior information about the solution, its corresponding opposite estimate is considered for getting a better solution. Here, the concept of opposite solution is given. Suppose that the current solution is x
(
x
1
,
x
2
,
…
,
x
i
,
…
,
x
n
) and the lower and upper bounds of x
i are v
i
,
min
⁡ and v
i
,
max
⁡, respectively. Then the opposite solution x
′
(
x
1
′
,
x
2
′
,
,
…
,
x
i
′
,
…
,
x
n
′
,
) of x is calculated as x
i
′
=
v
i
,
min
⁡
+
v
i
,
max
⁡
-
x
i [23].The advantage of doing this is that we can choose the better one between the solution and the opposite one as the initial solution. The procedure of initialization in the paper is given as follows. A population withN solutions is generated using uniform random distribution and a population containing their opposite solutions is constructed too. Next, a combined population with the population and the opposite population is constructed and the best N solutions of the combined population are selected by using the nondominated sorting in NSGA-II [2] and the diversity maintenance mechanism based on vicinity distance [30].


## 3.2. Differential Evolution (DE)

Differential Evolution (DE) presented by Storn and Price [8, 9] is an easy and efficient evolutionary algorithm. DE maintains a population of N solutions and contains three main operators: mutation, crossover, and selection.The mutation operation is defined as follows:(4)
z
=
x
r
1
+
f
*
(
x
r
2
-
x
r
3
)
,
where x
r
1, x
r
2, and x
r
3 are randomly selected from the current population and r
1
≠
r
2
≠
r
3
≠
i. f is a control parameter which is set in (
0,2
] in Storn and Price [8, 9]. For the sake of convenient description, the solution which is selected to plus a perturbed value is called as the base vector. x
r
1 (the base vector) is perturbed by adding to it the product of the control parameter and the difference of x
r
2 between x
r
3. z
(
z
1
,
z
2
,
…
,
z
n
) is the generated perturbed solution. It is important to note that if z
i is out of bounds, it must be repaired. As mentioned above, the lower and upper bounds of variable z
i are v
i
,
min
⁡ and v
i
,
max
⁡, respectively. In this study, if z
i is less than v
i
,
min
⁡, z
i is set to v
i
,
min
⁡; if z
i exceeds v
i
,
max
⁡, z
i is set to v
i
,
max
⁡. Ali et al. [23] have proved that using the one having the best fitness value in these three randomly selected individuals as the base vector has two advantages. The first advantage is that it is neither merely random nor merely greedy. The other one is that it can offer a localized effect to help in finding more different regions of the search space around the underlying solutions [23]. Therefore, in MODESA, after choosing the three individuals, we select the best one of them to be the base vector.After the perturbed solution is generated, a trial solutionw is produced by the crossover operation which can be defined as follows. For each variable j,
(5)
w
j
=
{
z
j
,
rand
j
≤
cr
∨
j
=
j
j
x
j
,
otherwise
,
where j
j
∈
{
1
,
…
,
n
} is a random number and cr
∈
[
0,1
] is the crossover rate. In general DE, the selection operator is utilized to decide whether the trial solution w replaces the target solution x
.The selection operation can be defined according to the following:
(6)
x
=
{
w
,
w
≺
x
x
,
otherwise
.In the proposed MODESA, the target individualx is not simply replaced by the trial solution w. A simulated annealing procedure introduced in Section 3.3 is employed to determine the acceptance of the trial solution.


## 3.3. Simulated Annealing Procedure

After the mutation phase, the crossover operator is performed, and the algorithm enters the selection phase. In MODESA, we incorporate the principle of simulated annealing into the selection phase. Before introducing the proposed selection phase, we first give a concept of domination in simulated annealing.


### 3.3.1. The Concept of Domination

To calculate the acceptance probability of a new solution, the amount of domination used in AMOSA [28] is adopted here. Given two solutions x and y, the amount of domination is calculated by the following equation:
(7)
Δ
dom
⁡
a
,
b
=
∏
i
=
1
,
f
i
(
x
)
≠
f
i
(
y
)
m
(
|
f
i
(
x
)
-
f
i
(
y
)
|
R
i
)
,
where m is the number of objectives and R
i is the range of the ith objective. When R
i is unknown, the solutions in the current population are used for calculating it.


### 3.3.2. Simulated Annealing Used in the Selection Operator of the Proposed MODESA

As (6), if the trial solution dominates the current solution, the general DE simply replaces the current solution by the trial solution. However, doing this may make the algorithm over converge and easy to trap into local optimal. As is known, there are three possibilities when comparing two solutions, the current solution x and the new solution y, as follows:(1)
xdominates y;
(2)
x is dominated by y;
(3)
x and y are incomparable.Ifx dominates y, it is obvious that y is not better than x, so that we can simply replace y by x. If x is dominated by y, we accept the new solution y distinctly. If x and y are incomparable, the solution y may have the potential to guide the solutions towards an unreached area of the search space. Therefore, in the DE process, we use the simulated annealing to accept y under a certain probability. Meanwhile, we need to control the acceptance probability, and the reason is that if we do not control the acceptance probability, all these potential solutions which are not absolutely good will enter to the next generation which will confuse the algorithm and make it unable to have a right direction towards the Pareto-optimal front. That is to say, we must decrease the acceptance probability. In simulated annealing, with the decreasing of the temperature, the acceptance probability is decreased too, which will avoid the algorithm producing too much potential solutions in the end. It is worth mentioning that if we allow some potential solutions existing during the procedure, the next nondominated sorting procedure may again remove these potential solutions, which will make the purpose of allowing some potential solutions to search more areas of the searching space fail. Therefore, in the proposed MODESA, each individual is assigned a prior life cycle. The number of prior life cycles means how many times the individual can be selected preferentially in the nondominated sorting. For the sake of convenience, the life cycle value of solution x
i is denoted as l
i
f
e
i. Once the individual is selected preferentially in the nondominated sorting, the prior life cycle of the individual is minus one when it is greater than zero. Note that, in the population initialization, the prior life cycle of each individual is set to zero.In this paper, given a temperaturetemp, the new solution is selected with a probability prob as follows:
(8)
prob
=
exp
⁡
⁡
(
-
Δ
dom
⁡
avg
temp
)
,
(9)
Δ
dom
⁡
avg
=
∑
i
=
1
k
Δ
dom
⁡
i
,
new-solution
k
,
where k is the number of solutions that dominates the new solution. It is worth mentioning that Δ
dom
⁡
avg represents the average amount of domination of the new solution and the k solutions that dominate it. Suppose that Δ
dom
⁡
avg is constant, the acceptance probability decreases when the temperature decreases. The simulated annealing used in the selection operator of MODESA is as shown in Algorithm 1.Algorithm 1:The proposed simulate annealing algorithm.
input: x
i, w
i, t
e
m
p
e
r
a
t
u
r
e, T
min
⁡, T
max
⁡, P
t, P
temp, m
a
x
_
l
i
f
e
_
c
i
c
l
e
If (t
e
m
p
e
r
a
t
u
r
e
>
T
min
⁡) {
Compute the value of Δ
dom
⁡
avg by (9).
Calculate the acceptance probability prob:
prob
=
ex
p
⁡
(
-
Δ
dom
⁡
avg
/
t
e
m
p
e
r
a
t
u
r
e
)
With a probability prob Do
Put the current individual x
i to the temporary population: P
temp
=
P
temp
∪
x
i.
Replace the current individual x
i by the trial individual w
i: x
i
=
w
i.
Assign the life value m
a
x
_
l
i
f
e
_
c
y
c
l
e to the new individual x
i:  l
i
f
e
i
=
m
a
x
_
l
i
f
e
_
c
y
c
l
e,
wherel
i
f
e
i represents the life cycle value of solution x
i.
otherwise
Put the current individual x
i to the temporary population: P
temp
=
P
temp
∪
x
i.
}// End If


## 3.3.1. The Concept of Domination

To calculate the acceptance probability of a new solution, the amount of domination used in AMOSA [28] is adopted here. Given two solutions x and y, the amount of domination is calculated by the following equation:
(7)
Δ
dom
⁡
a
,
b
=
∏
i
=
1
,
f
i
(
x
)
≠
f
i
(
y
)
m
(
|
f
i
(
x
)
-
f
i
(
y
)
|
R
i
)
,
where m is the number of objectives and R
i is the range of the ith objective. When R
i is unknown, the solutions in the current population are used for calculating it.


## 3.3.2. Simulated Annealing Used in the Selection Operator of the Proposed MODESA

As (6), if the trial solution dominates the current solution, the general DE simply replaces the current solution by the trial solution. However, doing this may make the algorithm over converge and easy to trap into local optimal. As is known, there are three possibilities when comparing two solutions, the current solution x and the new solution y, as follows:(1)
xdominates y;
(2)
x is dominated by y;
(3)
x and y are incomparable.Ifx dominates y, it is obvious that y is not better than x, so that we can simply replace y by x. If x is dominated by y, we accept the new solution y distinctly. If x and y are incomparable, the solution y may have the potential to guide the solutions towards an unreached area of the search space. Therefore, in the DE process, we use the simulated annealing to accept y under a certain probability. Meanwhile, we need to control the acceptance probability, and the reason is that if we do not control the acceptance probability, all these potential solutions which are not absolutely good will enter to the next generation which will confuse the algorithm and make it unable to have a right direction towards the Pareto-optimal front. That is to say, we must decrease the acceptance probability. In simulated annealing, with the decreasing of the temperature, the acceptance probability is decreased too, which will avoid the algorithm producing too much potential solutions in the end. It is worth mentioning that if we allow some potential solutions existing during the procedure, the next nondominated sorting procedure may again remove these potential solutions, which will make the purpose of allowing some potential solutions to search more areas of the searching space fail. Therefore, in the proposed MODESA, each individual is assigned a prior life cycle. The number of prior life cycles means how many times the individual can be selected preferentially in the nondominated sorting. For the sake of convenience, the life cycle value of solution x
i is denoted as l
i
f
e
i. Once the individual is selected preferentially in the nondominated sorting, the prior life cycle of the individual is minus one when it is greater than zero. Note that, in the population initialization, the prior life cycle of each individual is set to zero.In this paper, given a temperaturetemp, the new solution is selected with a probability prob as follows:
(8)
prob
=
exp
⁡
⁡
(
-
Δ
dom
⁡
avg
temp
)
,
(9)
Δ
dom
⁡
avg
=
∑
i
=
1
k
Δ
dom
⁡
i
,
new-solution
k
,
where k is the number of solutions that dominates the new solution. It is worth mentioning that Δ
dom
⁡
avg represents the average amount of domination of the new solution and the k solutions that dominate it. Suppose that Δ
dom
⁡
avg is constant, the acceptance probability decreases when the temperature decreases. The simulated annealing used in the selection operator of MODESA is as shown in Algorithm 1.Algorithm 1:The proposed simulate annealing algorithm.
input: x
i, w
i, t
e
m
p
e
r
a
t
u
r
e, T
min
⁡, T
max
⁡, P
t, P
temp, m
a
x
_
l
i
f
e
_
c
i
c
l
e
If (t
e
m
p
e
r
a
t
u
r
e
>
T
min
⁡) {
Compute the value of Δ
dom
⁡
avg by (9).
Calculate the acceptance probability prob:
prob
=
ex
p
⁡
(
-
Δ
dom
⁡
avg
/
t
e
m
p
e
r
a
t
u
r
e
)
With a probability prob Do
Put the current individual x
i to the temporary population: P
temp
=
P
temp
∪
x
i.
Replace the current individual x
i by the trial individual w
i: x
i
=
w
i.
Assign the life value m
a
x
_
l
i
f
e
_
c
y
c
l
e to the new individual x
i:  l
i
f
e
i
=
m
a
x
_
l
i
f
e
_
c
y
c
l
e,
wherel
i
f
e
i represents the life cycle value of solution x
i.
otherwise
Put the current individual x
i to the temporary population: P
temp
=
P
temp
∪
x
i.
}// End If


## 3.4. Diversity Maintenance Mechanism

The diversity maintenance mechanism adopted in the paper is presented by Kukkonen and Deb [30] for pruning of nondominated solutions. Instead of using the crowding distance in Deb et al. [2], the method proposed a new crowding estimation called vicinity distance using nearest neighbors. In this paper, we call this method as the vicinity distance based pruning method for convenience. The basic process of the method is to remove the most crowded solutions in the nondominated solution one by one. Moreover, after each removal, the vicinity distance values of the remaining solutions must be updated again. The vicinity distance is used to get an estimate of the density of solutions surrounding a particular solution, it calculates the product of k nearest neighbors of the particular solution. The solutions with the smallest value of vicinity distance are considered to be the most crowded ones. In [30], the value of k is assigned to the number of objectives; that is, k
=
m. The reader can refer to [30] to get more information about this method.


## 3.5. The Framework of the Proposed MODESA Algorithm

### 3.5.1. The Framework of the Proposed MODESA Algorithm

As mentioned before, the proposed algorithm first uses the population initialization phase to generate a population withN solutions. After the nondominated sorting, each solution has a ranking value. Next, the main loop begins. We introduce the tth generation of MODESA in the following.At thetth generation, an empty temporary population is created first. The temporary population is going to store the new good solutions generated by the DE in the loop. Then for each individual in the current population, the DE operators are used to generate a trial individual. Next, there are three possibilities when comparing the current individual and the trail individual. If the trial individual dominates the current individual, we simply replace the current solution by the trial solution. If the trial individual and the current solution are incomparable, the proposed simulated annealing process is adopted to determine whether the trial individual is accepted. If the trial individual is dominated by the current solution, we also put the trial individual into the temporary population as MODEA does. After all the individuals are handled, the current population and the temporary population are combined as a union population. Note that, if the current generation is the last generation, the members that are dominated by any other individuals must be removed in the union population. The reason of doing this is that, because we allow some potential solutions which are not absolutely good and have a positive life cycle to exist in the population; therefore, at the end of the generation, these solutions must be removed or they will influence the effectiveness of the results. Finally, the nondominated sorting and the vicinity distance based pruning approach are used for getting the best N solutions.In this paper, the parameters that need to be set a priori are listed in Table1 and the Pseudocode of MODESA is shown in Pseudocode 1.Table 1
The parameters used in the proposed algorithm.
Parameter's name
Description
Value
N
The population size
100
m
Number of objectives
As the test problem
n
Number of variables in a solution
As the test problem
T
Maximum number of generations
250
cr
The crossover rate used in DE
0.3
f
The control parameter used in DE
0.5
T
max
⁡
Maximum temperature in simulated annealing
100
T
min
⁡
Minimal temperature in simulated annealing
1
e
-
7
γ
The cooling rate in simulated annealing
0.6
m
a
x
_
l
i
f
e
_
c
i
c
l
e
Maximum prior life cycle value of the individual
1Pseudocode 1:Pseudocode of MODESA.
Step  1. Initialization
Initialize the values of the parameters in Table1.
Generate the initial populationP
0 using uniform distribution.
GenerateN opposite solutions denoted as P
obl by using OBL operation on P
0.
SelectN best individuals as the current population by using non-dominated sorting and the
vicinity distance based pruning method onP
0
∪
P
obl.
For (i
=
0; i
<
N; ++i) {
l
i
f
e
i
=
0.
}
/
/ End For
t
=
0
Step  2. Main loop
While (t
≤
T) {
t
e
m
p
e
r
a
t
u
r
e
=
T
max
⁡.
Construct a temporary population which is denoted asP
temp and set P
temp
=
∅.
For (i
=
0; i
<
N; ++i) {
Randomly select three distinct individuals:x
1, x
2 and x
3, which are different
from the current individualx
i from the current population.
Select the best individualx
best from the three distinct individuals as the base vector
for the next mutation operation.
Produce a trial individualw
i:
generate a random numberj
j
∈
{
1
,
…
,
n
}
For each variablej of w
i
{
With a probability cr or j is equal to j
j
do
w
j
,
i
=
x
j
,
1
+
f
*
(
x
j
,
2
-
x
j
,
3
).
otherwise
w
j
,
i
=
x
j
,
best.
If (w
j
,
i
>
v
j
,
max
⁡) w
j
,
i
=
v
j
,
max
⁡.
If (w
j
,
i
<
v
j
,
min
⁡) w
j
,
i
=
v
j
,
min
⁡.
}// End For
Check the domination status of the trial individual w
i and the current individual x
i.
If (w
i dominates x
i) {
Replace the current individual x
i by the trial individual w
i: x
i
=
w
i.
}
Else if (w
i and x
i are incomparable) {
Use the proposed simulated annealing to allocate the two individuals: w
i and x
i.
SimulateAnnealing (x
i, w
i, t
e
m
p
e
r
a
t
u
r
e, T
min
⁡, T
max
⁡, P
t, P
temp,  m
a
x
_
l
i
f
e
_
c
i
c
l
e).
t
e
m
p
e
r
a
t
u
r
e
=
γ
*
t
e
m
p
e
r
a
t
u
r
e.
}
Else if (w
i is dominated by x
i) {
Put the trial individual w
i to the temporary population.
}// End If
}// End For
Union the two population:P
t
∪
P
adv.
If (t
=
=
T) {
Use the non-dominated sorting to the union population to remove the individuals that
are dominated by any other solutions.
}// End If
SelectN individuals as the next population from the union population by using the
non-dominated sorting.
For (i
=
0; i
<
N; ++i) {
If (l
i
f
e
i
>
0)
l
i
f
e
i
-
-;
}// End For
Apply the vicinity distance based pruning method.
}// End While


### 3.5.2. The Computational Complexity Analysis

The proposed MODESA approach is a simple algorithm to be implemented. In the following, the overall computational complexity of MODESA is analyzed. Basic operations and their worst case complexities are given as follows.(1)
The initialization phase is selectingN solutions out of 2
N solutions by using nondominated sorting and the vicinity distance base pruning method: O
(
m
×
4
N
2
)
+
O
(
m
2
×
2
N
log
⁡
⁡
(
2
N
)
).
(2)
The procedure of mutation isO
(
m
×
2
N
).
(3)
The procedure of checking the domination status of the trial solution and the current solution isO
(
m
×
N
), and the procedure of simulated annealing is O
(
(
m
+
m
×
N
)
×
N
).
(4)
The procedure of selectingN solutions out of 2
N solutions by using nondominated sorting and the vicinity distance base pruning method: O
(
m
×
4
N
2
)
+
O
(
m
2
×
2
N
log
⁡
⁡
(
2
N
)
).Here,m and N represent the number of objectives and the population size, respectively. We can see that the overall worst case complexity of the MODESA is O
(
m
2
N
log
⁡
⁡
N
+
m
N
2
).


### 3.5.3. The Difference between MODESA and MODEA

Considering that the proposed algorithm MODESA is an improvement of MODEA, it is necessary to provide the difference between MODESA and MODEA. There are two main differences between MODESA and MODEA as follows.(1)
When comparing the new solution and the current solution, there exist two processing ways according to two different results of comparison. If the new solution dominates the current solution, the MODEA replaces the current solution with the new solution and puts the current solution into the temporary population; otherwise, the new solution is taken into the temporary population. After each solution in the population is handled, the two populations are combined. It indicates that the new solutions are all retained. However, MODESA works in three kinds of ways according to different results. If the new solution dominates the current solution, the current solution is replaced with the new solution and the current solution is removed directly, which means that if the current solution is worse than the new solution, it would not be retained. If the new solution and the current solution are incomparable, the proposed simulated annealing algorithm is adopted to determine the acceptance of the current solution; moreover, the usage of a prior life cycle enables some potential solutions preferentially to be selected in the next generation. Otherwise, the new solution is added to the temporary population and the current solution is retained in the current population.
(2)
The diversity maintenance mechanism is different. In MODEA, the diversity preserving strategy based on crowding distance in NSGA-II is utilized, while in MODESA, the vicinity distance based pruning method in [30] is adopted.


## 3.5.1. The Framework of the Proposed MODESA Algorithm

As mentioned before, the proposed algorithm first uses the population initialization phase to generate a population withN solutions. After the nondominated sorting, each solution has a ranking value. Next, the main loop begins. We introduce the tth generation of MODESA in the following.At thetth generation, an empty temporary population is created first. The temporary population is going to store the new good solutions generated by the DE in the loop. Then for each individual in the current population, the DE operators are used to generate a trial individual. Next, there are three possibilities when comparing the current individual and the trail individual. If the trial individual dominates the current individual, we simply replace the current solution by the trial solution. If the trial individual and the current solution are incomparable, the proposed simulated annealing process is adopted to determine whether the trial individual is accepted. If the trial individual is dominated by the current solution, we also put the trial individual into the temporary population as MODEA does. After all the individuals are handled, the current population and the temporary population are combined as a union population. Note that, if the current generation is the last generation, the members that are dominated by any other individuals must be removed in the union population. The reason of doing this is that, because we allow some potential solutions which are not absolutely good and have a positive life cycle to exist in the population; therefore, at the end of the generation, these solutions must be removed or they will influence the effectiveness of the results. Finally, the nondominated sorting and the vicinity distance based pruning approach are used for getting the best N solutions.In this paper, the parameters that need to be set a priori are listed in Table1 and the Pseudocode of MODESA is shown in Pseudocode 1.Table 1
The parameters used in the proposed algorithm.
Parameter's name
Description
Value
N
The population size
100
m
Number of objectives
As the test problem
n
Number of variables in a solution
As the test problem
T
Maximum number of generations
250
cr
The crossover rate used in DE
0.3
f
The control parameter used in DE
0.5
T
max
⁡
Maximum temperature in simulated annealing
100
T
min
⁡
Minimal temperature in simulated annealing
1
e
-
7
γ
The cooling rate in simulated annealing
0.6
m
a
x
_
l
i
f
e
_
c
i
c
l
e
Maximum prior life cycle value of the individual
1Pseudocode 1:Pseudocode of MODESA.
Step  1. Initialization
Initialize the values of the parameters in Table1.
Generate the initial populationP
0 using uniform distribution.
GenerateN opposite solutions denoted as P
obl by using OBL operation on P
0.
SelectN best individuals as the current population by using non-dominated sorting and the
vicinity distance based pruning method onP
0
∪
P
obl.
For (i
=
0; i
<
N; ++i) {
l
i
f
e
i
=
0.
}
/
/ End For
t
=
0
Step  2. Main loop
While (t
≤
T) {
t
e
m
p
e
r
a
t
u
r
e
=
T
max
⁡.
Construct a temporary population which is denoted asP
temp and set P
temp
=
∅.
For (i
=
0; i
<
N; ++i) {
Randomly select three distinct individuals:x
1, x
2 and x
3, which are different
from the current individualx
i from the current population.
Select the best individualx
best from the three distinct individuals as the base vector
for the next mutation operation.
Produce a trial individualw
i:
generate a random numberj
j
∈
{
1
,
…
,
n
}
For each variablej of w
i
{
With a probability cr or j is equal to j
j
do
w
j
,
i
=
x
j
,
1
+
f
*
(
x
j
,
2
-
x
j
,
3
).
otherwise
w
j
,
i
=
x
j
,
best.
If (w
j
,
i
>
v
j
,
max
⁡) w
j
,
i
=
v
j
,
max
⁡.
If (w
j
,
i
<
v
j
,
min
⁡) w
j
,
i
=
v
j
,
min
⁡.
}// End For
Check the domination status of the trial individual w
i and the current individual x
i.
If (w
i dominates x
i) {
Replace the current individual x
i by the trial individual w
i: x
i
=
w
i.
}
Else if (w
i and x
i are incomparable) {
Use the proposed simulated annealing to allocate the two individuals: w
i and x
i.
SimulateAnnealing (x
i, w
i, t
e
m
p
e
r
a
t
u
r
e, T
min
⁡, T
max
⁡, P
t, P
temp,  m
a
x
_
l
i
f
e
_
c
i
c
l
e).
t
e
m
p
e
r
a
t
u
r
e
=
γ
*
t
e
m
p
e
r
a
t
u
r
e.
}
Else if (w
i is dominated by x
i) {
Put the trial individual w
i to the temporary population.
}// End If
}// End For
Union the two population:P
t
∪
P
adv.
If (t
=
=
T) {
Use the non-dominated sorting to the union population to remove the individuals that
are dominated by any other solutions.
}// End If
SelectN individuals as the next population from the union population by using the
non-dominated sorting.
For (i
=
0; i
<
N; ++i) {
If (l
i
f
e
i
>
0)
l
i
f
e
i
-
-;
}// End For
Apply the vicinity distance based pruning method.
}// End While


## 3.5.2. The Computational Complexity Analysis

The proposed MODESA approach is a simple algorithm to be implemented. In the following, the overall computational complexity of MODESA is analyzed. Basic operations and their worst case complexities are given as follows.(1)
The initialization phase is selectingN solutions out of 2
N solutions by using nondominated sorting and the vicinity distance base pruning method: O
(
m
×
4
N
2
)
+
O
(
m
2
×
2
N
log
⁡
⁡
(
2
N
)
).
(2)
The procedure of mutation isO
(
m
×
2
N
).
(3)
The procedure of checking the domination status of the trial solution and the current solution isO
(
m
×
N
), and the procedure of simulated annealing is O
(
(
m
+
m
×
N
)
×
N
).
(4)
The procedure of selectingN solutions out of 2
N solutions by using nondominated sorting and the vicinity distance base pruning method: O
(
m
×
4
N
2
)
+
O
(
m
2
×
2
N
log
⁡
⁡
(
2
N
)
).Here,m and N represent the number of objectives and the population size, respectively. We can see that the overall worst case complexity of the MODESA is O
(
m
2
N
log
⁡
⁡
N
+
m
N
2
).


## 3.5.3. The Difference between MODESA and MODEA

Considering that the proposed algorithm MODESA is an improvement of MODEA, it is necessary to provide the difference between MODESA and MODEA. There are two main differences between MODESA and MODEA as follows.(1)
When comparing the new solution and the current solution, there exist two processing ways according to two different results of comparison. If the new solution dominates the current solution, the MODEA replaces the current solution with the new solution and puts the current solution into the temporary population; otherwise, the new solution is taken into the temporary population. After each solution in the population is handled, the two populations are combined. It indicates that the new solutions are all retained. However, MODESA works in three kinds of ways according to different results. If the new solution dominates the current solution, the current solution is replaced with the new solution and the current solution is removed directly, which means that if the current solution is worse than the new solution, it would not be retained. If the new solution and the current solution are incomparable, the proposed simulated annealing algorithm is adopted to determine the acceptance of the current solution; moreover, the usage of a prior life cycle enables some potential solutions preferentially to be selected in the next generation. Otherwise, the new solution is added to the temporary population and the current solution is retained in the current population.
(2)
The diversity maintenance mechanism is different. In MODEA, the diversity preserving strategy based on crowding distance in NSGA-II is utilized, while in MODESA, the vicinity distance based pruning method in [30] is adopted.


## 4. Results and Discussion

### 4.1. Test Problems

The performance of the proposed algorithm is tested on a set of five biobjective and two triobjective optimization problems. The first five problems are ZDT1, ZDT2, ZDT3, ZDT4, and ZDT6, which were described in [31]. The two triobjective optimization problems: DTLZ1 and DTLZ2 are taken from Deb [32]. None of these problems has any constraint. All objective functions are to be minimized. The number of variables, their bounds, the Pareto-optimal solutions, and the characteristics of the Pareto-optimal front of these test problems are also shown in Table 2 [2]. For the test problems DTLZ1 and DTLZ2, the number of objectives is set to 3. Therefore, the number of variables for DTLZ1 and DTLZ2 is equal to 7 and 12, respectively.Table 2
Test problems used in the paper.
Problem
n
Variable bounds
Objective functions
Optimal solutions
Comments
ZDT1
30
[
0,1
]
f
1
=
x
1
x
1
∈
[
0,1
]
Convex
f
2
=
g
(
x
)
[
1
-
x
1
g
(
x
)
]
x
=
0
i
=
2
,
…
,
n
g
(
x
)
=
1
+
9
(
∑
i
=
2
n
x
i
)
n
-
1
ZDT2
30
[
0,1
]
f
1
=
x
1
x
1
∈
[
0,1
]
Nonconvex
f
2
=
g
(
x
)
[
1
-
x
1
g
(
x
)
]
2
x
=
0
i
=
2
,
…
,
n
g
(
x
)
=
1
+
9
(
∑
i
=
2
n
x
i
)
n
-
1
ZDT3
30
[
0,1
]
f
1
=
x
1
x
1
∈
[
0,1
]
Convex disconnected
f
2
=
g
(
x
)
[
1
-
x
1
g
(
x
)
-
x
1
g
(
x
)
sin
⁡
⁡
(
10
π
x
1
)
]
x
=
0
g
(
x
)
=
1
+
9
∑
i
=
2
n
x
i
n
-
1
i
=
2
,
…
,
n
ZDT4
10
x
1
∈
[
0,1
]
f
1
=
x
1
x
1
∈
[
0,1
]
Nonconvex
x
i
∈
[
-
5,5
]
f
2
=
g
(
x
)
[
1
-
x
1
g
(
x
)
]
x
=
0
g
(
x
)
=
1
+
10
(
n
-
1
)
+
∑
i
=
2
n
[
x
i
2
-
10
cos
⁡
⁡
(
4
π
x
i
)
]
i
=
2
,
…
,
n
ZDT6
10
[
0,1
]
f
1
=
1
-
exp
⁡
(
-
4
x
i
)
sin
⁡
6
(
6
π
x
1
)
x
1
∈
[
0,1
]
Nonconvex
f
2
=
g
(
x
)
[
1
-
(
f
1
g
(
x
)
)
2
]
x
=
0
g
(
x
)
=
1
+
9
(
∑
i
=
2
n
x
i
n
-
1
)
0.25
i
=
2
,
…
,
n
DTLZ1
m
+
4
[
0,1
]
f
1
=
1
2
x
1
x
2
…
x
m
-
1
(
1
+
g
(
x
m
)
)
x
i
=
0.5
i
=
1,2
,
…
,
n
l
l
l
l
f
2
=
1
2
x
1
x
2
…
x
m
-
1
(
1
-
x
m
-
1
)
(
1
+
g
(
x
m
)
)
⋮
f
m
-
1
=
1
2
x
1
(
1
-
x
2
)
(
1
+
g
(
x
m
)
)
f
m
=
1
2
(
1
-
x
1
)
(
1
+
g
(
x
m
)
)
g
(
x
m
)
=
100
[
|
x
m
|
+
∑
x
i
ϵ
x
m
(
x
i
-
0.5
)
2
-
cos
⁡
(
20
π
(
x
i
-
0.5
)
)
]
DTLZ2
m
+
9
[
0,1
]
l
l
l
l
l
f
1
=
(
1
+
g
(
x
m
)
)
cos
⁡
⁡
(
x
1
π
2
)
cos
⁡
⁡
(
x
2
π
2
)
⋯
cos
⁡
⁡
(
x
m
-
2
π
2
)
cos
⁡
⁡
(
x
m
-
1
π
2
)
x
i
=
0.5
i
=
1,2
,
…
,
n
l
l
l
l
l
f
2
=
(
1
+
g
(
x
m
)
)
cos
⁡
⁡
(
x
1
π
2
)
cos
⁡
⁡
(
x
2
π
2
)
⋯
cos
⁡
⁡
(
x
m
-
2
π
2
)
sin
⁡
⁡
(
x
m
-
1
π
2
)
l
l
l
l
l
f
3
=
(
1
+
g
(
x
m
)
)
cos
⁡
⁡
(
x
1
π
2
)
cos
⁡
⁡
(
x
2
π
2
)
⋯
sin
⁡
⁡
(
x
m
-
2
π
2
)
⋮
l
l
l
l
l
f
m
=
(
1
+
g
(
x
m
)
)
sin
⁡
⁡
(
x
1
π
2
)
l
l
l
l
l
g
(
x
m
)
=
∑
x
i
ϵ
x
m
(
x
i
-
0.5
)
2


### 4.2. Performance Measures

In this section, three performance measures will be introduced to evaluate the effectiveness of the proposed algorithm MODESA with the other algorithms. Here, letP
* and P denote the Pareto-optimal front and the Pareto front which are obtained by the algorithm, respectively. In this study, three performance measures are used for evaluating the convergence, the diversity, and both of them, separately. It is worth mentioning that all the three performance measures can be used only when the Pareto-optimal front of the particular test problem is known. In this paper, 500 evenly distributed solutions in PF are selected as P
* for each biobjective problem, and 990 solutions for triobjective problem, which is the same as that in [33]. The reader can download these solutions from the website http://dces.essex.ac.uk/staff/qzhang/index.html.


#### 4.2.1. Inverted Generational Distance (IGD)

The concept of inverted generational distance (IGD) is proposed by Coello and Cortés in [34]. The metric can be formally described as follows:
(10)
IGD
=
∑
x
∈
P
*
m
i
n
_
e
d
(
x
,
P
)
|
P
*
|
,
where m
i
n
_
e
d
(
x
,
P
) is the minimum Euclidean distance between the solution x and the solutions in P. The IGD-metric value can evaluate both the convergence and diversity of P to a certain degree. Algorithms with smaller value of IGD are satisfying.


#### 4.2.2. Generational Distance (GD)

van Veldhuizen and Lamont [35] proposed the concept of generational distance (GD). The metric calculates the distance between P and P
*. It can be defined as follows:
(11)
GD
=
∑
x
∈
P
m
i
n
_
e
d
(
x
,
p
*
)
2
|
P
|
,
where m
i
n
_
e
d
(
x
,
P
*
) is the minimum Euclidean distance of solution x and the solutions in P
*. GD evaluates the closeness of the obtained Pareto front and the Pareto-optimal front. Algorithms with smaller value of GD are desirable.


#### 4.2.3. Spread(
Δ
)

The spread metricΔ was firstly proposed by Deb et al. [2] to measure the spread effect of the nondominated solutions obtained by the algorithm. However, it is only suitable for biobjective problems. Wang et al. [19] further extended the metric Δ and made it suitable to measure problems that have more than two objectives. The improved Δ can be stated mathematically as follows:
(12)
Δ
=
∑
i
=
1
m
d
(
E
i
,
Ω
)
+
∑
Χ
∈
Ω
|
d
(
X
,
Ω
)
-
d
-
|
∑
i
=
1
m
d
(
E
i
,
Ω
)
+
(
|
Ω
|
-
m
)
d
-
d
(
X
,
Ω
)
=
min
⁡
Y
∈
Ω
,
Y
≠
X
⁡
∥
F
(
X
)
-
F
(
Y
)
∥
d
-
=
1
|
Ω
|
∑
X
∈
Ω
d
(
X
,
Ω
)
,
where Ω is a set of solutions, m is the number of objectives, and (E
1
,
…
,
E
m) are m extreme solutions in the set of true Pareto front (PF). The smaller the value of Δ, the better the distribution and diversity of the obtained nondominated solutions.


### 4.3. Experimental Setup

In this paper, the population size and the number of generations are set to 100 and 250, respectively. For the other parameters in NSGA-II and MODEA, they have remained the same with their original studies [3, 23]. For MODESA, there are four parameters: T
max
⁡, T
min
⁡, γ, and m
a
x
_
l
i
f
e
_
c
i
c
l
e that need to be set. Here we set T
max
⁡
=
100, T
min
⁡
=
1
e
-
7, γ
=
0.6, and m
a
x
⁡
_
l
i
f
e
_
c
y
c
l
e
=
1 for the study. Each algorithm runs 10 times independently for each problem. All the test problems have been executed on Microsoft Window XP Intel core 2 Duo CPU E8400 3.00 GHz, with 2 GB RAM. All methods have been implemented in C++.


### 4.4. Discussion of the Results

#### 4.4.1. Results of Biobjective Optimization Problems

Tables3 and 4 show the mean and standard deviation IGD-metric and Δ-metric values of the obtained nondominated solutions in 10 runs on the group of ZDTs. It is obvious from Tables 3 and 4 that the IGD-metric values of MODESA achieve the best value in the three algorithms except for the test problem ZDT4.Table 3
The mean and standard deviation IGD-metric values of nondominated solutions in 10 runs on ZDTs whereN
=
100 and T
=
250.
IGD
NSGA-II (mean ± std)
MODEA (mean ± std)
MODESA (mean ± std)
ZDT1
0.00564974
±
0.000254893
0.00436863
±
0.000243648
0.00403528
±
7.01666
e
-
05
ZDT2
0.00557228
±
0.000242512
0.00451929
±
0.000134666
0.00408871
±
4.65746
e
-
05
ZDT3
0.00592808
±
0.000205118
0.00521698
±
0.000159453
0.00494596
±
0.000126657
ZDT4
0.0069977
±
0.00102386
0.0301279
±
0.0752722
0.0171826
±
0.0362604
ZDT6
0.0137079
±
0.00170579
0.00304881
±
0.000125906
0.0022843
±
0.00014073Table 4
The mean and standard deviationΔ-metric values of nondominated solutions in 10 runs on ZDTs where N
=
100 and T
=
250.
Spread
NSGA-II (mean ± std)
MODEA (mean ± std)
MODESA (mean ± std)
ZDT1
0.372425
±
0.0513842
0.333973
±
0.0772084
0.297395
±
0.0308609
ZDT2
0.385209
±
0.0291848
0.336548
±
0.0181084
0.296842
±
0.0282846
ZDT3
0.389909
±
0.0557332
0.355518
±
0.0370526
0.266801
±
0.0211668
ZDT4
0.399555
±
0.0786079
0.377719
±
0.0574188
0.352871
±
0.107095
ZDT6
0.407279
±
0.0276428
0.407335
±
0.0487992
0.368972
±
0.04
63564The test problems ZDT1 and ZDT2 may be the most simple problems in the series of ZDTs. The Pareto-optimal front of ZDT1 is convex while that of ZDT2 is nonconvex. It is clear from Tables3 and 4 that both the values of the IGD metric and the spread metric of MODESA are better than that of NSGA-II and MODEA. Figures 1 and 2 demonstrate that MODESA converges to the Pareto-optimal front and has a good spread over the entire front on the test problems ZDT1 and ZDT2. It can be concluded that the proposed MODESA can well deal with these two kinds of test problems.Figure 1
Plots of the nondominated solutions with the lowest IGD-metric values found by MODESA in 10 runs in the objective space on ZDT1 whereN
=
100 and T
=
250.Figure 2
Plots of the nondominated solutions with the lowest IGD-metric values found by MODESA in 10 runs in the objective space on ZDT2 whereN
=
100 and T
=
250.The third test problem ZDT3 is different from the above two test problems. The Pareto-optimal front of ZDT3 consists of five disjoint curves. From Figure3, it is obvious that the nondominated solutions obtained by the MODESA almost converge to the Pareto-optimal front and have a good distribution over the entire front. From Tables 3 and 4, MODESA obtains the best values of IGD-metric and spread-metric between the three algorithms.Figure 3
Plots of the nondominated solutions with the lowest IGD-metric values found by MODESA in 10 runs in the objective space on ZDT3 whereN
=
100 and T
=
250.The fourth test problem is ZDT4 which is perhaps the most difficult problem in the five problems. It has21
9 different local Pareto-optimal fronts that confuse the multiobjective optimization algorithm [2]. From Tables 3 and 4, NSGA-II obtains the best values on the IGD-metric and the spread-metric in these three algorithms. In fact, in the 10 independent runs, some results of MODEA and MODESA are better than that of NSGA-II; however, MODEA and MODESA get bad results in a few runs which affects the mean value of the 10 runs. Indeed, the best IGD-metric value of MODESA in the 10 runs is 0.0044915, which is better than NSGA-II. Figure 4 demonstrates the plots of the nondominated solutions with the lowest IGD-metric values found by MODESA in the 10 runs. It is clear from Figure 4 that MODESA can get quite good convergence and distribution over the Pareto-optimal front, which reveals that the proposed algorithm is a promising but unstable algorithm for solving the test problem ZDT4. We will leave this problem for future work.Figure 4
Plots of the nondominated solutions with the lowest IGD-metric values found by MODESA in 10 runs in the objective space on ZDT4 whereN
=
100 and T
=
250.The last test problem ZDT6 is such a problem with thin density and nonuniformed spread of solutions [23]. We can observe from Tables 3 and 4 that MODESA finds the best convergence and spread in comparison to the other two algorithms. To show the superiority of MODESA, plots of the nondominated solutions found by MODESA and MODEA on ZDT 6 are shown in Figures 5 and 6, respectively. It is clear from Figures 5 and 6 that MODESA gets a better spread than MODEA.Figure 5
Plots of the nondominated solutions with the lowest IGD-metric values found by MODESA in 10 runs in the objective space on ZDT6 whereN
=
100 and T
=
250.Figure 6
Plots of the nondominated solutions with the lowest IGD-metric values found by MODEA in 10 runs in the objective space on ZDT6 whereN
=
100 and T
=
250.


#### 4.4.2. Results of Triobjective Optimization Problems

In this section, to show the effectiveness of the proposed MODESA for solving problems that have more than two objectives, we choose two triobjective optimization problems DTLZ1 and DTLZ2 for comparison. The results of NSGA-II, MODEA, and MODESA are given in Table5. It is obvious from Table 5 that MODESA performs better than the other two algorithms in terms of both the IGD-metric and spread-metric. Figures 7, 8, 9, and 10 show the nondominated solutions obtained by MODEA and MODESA, and it is obvious that MODESA outperforms MODEA on both of the two triobjective optimization problems. In particular, it is clear that MODESA has a better distribution than MODEA, the reason is that the vicinity distance based pruning method performs better than the crowding distance of Deb in maintaining diversity and distribution of the nondominated solutions.Table 5
The mean and standard deviation IGD-metric andΔ-metric values of nondominated solutions in 10 runs on DTLZ1 and DTLZ2 where N
=
100 and T
=
250.
Metric
IGD
Spread
Problem
DTLZ1
DTLZ2
DTLZ1
DTLZ2
NSGA-II (mean± std)
0.025008
±
0.001199
0.077375
±
0.021572
0.60056
±
0.076489
0.530832
±
0.029126
MODEA (mean± std)
0.022851
±
0.000557
0.060254
±
0.001664
0.41962
±
0.019617
0.425707
±
0.033528
MODESA (mean± std)
0.020602
±
0.00035
0.055011
±
0.000799
0.26489
±
0.045448
0.260335
±
0.02856Figure 7
Plots of the nondominated solutions with the lowest IGD-metric values found by MODEA in 10 runs in the objective space on DTLZ1 whereN
=
100 and T
=
250.Figure 8
Plots of the nondominated solutions with the lowest IGD-metric values found by MODESA in 10 runs in the objective space on DTLZ1 whereN
=
100 and T
=
250.Figure 9
Plots of the nondominated solutions with the lowest IGD-metric values found by MODEA in 10 runs in the objective space on DTLZ2 whereN
=
100 and T
=
250.Figure 10
Plots of the nondominated solutions with the lowest IGD-metric values found by MODESA in 10 runs in the objective space on DTLZ2 whereN
=
100 and T
=
250.


### 4.5. Effect of the Proposed Simulated Annealing Procedure

As mentioned in Section3.5.3, one of the main differences of MODEA and MODESA is the employment of the proposed simulated annealing procedure. To show the effect of the proposed simulated annealing procedure, a version of MODESA: MODESA-CD (multiobjective differential evolution with crowding distance) is proposed for comparison. The MODESA-CD is analogous to MODESA except that it replaces the vicinity distance based pruning method in MODESA with the crowding distance based approach in NSGA-II. In other words, the difference between MODEA and MODESA-CD is the using of the proposed simulated annealing procedure. Table 6 shows the mean and standard deviation GD-metric values of nondominated solutions in 10 runs where N
=
100 and T
=
250. It is obvious from Table 6 that all the GD-metric values of MODESA-CD are better than that of MODEA. The reason of the bad performance of MODEA is that MODEA always accepts the current and the new solutions to the next generation and that makes it chaotic and unable to have a right direction towards the Pareto-optimal front. In MODESA-CD, we remove the bad solutions and only accept some potential solutions in an alterable probability.Table 6
The mean and standard deviation GD-metric values of nondominated solutions in 10 runs whereN
=
100 and T
=
250.
GD
MODEA (mean ± std)
MODESA-CD (mean ± std)
ZDT1
0.000124
±
8.236
e
-
06
0.000116
±
1.437
e
-
05
ZDT2
4.819
e
-
05
±
1.642
e
-
06
4.771
e
-
05
±
2.596
e
-
06
ZDT3
0.00035
±
6.128
e
-
06
0.000345
±
2.714
e
-
05
ZDT4
0.002812
±
0.008057
0.001365
±
0.003738
ZDT6
0.000429
±
1.78
e
-
05
0.000404
±
2.614
e
-
05
DTLZ1
0.000308
±
7.454
e
-
05
0.00029
±
2.129
e
-
05
DTLZ2
0.000677
±
0.000103
0.000633
±
3.871
e
-
05


## 4.1. Test Problems

The performance of the proposed algorithm is tested on a set of five biobjective and two triobjective optimization problems. The first five problems are ZDT1, ZDT2, ZDT3, ZDT4, and ZDT6, which were described in [31]. The two triobjective optimization problems: DTLZ1 and DTLZ2 are taken from Deb [32]. None of these problems has any constraint. All objective functions are to be minimized. The number of variables, their bounds, the Pareto-optimal solutions, and the characteristics of the Pareto-optimal front of these test problems are also shown in Table 2 [2]. For the test problems DTLZ1 and DTLZ2, the number of objectives is set to 3. Therefore, the number of variables for DTLZ1 and DTLZ2 is equal to 7 and 12, respectively.Table 2
Test problems used in the paper.
Problem
n
Variable bounds
Objective functions
Optimal solutions
Comments
ZDT1
30
[
0,1
]
f
1
=
x
1
x
1
∈
[
0,1
]
Convex
f
2
=
g
(
x
)
[
1
-
x
1
g
(
x
)
]
x
=
0
i
=
2
,
…
,
n
g
(
x
)
=
1
+
9
(
∑
i
=
2
n
x
i
)
n
-
1
ZDT2
30
[
0,1
]
f
1
=
x
1
x
1
∈
[
0,1
]
Nonconvex
f
2
=
g
(
x
)
[
1
-
x
1
g
(
x
)
]
2
x
=
0
i
=
2
,
…
,
n
g
(
x
)
=
1
+
9
(
∑
i
=
2
n
x
i
)
n
-
1
ZDT3
30
[
0,1
]
f
1
=
x
1
x
1
∈
[
0,1
]
Convex disconnected
f
2
=
g
(
x
)
[
1
-
x
1
g
(
x
)
-
x
1
g
(
x
)
sin
⁡
⁡
(
10
π
x
1
)
]
x
=
0
g
(
x
)
=
1
+
9
∑
i
=
2
n
x
i
n
-
1
i
=
2
,
…
,
n
ZDT4
10
x
1
∈
[
0,1
]
f
1
=
x
1
x
1
∈
[
0,1
]
Nonconvex
x
i
∈
[
-
5,5
]
f
2
=
g
(
x
)
[
1
-
x
1
g
(
x
)
]
x
=
0
g
(
x
)
=
1
+
10
(
n
-
1
)
+
∑
i
=
2
n
[
x
i
2
-
10
cos
⁡
⁡
(
4
π
x
i
)
]
i
=
2
,
…
,
n
ZDT6
10
[
0,1
]
f
1
=
1
-
exp
⁡
(
-
4
x
i
)
sin
⁡
6
(
6
π
x
1
)
x
1
∈
[
0,1
]
Nonconvex
f
2
=
g
(
x
)
[
1
-
(
f
1
g
(
x
)
)
2
]
x
=
0
g
(
x
)
=
1
+
9
(
∑
i
=
2
n
x
i
n
-
1
)
0.25
i
=
2
,
…
,
n
DTLZ1
m
+
4
[
0,1
]
f
1
=
1
2
x
1
x
2
…
x
m
-
1
(
1
+
g
(
x
m
)
)
x
i
=
0.5
i
=
1,2
,
…
,
n
l
l
l
l
f
2
=
1
2
x
1
x
2
…
x
m
-
1
(
1
-
x
m
-
1
)
(
1
+
g
(
x
m
)
)
⋮
f
m
-
1
=
1
2
x
1
(
1
-
x
2
)
(
1
+
g
(
x
m
)
)
f
m
=
1
2
(
1
-
x
1
)
(
1
+
g
(
x
m
)
)
g
(
x
m
)
=
100
[
|
x
m
|
+
∑
x
i
ϵ
x
m
(
x
i
-
0.5
)
2
-
cos
⁡
(
20
π
(
x
i
-
0.5
)
)
]
DTLZ2
m
+
9
[
0,1
]
l
l
l
l
l
f
1
=
(
1
+
g
(
x
m
)
)
cos
⁡
⁡
(
x
1
π
2
)
cos
⁡
⁡
(
x
2
π
2
)
⋯
cos
⁡
⁡
(
x
m
-
2
π
2
)
cos
⁡
⁡
(
x
m
-
1
π
2
)
x
i
=
0.5
i
=
1,2
,
…
,
n
l
l
l
l
l
f
2
=
(
1
+
g
(
x
m
)
)
cos
⁡
⁡
(
x
1
π
2
)
cos
⁡
⁡
(
x
2
π
2
)
⋯
cos
⁡
⁡
(
x
m
-
2
π
2
)
sin
⁡
⁡
(
x
m
-
1
π
2
)
l
l
l
l
l
f
3
=
(
1
+
g
(
x
m
)
)
cos
⁡
⁡
(
x
1
π
2
)
cos
⁡
⁡
(
x
2
π
2
)
⋯
sin
⁡
⁡
(
x
m
-
2
π
2
)
⋮
l
l
l
l
l
f
m
=
(
1
+
g
(
x
m
)
)
sin
⁡
⁡
(
x
1
π
2
)
l
l
l
l
l
g
(
x
m
)
=
∑
x
i
ϵ
x
m
(
x
i
-
0.5
)
2


## 4.2. Performance Measures

In this section, three performance measures will be introduced to evaluate the effectiveness of the proposed algorithm MODESA with the other algorithms. Here, letP
* and P denote the Pareto-optimal front and the Pareto front which are obtained by the algorithm, respectively. In this study, three performance measures are used for evaluating the convergence, the diversity, and both of them, separately. It is worth mentioning that all the three performance measures can be used only when the Pareto-optimal front of the particular test problem is known. In this paper, 500 evenly distributed solutions in PF are selected as P
* for each biobjective problem, and 990 solutions for triobjective problem, which is the same as that in [33]. The reader can download these solutions from the website http://dces.essex.ac.uk/staff/qzhang/index.html.


### 4.2.1. Inverted Generational Distance (IGD)

The concept of inverted generational distance (IGD) is proposed by Coello and Cortés in [34]. The metric can be formally described as follows:
(10)
IGD
=
∑
x
∈
P
*
m
i
n
_
e
d
(
x
,
P
)
|
P
*
|
,
where m
i
n
_
e
d
(
x
,
P
) is the minimum Euclidean distance between the solution x and the solutions in P. The IGD-metric value can evaluate both the convergence and diversity of P to a certain degree. Algorithms with smaller value of IGD are satisfying.


### 4.2.2. Generational Distance (GD)

van Veldhuizen and Lamont [35] proposed the concept of generational distance (GD). The metric calculates the distance between P and P
*. It can be defined as follows:
(11)
GD
=
∑
x
∈
P
m
i
n
_
e
d
(
x
,
p
*
)
2
|
P
|
,
where m
i
n
_
e
d
(
x
,
P
*
) is the minimum Euclidean distance of solution x and the solutions in P
*. GD evaluates the closeness of the obtained Pareto front and the Pareto-optimal front. Algorithms with smaller value of GD are desirable.


### 4.2.3. Spread(
Δ
)

The spread metricΔ was firstly proposed by Deb et al. [2] to measure the spread effect of the nondominated solutions obtained by the algorithm. However, it is only suitable for biobjective problems. Wang et al. [19] further extended the metric Δ and made it suitable to measure problems that have more than two objectives. The improved Δ can be stated mathematically as follows:
(12)
Δ
=
∑
i
=
1
m
d
(
E
i
,
Ω
)
+
∑
Χ
∈
Ω
|
d
(
X
,
Ω
)
-
d
-
|
∑
i
=
1
m
d
(
E
i
,
Ω
)
+
(
|
Ω
|
-
m
)
d
-
d
(
X
,
Ω
)
=
min
⁡
Y
∈
Ω
,
Y
≠
X
⁡
∥
F
(
X
)
-
F
(
Y
)
∥
d
-
=
1
|
Ω
|
∑
X
∈
Ω
d
(
X
,
Ω
)
,
where Ω is a set of solutions, m is the number of objectives, and (E
1
,
…
,
E
m) are m extreme solutions in the set of true Pareto front (PF). The smaller the value of Δ, the better the distribution and diversity of the obtained nondominated solutions.


## 4.2.1. Inverted Generational Distance (IGD)

The concept of inverted generational distance (IGD) is proposed by Coello and Cortés in [34]. The metric can be formally described as follows:
(10)
IGD
=
∑
x
∈
P
*
m
i
n
_
e
d
(
x
,
P
)
|
P
*
|
,
where m
i
n
_
e
d
(
x
,
P
) is the minimum Euclidean distance between the solution x and the solutions in P. The IGD-metric value can evaluate both the convergence and diversity of P to a certain degree. Algorithms with smaller value of IGD are satisfying.


## 4.2.2. Generational Distance (GD)

van Veldhuizen and Lamont [35] proposed the concept of generational distance (GD). The metric calculates the distance between P and P
*. It can be defined as follows:
(11)
GD
=
∑
x
∈
P
m
i
n
_
e
d
(
x
,
p
*
)
2
|
P
|
,
where m
i
n
_
e
d
(
x
,
P
*
) is the minimum Euclidean distance of solution x and the solutions in P
*. GD evaluates the closeness of the obtained Pareto front and the Pareto-optimal front. Algorithms with smaller value of GD are desirable.


## 4.2.3. Spread(
Δ
)

The spread metricΔ was firstly proposed by Deb et al. [2] to measure the spread effect of the nondominated solutions obtained by the algorithm. However, it is only suitable for biobjective problems. Wang et al. [19] further extended the metric Δ and made it suitable to measure problems that have more than two objectives. The improved Δ can be stated mathematically as follows:
(12)
Δ
=
∑
i
=
1
m
d
(
E
i
,
Ω
)
+
∑
Χ
∈
Ω
|
d
(
X
,
Ω
)
-
d
-
|
∑
i
=
1
m
d
(
E
i
,
Ω
)
+
(
|
Ω
|
-
m
)
d
-
d
(
X
,
Ω
)
=
min
⁡
Y
∈
Ω
,
Y
≠
X
⁡
∥
F
(
X
)
-
F
(
Y
)
∥
d
-
=
1
|
Ω
|
∑
X
∈
Ω
d
(
X
,
Ω
)
,
where Ω is a set of solutions, m is the number of objectives, and (E
1
,
…
,
E
m) are m extreme solutions in the set of true Pareto front (PF). The smaller the value of Δ, the better the distribution and diversity of the obtained nondominated solutions.


## 4.3. Experimental Setup

In this paper, the population size and the number of generations are set to 100 and 250, respectively. For the other parameters in NSGA-II and MODEA, they have remained the same with their original studies [3, 23]. For MODESA, there are four parameters: T
max
⁡, T
min
⁡, γ, and m
a
x
_
l
i
f
e
_
c
i
c
l
e that need to be set. Here we set T
max
⁡
=
100, T
min
⁡
=
1
e
-
7, γ
=
0.6, and m
a
x
⁡
_
l
i
f
e
_
c
y
c
l
e
=
1 for the study. Each algorithm runs 10 times independently for each problem. All the test problems have been executed on Microsoft Window XP Intel core 2 Duo CPU E8400 3.00 GHz, with 2 GB RAM. All methods have been implemented in C++.


## 4.4. Discussion of the Results

### 4.4.1. Results of Biobjective Optimization Problems

Tables3 and 4 show the mean and standard deviation IGD-metric and Δ-metric values of the obtained nondominated solutions in 10 runs on the group of ZDTs. It is obvious from Tables 3 and 4 that the IGD-metric values of MODESA achieve the best value in the three algorithms except for the test problem ZDT4.Table 3
The mean and standard deviation IGD-metric values of nondominated solutions in 10 runs on ZDTs whereN
=
100 and T
=
250.
IGD
NSGA-II (mean ± std)
MODEA (mean ± std)
MODESA (mean ± std)
ZDT1
0.00564974
±
0.000254893
0.00436863
±
0.000243648
0.00403528
±
7.01666
e
-
05
ZDT2
0.00557228
±
0.000242512
0.00451929
±
0.000134666
0.00408871
±
4.65746
e
-
05
ZDT3
0.00592808
±
0.000205118
0.00521698
±
0.000159453
0.00494596
±
0.000126657
ZDT4
0.0069977
±
0.00102386
0.0301279
±
0.0752722
0.0171826
±
0.0362604
ZDT6
0.0137079
±
0.00170579
0.00304881
±
0.000125906
0.0022843
±
0.00014073Table 4
The mean and standard deviationΔ-metric values of nondominated solutions in 10 runs on ZDTs where N
=
100 and T
=
250.
Spread
NSGA-II (mean ± std)
MODEA (mean ± std)
MODESA (mean ± std)
ZDT1
0.372425
±
0.0513842
0.333973
±
0.0772084
0.297395
±
0.0308609
ZDT2
0.385209
±
0.0291848
0.336548
±
0.0181084
0.296842
±
0.0282846
ZDT3
0.389909
±
0.0557332
0.355518
±
0.0370526
0.266801
±
0.0211668
ZDT4
0.399555
±
0.0786079
0.377719
±
0.0574188
0.352871
±
0.107095
ZDT6
0.407279
±
0.0276428
0.407335
±
0.0487992
0.368972
±
0.04
63564The test problems ZDT1 and ZDT2 may be the most simple problems in the series of ZDTs. The Pareto-optimal front of ZDT1 is convex while that of ZDT2 is nonconvex. It is clear from Tables3 and 4 that both the values of the IGD metric and the spread metric of MODESA are better than that of NSGA-II and MODEA. Figures 1 and 2 demonstrate that MODESA converges to the Pareto-optimal front and has a good spread over the entire front on the test problems ZDT1 and ZDT2. It can be concluded that the proposed MODESA can well deal with these two kinds of test problems.Figure 1
Plots of the nondominated solutions with the lowest IGD-metric values found by MODESA in 10 runs in the objective space on ZDT1 whereN
=
100 and T
=
250.Figure 2
Plots of the nondominated solutions with the lowest IGD-metric values found by MODESA in 10 runs in the objective space on ZDT2 whereN
=
100 and T
=
250.The third test problem ZDT3 is different from the above two test problems. The Pareto-optimal front of ZDT3 consists of five disjoint curves. From Figure3, it is obvious that the nondominated solutions obtained by the MODESA almost converge to the Pareto-optimal front and have a good distribution over the entire front. From Tables 3 and 4, MODESA obtains the best values of IGD-metric and spread-metric between the three algorithms.Figure 3
Plots of the nondominated solutions with the lowest IGD-metric values found by MODESA in 10 runs in the objective space on ZDT3 whereN
=
100 and T
=
250.The fourth test problem is ZDT4 which is perhaps the most difficult problem in the five problems. It has21
9 different local Pareto-optimal fronts that confuse the multiobjective optimization algorithm [2]. From Tables 3 and 4, NSGA-II obtains the best values on the IGD-metric and the spread-metric in these three algorithms. In fact, in the 10 independent runs, some results of MODEA and MODESA are better than that of NSGA-II; however, MODEA and MODESA get bad results in a few runs which affects the mean value of the 10 runs. Indeed, the best IGD-metric value of MODESA in the 10 runs is 0.0044915, which is better than NSGA-II. Figure 4 demonstrates the plots of the nondominated solutions with the lowest IGD-metric values found by MODESA in the 10 runs. It is clear from Figure 4 that MODESA can get quite good convergence and distribution over the Pareto-optimal front, which reveals that the proposed algorithm is a promising but unstable algorithm for solving the test problem ZDT4. We will leave this problem for future work.Figure 4
Plots of the nondominated solutions with the lowest IGD-metric values found by MODESA in 10 runs in the objective space on ZDT4 whereN
=
100 and T
=
250.The last test problem ZDT6 is such a problem with thin density and nonuniformed spread of solutions [23]. We can observe from Tables 3 and 4 that MODESA finds the best convergence and spread in comparison to the other two algorithms. To show the superiority of MODESA, plots of the nondominated solutions found by MODESA and MODEA on ZDT 6 are shown in Figures 5 and 6, respectively. It is clear from Figures 5 and 6 that MODESA gets a better spread than MODEA.Figure 5
Plots of the nondominated solutions with the lowest IGD-metric values found by MODESA in 10 runs in the objective space on ZDT6 whereN
=
100 and T
=
250.Figure 6
Plots of the nondominated solutions with the lowest IGD-metric values found by MODEA in 10 runs in the objective space on ZDT6 whereN
=
100 and T
=
250.


### 4.4.2. Results of Triobjective Optimization Problems

In this section, to show the effectiveness of the proposed MODESA for solving problems that have more than two objectives, we choose two triobjective optimization problems DTLZ1 and DTLZ2 for comparison. The results of NSGA-II, MODEA, and MODESA are given in Table5. It is obvious from Table 5 that MODESA performs better than the other two algorithms in terms of both the IGD-metric and spread-metric. Figures 7, 8, 9, and 10 show the nondominated solutions obtained by MODEA and MODESA, and it is obvious that MODESA outperforms MODEA on both of the two triobjective optimization problems. In particular, it is clear that MODESA has a better distribution than MODEA, the reason is that the vicinity distance based pruning method performs better than the crowding distance of Deb in maintaining diversity and distribution of the nondominated solutions.Table 5
The mean and standard deviation IGD-metric andΔ-metric values of nondominated solutions in 10 runs on DTLZ1 and DTLZ2 where N
=
100 and T
=
250.
Metric
IGD
Spread
Problem
DTLZ1
DTLZ2
DTLZ1
DTLZ2
NSGA-II (mean± std)
0.025008
±
0.001199
0.077375
±
0.021572
0.60056
±
0.076489
0.530832
±
0.029126
MODEA (mean± std)
0.022851
±
0.000557
0.060254
±
0.001664
0.41962
±
0.019617
0.425707
±
0.033528
MODESA (mean± std)
0.020602
±
0.00035
0.055011
±
0.000799
0.26489
±
0.045448
0.260335
±
0.02856Figure 7
Plots of the nondominated solutions with the lowest IGD-metric values found by MODEA in 10 runs in the objective space on DTLZ1 whereN
=
100 and T
=
250.Figure 8
Plots of the nondominated solutions with the lowest IGD-metric values found by MODESA in 10 runs in the objective space on DTLZ1 whereN
=
100 and T
=
250.Figure 9
Plots of the nondominated solutions with the lowest IGD-metric values found by MODEA in 10 runs in the objective space on DTLZ2 whereN
=
100 and T
=
250.Figure 10
Plots of the nondominated solutions with the lowest IGD-metric values found by MODESA in 10 runs in the objective space on DTLZ2 whereN
=
100 and T
=
250.


## 4.4.1. Results of Biobjective Optimization Problems

Tables3 and 4 show the mean and standard deviation IGD-metric and Δ-metric values of the obtained nondominated solutions in 10 runs on the group of ZDTs. It is obvious from Tables 3 and 4 that the IGD-metric values of MODESA achieve the best value in the three algorithms except for the test problem ZDT4.Table 3
The mean and standard deviation IGD-metric values of nondominated solutions in 10 runs on ZDTs whereN
=
100 and T
=
250.
IGD
NSGA-II (mean ± std)
MODEA (mean ± std)
MODESA (mean ± std)
ZDT1
0.00564974
±
0.000254893
0.00436863
±
0.000243648
0.00403528
±
7.01666
e
-
05
ZDT2
0.00557228
±
0.000242512
0.00451929
±
0.000134666
0.00408871
±
4.65746
e
-
05
ZDT3
0.00592808
±
0.000205118
0.00521698
±
0.000159453
0.00494596
±
0.000126657
ZDT4
0.0069977
±
0.00102386
0.0301279
±
0.0752722
0.0171826
±
0.0362604
ZDT6
0.0137079
±
0.00170579
0.00304881
±
0.000125906
0.0022843
±
0.00014073Table 4
The mean and standard deviationΔ-metric values of nondominated solutions in 10 runs on ZDTs where N
=
100 and T
=
250.
Spread
NSGA-II (mean ± std)
MODEA (mean ± std)
MODESA (mean ± std)
ZDT1
0.372425
±
0.0513842
0.333973
±
0.0772084
0.297395
±
0.0308609
ZDT2
0.385209
±
0.0291848
0.336548
±
0.0181084
0.296842
±
0.0282846
ZDT3
0.389909
±
0.0557332
0.355518
±
0.0370526
0.266801
±
0.0211668
ZDT4
0.399555
±
0.0786079
0.377719
±
0.0574188
0.352871
±
0.107095
ZDT6
0.407279
±
0.0276428
0.407335
±
0.0487992
0.368972
±
0.04
63564The test problems ZDT1 and ZDT2 may be the most simple problems in the series of ZDTs. The Pareto-optimal front of ZDT1 is convex while that of ZDT2 is nonconvex. It is clear from Tables3 and 4 that both the values of the IGD metric and the spread metric of MODESA are better than that of NSGA-II and MODEA. Figures 1 and 2 demonstrate that MODESA converges to the Pareto-optimal front and has a good spread over the entire front on the test problems ZDT1 and ZDT2. It can be concluded that the proposed MODESA can well deal with these two kinds of test problems.Figure 1
Plots of the nondominated solutions with the lowest IGD-metric values found by MODESA in 10 runs in the objective space on ZDT1 whereN
=
100 and T
=
250.Figure 2
Plots of the nondominated solutions with the lowest IGD-metric values found by MODESA in 10 runs in the objective space on ZDT2 whereN
=
100 and T
=
250.The third test problem ZDT3 is different from the above two test problems. The Pareto-optimal front of ZDT3 consists of five disjoint curves. From Figure3, it is obvious that the nondominated solutions obtained by the MODESA almost converge to the Pareto-optimal front and have a good distribution over the entire front. From Tables 3 and 4, MODESA obtains the best values of IGD-metric and spread-metric between the three algorithms.Figure 3
Plots of the nondominated solutions with the lowest IGD-metric values found by MODESA in 10 runs in the objective space on ZDT3 whereN
=
100 and T
=
250.The fourth test problem is ZDT4 which is perhaps the most difficult problem in the five problems. It has21
9 different local Pareto-optimal fronts that confuse the multiobjective optimization algorithm [2]. From Tables 3 and 4, NSGA-II obtains the best values on the IGD-metric and the spread-metric in these three algorithms. In fact, in the 10 independent runs, some results of MODEA and MODESA are better than that of NSGA-II; however, MODEA and MODESA get bad results in a few runs which affects the mean value of the 10 runs. Indeed, the best IGD-metric value of MODESA in the 10 runs is 0.0044915, which is better than NSGA-II. Figure 4 demonstrates the plots of the nondominated solutions with the lowest IGD-metric values found by MODESA in the 10 runs. It is clear from Figure 4 that MODESA can get quite good convergence and distribution over the Pareto-optimal front, which reveals that the proposed algorithm is a promising but unstable algorithm for solving the test problem ZDT4. We will leave this problem for future work.Figure 4
Plots of the nondominated solutions with the lowest IGD-metric values found by MODESA in 10 runs in the objective space on ZDT4 whereN
=
100 and T
=
250.The last test problem ZDT6 is such a problem with thin density and nonuniformed spread of solutions [23]. We can observe from Tables 3 and 4 that MODESA finds the best convergence and spread in comparison to the other two algorithms. To show the superiority of MODESA, plots of the nondominated solutions found by MODESA and MODEA on ZDT 6 are shown in Figures 5 and 6, respectively. It is clear from Figures 5 and 6 that MODESA gets a better spread than MODEA.Figure 5
Plots of the nondominated solutions with the lowest IGD-metric values found by MODESA in 10 runs in the objective space on ZDT6 whereN
=
100 and T
=
250.Figure 6
Plots of the nondominated solutions with the lowest IGD-metric values found by MODEA in 10 runs in the objective space on ZDT6 whereN
=
100 and T
=
250.


## 4.4.2. Results of Triobjective Optimization Problems

In this section, to show the effectiveness of the proposed MODESA for solving problems that have more than two objectives, we choose two triobjective optimization problems DTLZ1 and DTLZ2 for comparison. The results of NSGA-II, MODEA, and MODESA are given in Table5. It is obvious from Table 5 that MODESA performs better than the other two algorithms in terms of both the IGD-metric and spread-metric. Figures 7, 8, 9, and 10 show the nondominated solutions obtained by MODEA and MODESA, and it is obvious that MODESA outperforms MODEA on both of the two triobjective optimization problems. In particular, it is clear that MODESA has a better distribution than MODEA, the reason is that the vicinity distance based pruning method performs better than the crowding distance of Deb in maintaining diversity and distribution of the nondominated solutions.Table 5
The mean and standard deviation IGD-metric andΔ-metric values of nondominated solutions in 10 runs on DTLZ1 and DTLZ2 where N
=
100 and T
=
250.
Metric
IGD
Spread
Problem
DTLZ1
DTLZ2
DTLZ1
DTLZ2
NSGA-II (mean± std)
0.025008
±
0.001199
0.077375
±
0.021572
0.60056
±
0.076489
0.530832
±
0.029126
MODEA (mean± std)
0.022851
±
0.000557
0.060254
±
0.001664
0.41962
±
0.019617
0.425707
±
0.033528
MODESA (mean± std)
0.020602
±
0.00035
0.055011
±
0.000799
0.26489
±
0.045448
0.260335
±
0.02856Figure 7
Plots of the nondominated solutions with the lowest IGD-metric values found by MODEA in 10 runs in the objective space on DTLZ1 whereN
=
100 and T
=
250.Figure 8
Plots of the nondominated solutions with the lowest IGD-metric values found by MODESA in 10 runs in the objective space on DTLZ1 whereN
=
100 and T
=
250.Figure 9
Plots of the nondominated solutions with the lowest IGD-metric values found by MODEA in 10 runs in the objective space on DTLZ2 whereN
=
100 and T
=
250.Figure 10
Plots of the nondominated solutions with the lowest IGD-metric values found by MODESA in 10 runs in the objective space on DTLZ2 whereN
=
100 and T
=
250.


## 4.5. Effect of the Proposed Simulated Annealing Procedure

As mentioned in Section3.5.3, one of the main differences of MODEA and MODESA is the employment of the proposed simulated annealing procedure. To show the effect of the proposed simulated annealing procedure, a version of MODESA: MODESA-CD (multiobjective differential evolution with crowding distance) is proposed for comparison. The MODESA-CD is analogous to MODESA except that it replaces the vicinity distance based pruning method in MODESA with the crowding distance based approach in NSGA-II. In other words, the difference between MODEA and MODESA-CD is the using of the proposed simulated annealing procedure. Table 6 shows the mean and standard deviation GD-metric values of nondominated solutions in 10 runs where N
=
100 and T
=
250. It is obvious from Table 6 that all the GD-metric values of MODESA-CD are better than that of MODEA. The reason of the bad performance of MODEA is that MODEA always accepts the current and the new solutions to the next generation and that makes it chaotic and unable to have a right direction towards the Pareto-optimal front. In MODESA-CD, we remove the bad solutions and only accept some potential solutions in an alterable probability.Table 6
The mean and standard deviation GD-metric values of nondominated solutions in 10 runs whereN
=
100 and T
=
250.
GD
MODEA (mean ± std)
MODESA-CD (mean ± std)
ZDT1
0.000124
±
8.236
e
-
06
0.000116
±
1.437
e
-
05
ZDT2
4.819
e
-
05
±
1.642
e
-
06
4.771
e
-
05
±
2.596
e
-
06
ZDT3
0.00035
±
6.128
e
-
06
0.000345
±
2.714
e
-
05
ZDT4
0.002812
±
0.008057
0.001365
±
0.003738
ZDT6
0.000429
±
1.78
e
-
05
0.000404
±
2.614
e
-
05
DTLZ1
0.000308
±
7.454
e
-
05
0.00029
±
2.129
e
-
05
DTLZ2
0.000677
±
0.000103
0.000633
±
3.871
e
-
05


## 5. Conclusion

In this paper, an enhanced differential evolution based algorithm with simulated annealing, named MODESA, is presented for solving MOPs. During the selection, MODESA employs a procedure of simulated annealing to control the acceptance of every candidate solution, and a concept of the prior life circle is proposed to allow some potential solutions entering the next generation to escape from the local optimum. Finally, a fast and efficient diversity maintenance mechanism is adopted. The proposed algorithm is tested on five biobjective and two triobjective optimization problems in terms of IGD-metric and spread metric, and the experimental results show that MODESA outperforms the other two algorithms on these test problems except on ZDT4. Furthermore, the effect of convergence is also tested in comparison to MODEA, the computational results show the superiority of the proposed MODESA.Improvement can be made further in the performance of the proposed MODESA on the test problem ZDT4. In a word, MODESA is an easy and efficient method for solving MOPs. In the near future, applying MODESA for constrained MOPs and real-life application problems are our work.



---

*Source: 931630-2014-05-07.xml*