Self-adaptive Differential Evolution for Multi-objective Optimization with Local Search and Indicator-based Selection

1. Introduction

Multi-objective optimization problems (MOPs), characterized by their inherent complexity and conflicting objectives, necessitate the development of efficient strategies to balance these diverse criteria. In response to this challenge, a multitude of algorithms have been devised, each aiming to effectively address MOPs. Notably, Differential Evolution (DE) [1] has emerged as a promising contender within the realm of multi-objective optimization, exhibiting remarkable performance and demonstrating its potential to efficiently navigate the complex landscape of conflicting objectives. A detailed survey of some DE-based MOEAs has been conducted [2], discussing their advantages and disadvantages. Most of these methods employ the widely used DE/rand/1/bin variant. While they have demonstrated effectiveness in accelerating the search and enhancing results, the review also reveals a latent issue with DE: the tendency to become trapped in local optimum fronts. To address this, diversity maintenance schemes must be incorporated into DE-based MOEAs. Additionally, handling the non-dominance relation between the trial individual and the target individual poses a challenge in these algorithms.

To tackle the aforementioned challenges, a novel algorithm, ILSDEMO, is proposed in this paper. It is built upon the "DE/best/1/bin" and "DE/best/1/exp" variants [1] of DE. Prior experimental validation by E. Mezura-Montes [3] revealed that, regardless of the problem characteristics, the "DE/best/1/bin" variant exhibited the highest competitiveness in terms of both quality and robustness of results among eight DE variants. Different from the "DE/best/1/bin" variant, the "DE/best/1/exp" variant allows larger mutation blocks to be inherited, which helps the algorithm escape from local optimal solutions on certain problems and enhances its global search capability. While several researchers have attempted to leverage these variants for solving MOPs [4,5], our approach distinguishes itself in several key aspects. Firstly, ILSDEMO initializes the parent population using the orthogonal design method [6], which ensures a diverse initial population. Secondly, an archive population is employed to store non-dominated solutions, thereby maintaining a repository of promising candidates. Additionally, a k-nearest-neighbor rule is adopted to guarantee a uniformly distributed archive population in the objective space, enhancing the algorithm's ability to explore diverse regions of the search space. To generate a new parent population, an indicator-based selection mechanism is utilized, which leverages the performance indicators to identify promising individuals. Furthermore, the two DE variants and their corresponding parameters are selected in a self-adaptive manner, enabling the algorithm to adapt to the specific characteristics of the optimization problem. Lastly, a local search operator, specifically the q-Gaussian mutation[7,8], is applied to refine trial individuals that demonstrate superior performance compared to target individuals. This refinement step further enhances the algorithm's ability to find optimal solutions.

The remainder of the paper is organized as follows. Section 2 defines several notable concepts pertaining to multi-objective optimization problems (MOPs) and briefly introduces the principles of indicator-based selection, differential evolution, and q-Gaussian mutation. In Section 3, a comprehensive account of the ILSDEMO algorithm is provided, including a detailed description of its key components. Section 4 presents an experimental investigation of the algorithm on the ZDT-series and DTLZ-series benchmarks, along with a comparative analysis of the results obtained with selected multi-objective evolutionary algorithms. Section 5 summarizes the key findings and conclusions of the paper.

2. Preliminaries

Before delving into the details of the ILSDEMO algorithm, it is essential to introduce some fundamental concepts and theories that provide the necessary background for the algorithm's comprehension. These preliminary sections not only lay the foundation for understanding the principles and design of ILSDEMO but also aid in appreciating the algorithm's strengths and characteristics in solving multi-objective optimization problems. In the following sections, multi-objective optimization, orthogonal design, indicator-based selection, differential evolution, and Q-Gaussian mutation are introduced sequentially.

2.1. Multi-Objective Optimization

In contrast to single-objective optimization problems (SOPs), MOPs involves the optimization of multiple objectives. For the purpose of this discussion, it is assumed that all objectives are to be minimized. The problem can be formulated using Formula 1.

$$\begin{array}{l}Minimize\hspace{0.33em}Z(\stackrel{\rightharpoonup }{X})=\{z_1=f_1(\stackrel{\rightharpoonup }{X}),\hspace{0.33em}{z}_2=f_2(\stackrel{\rightharpoonup }{X}),\hspace{0.33em}\mathrm{...},\hspace{0.33em}{z}_m=f_m(\stackrel{\rightharpoonup }{X})\}\\ s.t.\hspace{0.33em}{g}_i(\stackrel{\rightharpoonup }{X})\le 0,\mathrm{i}=1,2,\mathrm{...},r\end{array}$$

(1)

where $\stackrel{\rightharpoonup }{X}\in R^n$ is a vector with n decision variables, $f_1(\stackrel{\rightharpoonup }{X})$, $f_2(\stackrel{\rightharpoonup }{X})$, …, $f_m(\stackrel{\rightharpoonup }{X})$ are objective functions, $g_1(\stackrel{\rightharpoonup }{X})$, $g_2(\stackrel{\rightharpoonup }{X})$, …, $g_r(\stackrel{\rightharpoonup }{X})$ are constraint functions. The decision space can be denoted as $S=\{\stackrel{\rightharpoonup }{X}\in R^n|g_i(\stackrel{\rightharpoonup }{X})\le 0,i=1,\hspace{0.33em}2,\hspace{0.33em}\mathrm{...},\hspace{0.33em}\hspace{0.33em}r\}$. A solution in S is called a feasible solution. The objective space can be represented as $\hspace{0.33em}Z=\{\stackrel{\rightharpoonup }{z}\in R^m|z_1=f_1(\stackrel{\rightharpoonup }{X}),\hspace{0.33em}{z}_1=f_2(\stackrel{\rightharpoonup }{X}),\hspace{0.33em}\mathrm{...},\hspace{0.33em}{z}_m=f_m(\stackrel{\rightharpoonup }{X}),\hspace{0.33em}\stackrel{\rightharpoonup }{X}\in S\}$.

For the sake of facilitating the forthcoming discussion, we hereby present a comprehensive compilation of mathematical terminologies in MOO below.

Definition 1 (Pareto Dominance,$\succ$)

$${\stackrel{\rightharpoonup }{X}}_1\succ {\stackrel{\rightharpoonup }{X}}_2\hspace{0.33em}({\stackrel{\rightharpoonup }{X}}_1\in S,\hspace{0.33em}{\stackrel{\rightharpoonup }{X}}_2\in S)\hspace{0.33em}\hspace{0.33em}iff\hspace{0.33em}\forall i\in \{1,\hspace{0.33em}2,\hspace{0.33em}\mathrm{...},\hspace{0.33em}m\}:f_i({\stackrel{\rightharpoonup }{X}}_1)\le f_i({\stackrel{\rightharpoonup }{X}}_2)\hspace{0.33em}\wedge \exists j\in \{1,\hspace{0.33em}2,\hspace{0.33em}\mathrm{...},\hspace{0.33em}m\}:\hspace{0.33em}{f}_j({\stackrel{\rightharpoonup }{X}}_1)<f_j({\stackrel{\rightharpoonup }{X}}_2)$$

Definition 2 (Non-dominated, ~)

$${\stackrel{\rightharpoonup }{X}}_1~{\stackrel{\rightharpoonup }{X}}_2\hspace{0.33em}({\stackrel{\rightharpoonup }{X}}_1\in S,\hspace{0.33em}{\stackrel{\rightharpoonup }{X}}_2\in S)\hspace{0.33em}\hspace{0.33em}iff\hspace{0.33em}\exists i,\hspace{0.33em}j\in \{1,\hspace{0.33em}2,\hspace{0.33em}\mathrm{...},\hspace{0.33em}m\}:f_i({\stackrel{\rightharpoonup }{X}}_1)<f_i({\stackrel{\rightharpoonup }{X}}_2)\hspace{0.33em}\wedge f_j({\stackrel{\rightharpoonup }{X}}_1)>f_j({\stackrel{\rightharpoonup }{X}}_2)$$

Definition 3 (Non-dominated set, NDS)

The non-dominated set, denoted as $NDS=\left\{{\stackrel{\rightharpoonup }{X}}_{*}\right|iff\hspace{0.33em}\neg \exists \stackrel{\rightharpoonup }{X}\in NDS:\stackrel{\rightharpoonup }{X}\succ {\stackrel{\rightharpoonup }{X}}_{*}\}$

Definition 4 (Pareto optimal solution)

A feasible solution ${\stackrel{\rightharpoonup }{X}}_{*}\in S$ is called a Pareto optimal solution, $iff\hspace{0.33em}\neg \exists \stackrel{\rightharpoonup }{X}\in S:\stackrel{\rightharpoonup }{X}\succ {\stackrel{\rightharpoonup }{X}}_{*}$

Definition 5 (Pareto optimal set, PS)

The Pareto optimal set, denoted as$PS=\{{\stackrel{\rightharpoonup }{X}}_{*}\in S|iff\hspace{0.33em}\neg \exists \stackrel{\rightharpoonup }{X}\in S:\stackrel{\rightharpoonup }{X}\succ {\stackrel{\rightharpoonup }{X}}_{*}\}$, is the set of all Pareto optimal solutions.

Definition 6 (Pareto front, PF)

The Pareto front is defined as$PF=\{\stackrel{\rightharpoonup }{z}(\stackrel{\rightharpoonup }{X})\in Z|\stackrel{\rightharpoonup }{X}\in PS\}$.

In accordance with definition 5, PS is the ultimate goal of optimization. However, the PS, especially for continuous MOPs, cannot be completely represented owing to the enormous amount of solutions. Moreover, too many solutions make no sense to help decision makers. Therefore, an appropriate selection for MOPs is to find a non-dominated set representing the Pareto optimal set as far as possible.

2.2. Orthogonal Design

The experimental design is an efficient approach for organizing experiments to enable the analysis of obtained data, leading to valid and objective conclusions [9,10]. In experimental designs, the decision variables are designated as factors, and the range of each factor is partitioned into multiple levels. Subsequently, all possible combinations of factors across all levels are tested, ultimately leading to the selection of the optimal combination. Given n factors at q levels, the total number of combinations is calculated as $q^n$. When n and q are small, evaluating all combinations may be feasible. However, for a large number of factors and levels, testing all combinations becomes impractical. Therefore, a representative subset of combinations is selected for evaluation. Orthogonal design (OD) was developed for this purpose and has been successfully employed in evolutionary algorithms [11,12].

In OD, an orthogonal table provides some representative combinations of the factors at different levels. Let L_m(q, c) be an orthogonal table for c factors and q levels, where L denotes a Latin square and m is the number of the constructed combinations of the levels ($m=q^j,c=\frac{q^j-1}{q-1},j>0$). In the orthogonal table with m rows and c columns, each column denotes a factor and each row indicates a combination of the levels regarding the factors. The orthogonal table has some characteristics: 1) every level has the same occurrence number, i.e., m/q times, in any column; 2) every combination of two levels with respect to two factors has the same occurrence number, i.e., m/q²; 3) the selected combinations are uniformly distributed over the whole space of all the possible combinations; 4) if any two columns are swapped, the altered table is still an orthogonal table; 5) if some columns are removed from the orthogonal, the modified table is still an orthogonal table. A positive integer j should be found to minimize m with a constraint c ≥ n (n is the number of decision variables).

In this study, the algorithm described in [13] is employed to produce an orthogonal table L_m(q, c). Each element of the table, a (i, j), indicates the level of the jth factor in the ith combination in L_m(q, c). If c>n, the last c-n columns will be deleted to get an orthogonal table with n factors. An element of the table is calculated as formula 2.

$$a(i,j)=\left\{\begin{array}{c}\hspace{0.33em}{l}_j\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}i=1\\ l_j+(i-1)(\frac{u_j-l_j}{q-1})\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}2\le i\le q-1\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\\ u_j \hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}q\hfill \end{array}\right.$$

(2)

Here, l_j and u_j are the lower bound and the upper bound of the jth decision variable respectively, and the domain is quantized q-1 fractions.

2.3. Indicator-Based Selection

Indicator-based selection emerged from IBEA was first proposed by Zitzler [13]. IBEA is based on quality indicators where a function I assigns each Pareto set approximation a real value reflecting its quality, and thus the optimization goal becomes the identification of a Pareto set approximation that minimizes (or maximizes) I. IBEA only compares pairs of individuals instead of entire approximation sets. The main advantage of the indicator concept is that no additional diversity preservation mechanisms are required.

A binary quality indicator is a function that maps multiple Pareto set approximations to a real number, and can be used to compare the quality of two Pareto set approximations. The indicator also can be used to compare two single solutions, and thus can serve for the selection process of evolutionary algorithms. As suggested in the literature [18,14], the majority of the multi-objective ranking techniques in the literature are easily adaptable into binary quality indicators. Here, the epsilon indicator $I_{\epsilon }$, which will be utilized in the algorithm, is solely introduced. The indicator quantifies the difference in quality between two solutions. It is defined as formula 3.

$$I_{\epsilon }({\stackrel{\rightharpoonup }{X}}_1,{\stackrel{\rightharpoonup }{X}}_2)=\max (f_i({\stackrel{\rightharpoonup }{X}}_1)-f_i({\stackrel{\rightharpoonup }{X}}_2)),\hspace{0.33em}i\in \{1,\mathrm{...},m\}$$

(3)

In order to assign fitness to the individuals, one of the following two schemes can be chosen: 1) One approach is to simply sum up the indicator values for each individual paired with the rest of the population as shown in formula 4. 2) Another approach is to amplify the influence of dominating individuals over dominated ones as shown in formula 5. Here, the parameter k is a scaling factor depending on the indicator and the problem.

$$F({\stackrel{\rightharpoonup }{X}}_1)={\displaystyle {\sum }_{{ }_{{\stackrel{\rightharpoonup }{X}}_2\in P\backslash \left\{{\stackrel{\rightharpoonup }{X}}_1\right\}}}I({\stackrel{\rightharpoonup }{X}}_2,\hspace{0.33em}{\stackrel{\rightharpoonup }{X}}_1)}$$

(4)

$$F({\stackrel{\rightharpoonup }{X}}_1)={\displaystyle {\sum }_{{ }_{{\stackrel{\rightharpoonup }{X}}_2\in P\backslash \left\{{\stackrel{\rightharpoonup }{X}}_1\right\}}}-e^{-I({\stackrel{\rightharpoonup }{X}}_2,\hspace{0.33em}{\stackrel{\rightharpoonup }{X}}_1)/k}}$$

(5)

The fitness values of the individuals produce a total ordering relation of the whole population P. As proved by zitzler [18], $I_{\epsilon }$ is dominance preserving. Thus, selection can be performed on the population P according to the orders of the fitness values.

Definition 7 (Dominance preserving) A binary quality indicator I is denoted as dominance preserving if:

1)$\forall {\stackrel{\rightharpoonup }{X}}_1,\hspace{0.33em}{\stackrel{\rightharpoonup }{X}}_2\in S,\hspace{0.33em}{\stackrel{\rightharpoonup }{X}}_1\succ {\stackrel{\rightharpoonup }{X}}_2\Rightarrow I({\stackrel{\rightharpoonup }{X}}_1,{\stackrel{\rightharpoonup }{X}}_2)<I({\stackrel{\rightharpoonup }{X}}_2,{\stackrel{\rightharpoonup }{X}}_1)$, and

2)$\forall {\stackrel{\rightharpoonup }{X}}_1,\hspace{0.33em}{\stackrel{\rightharpoonup }{X}}_2,{\stackrel{\rightharpoonup }{X}}_3\in S,\hspace{0.33em}{\stackrel{\rightharpoonup }{X}}_1\succ {\stackrel{\rightharpoonup }{X}}_2\Rightarrow I({\stackrel{\rightharpoonup }{X}}_3,{\stackrel{\rightharpoonup }{X}}_1)\ge I({\stackrel{\rightharpoonup }{X}}_3,{\stackrel{\rightharpoonup }{X}}_2)$

2.4. Differential Evolution

Differential evolution, proposed by Storn and Price [1], is one of the most popular meta-heuristic optimizers, capable of handling non-differentiable, nonlinear, high-dimensional, multimodal continuous optimization problems. It has a concise idea but a strong global convergence. DE generates offspring by perturbing solutions with scaled differences of randomly selected vectors instead of the operations of crossover and mutation used in conventional evolutionary algorithms [15].

The proposed algorithm will utilize two variants of DE, namely DE/best/1/bin and DE/best/1/exp, which use binomial crossover and exponential crossover respectively [1]. The main difference between binomial crossover and exponential crossover is the fact that the components inherited from the mutant vector are arbitrarily selected and usually scattered in binomial crossover while one gene block of the mutant vector is inherited in exponential crossover. The following two key program segments illustrate the procedures of them. Here, D is the dimension of the decision variable; randint(1, D) denotes generating a random integer number between 1 and D; r1 and r2 are two randomly generated integers with uniform distribution and their values are lower than or equal to the population size, and mutually different; best is a randomly generated integer with uniform distribution and less than the size of the non-dominated solution set; F and CR represent the scaling factor and the crossover probability respectively.

Program segment 1

jrand = randint (1, D)

for j=1 to D do

if rand(0,1)<CR || j == jrand then

Vji = Xjbest + F*(Xjr1-Xjr2)

else

Vji=Xji

endif

endfor

Program segment 2

j=randint (1, D)

k=0

do

Vji = Xjbest + F*(Xjr1-Xjr2)

j = (j+1) % D

k++

while (rand(0,1)<CR && k<D)

2.5. Q-Gaussian Mutation

The q-Gaussian distribution [16] allows to control the shape of the distribution by setting a real parameter q and can reproduce either finite second moment distributions, like the Gaussian distribution, or infinite second moment distributions, like the heavy tail Lévy distribution.

One of the most interesting properties of the Gaussian distribution is that it maximizes, under certain constraints, the Boltzman-Gibbs entropy which is defined as formula 6.

$$S={\displaystyle {\int }_{-\infty }^{+\infty }p(x)\ln (p(x))dx}$$

(6)

Because the Gaussian distribution cannot represent well correlated systems with an infinite second moment, Tsallis[17] proposed a generalized entropy form as formula 7.

$$S_q={\displaystyle \frac{1-{\displaystyle {\int }_{-\infty }^{+\infty }p{(x)}^{q}dx}}{q-1}}$$

(7)

The q-Gaussian distribution arises when maximizing the generalized entropy form given by formula 7. The parameter q controls the shape of the q-Gaussian distribution. The q-Gaussian distribution reproduces the usual Gaussian distribution for q=1. When q=(3+m)/(1+m) and 0<m<∞, it becomes a Student’s t-distribution with m degrees of freedom. The q-Gaussian reproduces the Cauchy distribution for q = 2. When q<1, the q-Gaussian distribution has a compact form. When -∞<q<3, the q-Gaussian distribution density is given as formula 8.

$$p(x;{\overline{\mu }}_q,{\overline{\sigma }}_q)=A_q\sqrt{B_q}{e}_q^{-B_q{(x-{\overline{\mu }}_q)}^2}$$

(8)

In formula 8, $e_q^x$ is defined as formula 9. The q-mean ${\overline{\mu }}_q$ and q-variance ${\overline{\sigma }}_q$ are defined as formula 10 and formula 11 respectively. $A_q$ is the normalization factor and $B_q$ controls the width of the q-Gaussian distribution and is given by $B_q={\left((3-q){\overline{\sigma }}_q^2\right)}^{-1}$. A random variable x following a q-Gaussian distribution with ${\overline{\mu }}_q$ and ${\overline{\sigma }}_q^2$ is denoted by $x~N_q({\overline{\mu }}_q,{\overline{\sigma }}_q^2)$. The generalized Box-Muller method can be used to generate q-Gaussian distribution for -∞<q<3. As described in [11], larger values of q result in longer tails of the q-Gaussian distribution. The algorithm given by [11] is employed to mutate a promising individual.

$$e_q^x\equiv \left\{\begin{array}{c}{(1+(1-q)x)}^{\frac{1}{1-q}}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}if(1+(1-q)x)\ge 0\\ 0\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}otherwise\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\end{array}\right.$$

(9)

$${\overline{\mu }}_q\equiv {\displaystyle \frac{{\displaystyle \int xp{(x)}^{q}dx}}{{\displaystyle \int p{(x)}^{q}dx}}}$$

(10)

$${\overline{\sigma }}_q^2\equiv {\displaystyle \frac{{\displaystyle \int (x-{\overline{\mu }}_q}{)}^{2}p{(x)}^{q}dx}{{\displaystyle \int p{(x)}^{q}dx}}}$$

(11)

3. The Proposed Method

3.1. The Motivation

In evolutionary algorithms, many metrics such as convergence, accuracy, speed, and stability are employed to measure the performance of the algorithm. Besides, uniformity of the PF also acts as an important performance metric in MOEAs. However, a single method cannot perform well in terms of all the metrics [18]. Take three popular algorithms for example. NSGAII performs well on two-objective problems in terms of convergence. Nevertheless, it usually cannot reach desirable results on three-objective problems [19]. SPEA2 also works well on two or three-objective problems, but has a deteriorated convergence and time performance with increasing number of objectives. PESAII[20] performs well in terms of convergence, but like NSGAII, it also relies on the quality of the initial population, and when dealing with large-scale or high-dimensional optimization problems, the computational cost and implementation difficulty of the algorithm may increase significantly.

In order to get an overall promotion in multiple metrics, some ways from different perspectives are worth trying. In Pareto rank based methods, the running time of non-dominated sorting and maintaining the diversity of the obtained non-dominated front all can be the starting points to improve the performance. Moreover, it is preferable to employ simple and efficient variation operators, such as DE and particle swarm optimization, to improve the time performance of the algorithm. In addition, using different selection operators to adjust the selection pressure and employing local search operators to improve the accuracy of the results are also popular in recent years [21]. Recently, self-adaptive parameters and strategies have received more and more attention due to their significant impacts on the performance of the algorithms.

Based on the above considerations, we propose a novel MOEA which mainly incorporates orthogonal design, self-adaptive DE, local search, indicator-based selection, and updating method using k-nearest-neighbor rule.

3.2. Population Initialization

The state of the initial population usually has an effect on the results. Experimental design methods such as uniform design and orthogonal design are statistically sound [12] and popular in many disciplines.

In our method, orthogonal design is utilized to generate an initial population where the individuals scattered uniformly over the feasible solution space, so that the algorithm can evenly scan the search space and locate the promising solutions quickly.

3.3. Variation Operators

3.3.1. Self-Adaptive Differential Evolution

In self-adaptive DE, not only the two control parameters but also the strategies are determined self-adaptively to avoid early stagnation. More specially, the control parameters and the strategy are encoded into the individuals. Initially, the scaling factor and crossover probability are randomly initialized within [0.2, 0.8] and [0.7, 1.0] respectively. One of the strategies, DE/best/1/bin or DE/best/1/exp, is randomly assigned to every individual. At the end of each iteration, the means and variances of two parameter values collected from the improved individuals are calculated to build two Gaussian models. The unimproved individuals will be assigned new parameter values by sampling the Gaussian models. Moreover, the strategies of these individuals will be reinitialized according to the probabilities of the two strategies based on the improved individuals.

3.3.2. Local Search

The purpose of local search methods is to enhance the exploitation capability of the algorithm. Previous researches have been focused on local search in MOO. In [22], a local search procedure was applied to new solutions but only a small number of neighborhood solutions were examined to prevent the local search procedure from spending almost all available computation time. In [23], a gradient-based mathematical programming method was employed to perform the local search on the offspring population. In [24], based on non-dominated sorting and crowding distance, local RBF (Radial Basis Function) surrogates were constructed around the sparse points, with the aim of accurately predicting the values of the objective function. In [25], a PF model-based local search method was proposed to accelerate the exploration and exploitation of the PF. Different from other local search methods in the literature of MOEA, q-Gaussian mutation is used to perform exploitation once the trial individual is better than the target individual. But the replacement scheme is noticeable in view of the fact that the comparison of the solutions for MOPs is distinct from that of SOPs. In this algorithm, if the individual obtained by local search dominates the target individual, the latter will be replaced by the former. Otherwise, the target individual remains unchanged. In either case, the trial individual should be used to update the archive population in order to avoid neglecting possible improvement on the distribution of the archive population.

3.4. The Selection Operator

In SOPs, the selection operator is usually employed to adjust the proportion of the best individuals involved in the variation of the next generation. In MOPs, it has a similar effect. However, besides the convergence, the diversity of the parent population in the objective space should also be concerned. In NSGAII, the uniformity of the parent population is only maintained by the crowding selection performed on the part of the population in the early stage. Thus, the crowding distance cannot effectively reflect the density of a solution especially in three or more objective problems. In fact, using the same selection operator on the archive population, NSGAII only achieves acceptable performance in terms of diversity in two-objective problems. Fortunately, the deficiency is compensated via propagating the extreme solutions to the next generation. In contrast, indicator-based selection can achieve better uniformity in the whole parent population. Nevertheless, indicator-based selection is prone to suffer a loss of the extreme solutions. Thus, the performance deteriorates with a worse coverage. Moreover, the algorithm cannot perform further exploitation to the extreme solutions.

In this study, indicator-based selection is employed to generate a new parent population based on a union of the old parent population and the archive population. First, the duplicated individuals are eliminated from the union. Then, the extreme solutions will be ensured to survive to the new parent population via assigning small enough indicator values to them. In addition, the formula 5 is converted into the formula 12 so as to keep the selection of indicator values in harmony with the objectives, i.e., minimization. Thus, if the number of the remaining individuals is not beyond the parent population size, all the individuals will be propagated into the next parent population. Moreover, the parent population may need to be supplemented by some new initialized individuals. Otherwise, the indicator-based selection will work. Specifically, an individual with a maximum indicator value is deleted from the union, and the indicator values of the other individuals are updated concurrently. The procedure is repeated until the number of the remaining individuals is equal to the parent population size.

$$F(x_{{ }^1})={\displaystyle {\sum }_{x_{{ }^2}\in P\backslash \left\{x_{{ }^1}\right\}}{e}^{-I(x_2,x_1)/({\max }_{x_{{ }_i},x_{{ }_j}\in S,\hspace{0.33em}i\ne j}|I(x_i,x_j)|\hspace{0.33em}*\hspace{0.33em}k)}}$$

(12)

3.5. The Updating Rule

Although indicator-based selection is an effective method in favor of finding a well-distributed solution set, it still cannot obtain ideal uniformity as well as some other algorithms such as SPEA2. Therefore, an updating criterion using the nearest neighbor will be introduced for the update of the archive population in both algorithms. In detail, when the number of non-dominated solutions exceeds the archive population size, the cumulative Euclidean distance (called density) to its k nearest neighbors for every non-dominated solution will be calculated. Here k is set to 2 * (m - 1), and m is the objective number. The individual with the smallest density will be removed. Meanwhile, the densities of other individuals are updated. Then, the procedure is repeated until the quantity of the rest is equal to the archive population size.

3.6. The Framework of the Proposed Algorithm

Based on the above design, the procedure of the proposed ILSDEMO algorithm for solving MOPs is summarized in Algorithm 1.

Algorithm 1 ILSDEMO

1: Set n=0. Initialize the population Pn{(Xi,Fi,CRi,Si),i∈[1,NP]} using OD 2: Evaluate Pn 3: Copy the non-dominated solutions to the archive population PA. 4: while the termination criteria are not satisfied do 5: Set n = n + 1 6: for i=1 to NP do 7:  select r1≠r2 ∈{1, 2, …, NP} randomly 8:  best∈{1,2,…,fsize} //fsize is the size of the non-dominated set in PA. 9:  Perform program segment 1 or 2 using Si,Fi,CRi 10:  Evaluate Vi 11:  if f(Vi) is non-dominated by f(Xi) then 12:   if f(Vi) dominates f(Xi) then 13:    Mark f(Xi) as an improved individual 14:   end if 15:   Perform q-Gaussian mutation on Vi 16:   Add Vi to PA using k-nearest neighbor rule 17:  end if 18:  if f(Vi) dominates f(Xi) then 19:   Xi = Vi 20:  end if 21: endfor 22: Build two Gaussian Models using the F and CR values collected from the individuals improved by DE 23: Calculate the probabilities of the strategies based on the individuals improved by DE 24: Generate Pn+1 based on Pn and PA using indicator-based selection. 25: Reinitialize F, CR, and S for the unimproved individuals according to the above models 26: endwhile

4. Experiments

In this section, the performance of ILSDEMO is experimentally investigated, compared with NSGAII, IBEA, and DEMO, which are renowned for their distinctive attributes, manifesting diverse strengths and applicability across various problem domains, commonly serving as benchmark algorithms for comparative analysis. The comparative assessment is conducted utilizing two comprehensive test suites and three rigorous performance metrics, with statistical analysis employed to discern the outcomes. All the experiments are implemented in VC6.0, Matlab 7.1, Sigmaplot 12.0 and run on Intel i5 760 2.8GHz machines with 4GB DDR3 RAM.

4.1. Performance Metrics

In general, the quality of the non-dominated sets should be evaluated according to the convergence and the diversity in the objective space. The convergence depicts the closeness of the final non-dominated solutions to the true PF, whereas the diversity aims at the uniformity and the coverage of the final solution set along the true PF. A number of performance metrics, such as hypervolume (HV) [26], generational distance (GD) [27], inverted generational distance (IGD) [31], spacing metric (Spacing) [28,29], set coverage (SC), among others, have been proposed. However, none of the performance metrics can reliably evaluate both performance goals [30]. Given that IGD and HV concurrently account for the convergence and diversity of the solution set[31], while the Spacing metric precisely captures the uniformity of the distribution of the solution set in the objective space, these three metrics are selected for comparing the non-dominated solutions obtained by the four algorithms. A brief description of the metrics is given as follows.

1) HV metric

HV metric, originally proposed by Zitzler and Thiele [32], computes the volume in the objective space covered by the non-dominated solutions and a reference point. Coello Coello, Van Veldhuizen and Lamont [33] described it as the Lebesgue measure $\mathrm{\Lambda }$ of the union of hypercubes enwrapped by a non-dominated set B and a bounding reference point z_r which is dominated by all points. The formalization is shown as formula 13. As elaborated by Knowles [34], the hypervolume metric is the only unary indicator that is capable of detecting that solution set A with a low indicator value is not better than B with a high value.

$$HV(B,\hspace{0.33em}{z}_r)=\mathrm{\Lambda }({\displaystyle \underset{z\in B}{\cup }\{z\text{'}|z\prec z\text{'}\prec z_r\})},\hspace{0.33em}B\subseteq R^m$$

(13)

In subsequent experiments, a fast algorithm termed HSO (Hypervolume by Slicing Objectives)[34] will be utilized for the computation of hypervolume, and the difference in hypervolume with respect to a reference set RS will be adopted as a performance indicator designated as HV*. The indicator value is defined as formula 14. Hence, a smaller value of HV* corresponds to a higher quality.

$$HV*=\hspace{0.33em}HV(RS,Zr)\hspace{0.33em}-\hspace{0.33em}HV(A,Zr)$$

(14)

2) IGD

Suppose P* and P represent a uniformly sampled solution set along the true PF and a non-dominated solution set obtained by an algorithm respectively. The mean distance from P* to P is calculated as formula 15:

$$D(P*,P)={\displaystyle \frac{{\displaystyle {\sum }_{{ }_{v\in P*}}d(v,P)}}{\left|P*\right|}}$$

(15)

where $d(v,P)$ indicates the minimum Euclidean distance between v and the points in P. If P* can represent the PF very well, the metric can measure both the convergence and the diversity while GD metric can only measure the convergence. A small value for IGD is desirable as it indicates better performance of the algorithm in terms of convergence and distribution.

3) Spacing

The Spacing metric measures how evenly the points in the approximation set are distributed in the objective space. An enhanced version of the metric will be employed, taking into account the Euclidean distances between the two extreme endpoints of the non-dominated front and their corresponding endpoints on the PF. A smaller spacing value indicates the solution set has a better spread. The metric is given by formula 16 where $d_f$ and $d_l$ are the Euclidian distances between the extreme solutions in PF and NDS, $d_i$ denotes the distance between the adjacent solutions via sorting one objective, and $\overline{d}$ denotes the average of $d_i$, $i\in [1,\hspace{0.33em}\left|NDS\right|-1].$Although the metric is utilized to measure the diversity of non-dominated set extensively, it only works well for two objective problems due to the calculation method of d_i.

$$\mathrm{\Delta }={\displaystyle \frac{d_f+d_l+{\displaystyle {\sum }_{i=1}^{\left|NDS\right|-1}\left|d_i-\overline{d}\right|}}{d_f+d_l+(\left|NDS\right|-1)\overline{d}}}$$

(16)

4.2. Parameters Settings

According to the preliminary experiments, the common parameters are set as follows: NP=50, NA=100 (150 for DTLZ-series), F∈[0.2, 0.8], CR∈[0.7,1.0], RUNs=30, Max_Evals=50000 (100000 for DTLZ-series). Here, NP represents the parent population size, NA denotes archive population size, F and CR are the scaling factor and crossover rate of DE respectively, RUNs and Max_Evals indicate the number of runs and evaluations in one run respectively. In addition, a factor k is set to 0.02 in indicator-based selection. In IBEA and NSGAII, NO (=50) indicates the offspring population size, PC (=0.8) and PM (=0.1) represent the probabilities of simulated binary crossover (SBX) and polynomial mutation respectively, eta_cross (=20) and eta_mut (=20) are the distribution indices of SBX and polynomial mutation respectively. In the computation of performance metrics, the objective values are normalized using the arctan function and subsequently mapped onto the interval [0, 1]. Additionally, the objective value of each reference point is set to 1.0.

4.3. Experimental Results and Discussions

In this section, the experimental results obtained by ILSDEMO, NSGAII, IBEA, and DEMO on the DTLZ-series and ZDT-series test instances are presented and discussed. The ZDT-series comprises two-objective problems, with the dimension of decision variables set to 30 for ZDT1, ZDT2, and ZDT3, and 10 for ZDT4 and ZDT6 in the conducted experiments. The DTLZ-series, on the other hand, represents scalable problems. For the DTLZ-series, the objective dimension is fixed at 3, with the dimension of decision variables set to 7 for DTLZ1, 12 for DTLZ2-6, and 22 for DTLZ7 respectively.

According to the HV* metric, Figure 1, Figure 2, and Figure 3 demonstrate the optimal approximations to the true PFs achieved by the four algorithms on the ZDT-series and DTLZ-series test instances.

Figure 1 reveals that the non-dominated solutions attained by ILSDEMO for the ZDT-series exhibit a more uniform distribution along the non-dominated front, indicating promising approximations across all algorithms. In contrast, NSGAII and DEMO exhibit seemingly inferior diversity, while IBEA maintains the worst distribution due to its inability to effectively retain extreme solutions, resulting in incomplete coverage of the non-dominated front.

As evident from Figure 2 and Figure 3, ILSDEMO achieves more accurate and uniform approximations to the true PFs of the DTLZ-series. Specifically, for DTLZ1, IBEA's solutions exhibit a more even distribution, while NSGAII maintains the worst distribution. For DTLZ2 and DTLZ4, NSGAII and DEMO display similar distributions, with IBEA exhibiting the worst distributions. While IBEA's solutions demonstrate a certain degree of uniformity in the centers of the fronts, significant gaps are observed between the centers and edges, along with the loss of some extreme points of the three objectives. For DTLZ3, NSGAII's approximate front is the worst due to outliers, while DEMO outperforms IBEA. For DTLZ5 and DTLZ6, ILSDEMO and DEMO yield similar distributions and accuracies, while IBEA produces the worst results. Finally, for DTLZ7, ILSDEMO approaches the true PF more evenly, while the non-dominated front obtained by NSGAII resembles that of DEMO.

To gain further insight into the performance of ILSDEMO, a comparison with the other three algorithms is conducted based on the convergence speed, which is quantified by the average HV* values as a function of the number of function evaluations. Figure 4 clearly demonstrates a consistent ranking in terms of convergence speed among the four algorithms: ILSDEMO exhibits the fastest convergence, followed by DEMO, NSGAII, and finally IBEA, which converges the slowest. Figure 5 also confirms that ILSDEMO maintains its position as the fastest algorithm. Specifically, for DTLZ1, DTLZ3, and DTLZ7, ILSDEMO converges significantly faster than the other algorithms. However, for DTLZ2, IBEA converges only slightly slower than ILSDEMO. Similarly, for DTLZ4 and DTLZ6, DEMO converges with a minor lag behind ILSDEMO. Lastly, for DTLZ5, NSGAII performs almost on par with ILSDEMO.

The box-plots presented in Figure 6, Figure 7, and Figure 8 offer statistical evidence of the variations in the performance of the four algorithms across three performance metrics. Notably, ILSDEMO's overall performance remains highly competitive compared to the other three methods in both test sets.

Analysis of Figure 6 reveals that, apart from ZDT1, IBEA demonstrates inferior performance on the ZDT-series, attributed to the loss of extreme solutions. For ZDT1, ZDT2, ZDT3, and ZDT4, NSGAII outperforms DEMO, albeit with some outliers exceeding the upper quartile in ZDT3 and ZDT4. Intriguingly, DEMO exhibits inferior IGD and Spacing values compared to NSGAII, yet achieves superior HV* metrics.

Examining Figure 7, it is evident that only IBEA demonstrates better spacing values than ILSDEMO on DTLZ2. For DTLZ1 and DTLZ2, IBEA outperforms NSGAII and DEMO, with the exception of inferior IGD values compared to DEMO on DTLZ2, where NSGAII performs the worst. DEMO exhibits superior performance over NSGAII and IBEA on DTLZ3 and DTLZ4. Specifically, NSGAII outperforms IBEA in terms of IGD on DTLZ3, but lags behind in spacing and HV* values. For DTLZ4, NSGAII demonstrates better IGD and Spacing values than IBEA, barring some outliers, but lags in HV* values.

Figure 8 indicates that, for DTLZ5, IBEA performs the worst, while NSGAII outperforms DEMO. For DTLZ6, NSGAII performs the worst, except for its superior spacing values compared to IBEA. DEMO achieves similar performance to ILSDEMO, albeit with inferior spacing values.

5. Conclusions

The quality of non-dominated solutions is influenced by several factors, including the initial population, variation operators, and selection operator. This paper presents an efficient multi-objective evolutionary algorithm based on differential evolution. Multiple approaches from different perspectives are employed to enhance the algorithm's performance. Orthogonal design is utilized to generate a statistically robust population. Additionally, self-adaptive parameters and strategies are incorporated to automatically adjust parameters, enabling adaptability to varying problems and stages. The q-Gaussian mutation is leveraged to strengthen the local search capability in the greedy selection of differential evolution. Furthermore, the k-nearest neighbor rule is used to prune crowded individuals, maintaining a fixed-size archive population. To ensure uniform distribution of individuals in the parent population, an indicator-based selection method is incorporated to form a new parent population based on the existing parent and archive populations. In comparisons with NSGAII, IBEA, and DEMO on ZDT-series and DTLZ-series test instances, the proposed approach demonstrates superior convergence to the true PF with greater accuracy and uniformity.