An Improved Multiobjective Particle Swarm Optimization Based on Culture Algorithms

Abstract

In this paper, we propose a new approach to raise the performance of multiobjective particle swam optimization. The personal guide and global guide are updated using three kinds of knowledge extracted from the population based on cultural algorithms. An epsilon domination criterion has been employed to enhance the convergence and diversity of the approximate Pareto front. Moreover, a simple polynomial mutation operator has been applied to both the population and the non-dominated archive. Experiments on two series of bench test suites have shown the effectiveness of the proposed approach. A comparison with several other algorithms that are considered good representatives of particle swarm optimization solutions has also been conducted, in order to verify the competitive performance of the proposed algorithm in solve multiobjective optimization problems.

1. Introduction

Since many real-life problems are in fact multiobjective optimization problems (MOPs), evolutionary algorithms used to solve MOPs have attracted much attention in recent years. Multiobjective optimization evolutionary algorithms (MOEAs) mainly include two branches: dominance-based and decomposition-based [1]. Nondominated sorting genetic algorithm II (NSGA-II) [2], strength Pareto evolutionary algorithm 2 (SPEA2) [3] and multiobjective particle swarm optimization (MOPSO) [4,5] are in the first branch; multiobjective evolutionary algorithm based on decomposition (MOEA/D) [6,7] is in the second branch. For those dominance-based algorithms, the main goals of multiobjective optimization are that the non-dominated set obtained by an algorithm is a good approximation to the true Pareto front, and that this non-dominated set is uniformly distributed on the Pareto front [8].

Particle swarm optimization (PSO) is an evolutionary computation technique based on the social behavior of species, such as a flock of birds or a school of fish.

To handle MOPs, two problems need to be taken into account. The determination of personal and global guides is the first, and second is to maintain the diversity of the population.

Some efforts have been made in selecting the global and local guides. Carlos introduced the secondary repository that is based on Pareto optimality to store the non-dominated particles in MOPSO, and the global guide is selected using the adaptive grid [4]. In crowding distance (CD)-MOPSO, the global guide is selected using crowding distance in MOPSO [9]. In clustering MOPSO, the global guide is selected based on the clustering technique [10]. Mostaghim and Teich [11] selected the local guide based on the sigma method, while Branke and Mostaghim [12] compared the results of different methods in selecting the local guide.

Coello proposed the use of a cultural algorithm to deal with multiobjective optimization problems, in which the knowledge based on the cultural algorithm will affect the selection of both the global guide and the local guide [13]. Cultural multiobjective quantum particle swarm optimization (MOQPSO) [1], Cultural MOPSO (CMOPSO) [14], constrained CMOPSO [15], and Cultural dynamic PSO (DPSO ) [16] are all based on a cultural evolution mechanism.

A special domination principle and genetic operator have been used to raise the diversity of the Pareto front. Laumanns defined ε-box dominance to combine convergence and diversity [17], while Zitzler introduced the elitism mechanism that performs crossover and mutation on individuals selected from the union of population and repository [18]. Simulated Binary Crossover (SBX) has been used in the work of Tang [19] and Yu [20]. The Borg MOEA represents a class of algorithms whose operators are adaptively selected based on the problem [21]. Diversity-loss-based selection method is used in [22]. Bare-bones (BB)-MOPSO use an approach based on particle diversity to update the global particle leaders [23]. Pareto entropy MOPSO (peMOPSO) using cell distance based individual density to select the global guide [24].

Although there has been considerable research conducted on PSO proposed to solve MOPs, to raise the convergence and diversity of the approximate Pareto front further still remains an issue that needs to be considered.

An improved multiobjective particle swam optimization has been developed in this paper. Situational, normative and topographical knowledge are to be used to update the personal and global guides. A polynomial mutation operator is applied to every particle in the swarm and in the repository to enhance the convergence and diversity of the Pareto front. Epsilon domination principle depends both on the situation in the adaptive grid and the objective value that has been used.

The paper is structured as follows. The concept of multiobjective problems, cultural algorithms, and particle swarm optimization are briefly reviewed in the next section. In Section 3, a new multiobjective particle swarm optimization is proposed to improve the convergence and spread to the true Pareto front. In Section 4, the results of the experiments and analysis are shown. Conclusions are given in Section 5.

2. Preliminary

2.1. Cultural Algorithms

A cultural algorithm usually includes five kinds of knowledge: situational, normative, topographical, domain, and historical knowledge [25].

Situational knowledge is a set of useful instances for the interpretation of the experiences of all individuals. Situational knowledge guides individuals to move toward exemplars.

Normative knowledge consists of a set of promising ranges of variables. It provides guidelines for individual adjustments. Normative knowledge leads individuals to plunge into a good range.

Topographical knowledge divides the whole functional landscape into cells representing different spatial characteristics, and each cell records the best individual in its specific ranges. Topographical knowledge leads individuals to catch up to the best cell.

Domain knowledge adopts information about the problem domain to lead the search. Domain knowledge is useful during the search process.

Historical knowledge keeps track of the history of the search process and records key events in the search. It might be either a considerable move in the search space or a discovery of a landscape change. Individuals use historical knowledge for guidance in selecting a direction in which to move.

Particles are evaluated by a performance function in each generation. An acceptance function determines which individuals will be used to update the belief space. Experiences of those chosen individuals are then added to the belief space by the update function. After that, knowledge from the belief space influences the selection of individuals for the next generation through the influence function. The framework of a cultural algorithm is shown in Figure 1.

2.2. Multiobjective Optimization Problems

Multiobjective problems (MOPs) are problems whose goal is to optimize multiple objective functions simultaneously. A MOP global minimum problem with n objectives is defined as:

min f(x) = [f₁(x), f₂(x), …, f_n(x)]
s.t.g(x) ≤ 0, h(x) = 0

(1)

where g(x) ≤ 0 and h(x) = 0 represent the constraints.

Definition 1.

(Pareto dominance) [26]: x is said to dominate y (denoted as x < y) if $\forall$i, f_i(x) ≤ f_i(y), $\exists$i, f_i(x) < f_i(y).

Definition 2.

(Pareto optimal set) [26]: P* = {x$|\neg \exists$y ∈ Ω, f(y) $\prec$ f(x)}.

Definition 3.

(Pareto front) [26]: PF* = {f(x) = [f₁(x), f₂(x), …, f_n(x)]|x∈P*}.

Definition 4.

(ε-Pareto dominance) [17]: x is said to ε-dominate y, denoted as x $\prec$ _εy, if f(x) · (1 + ε) $\prec$ f(y).

2.3. Particle Swam Optimization

The standard PSO [27] is:

$$\begin{gathered} {\mathrm{v}}_{id}^{(k+1)}=\mathsf{\omega }·{\mathrm{v}}_{id}^{\left(k\right)}+{\mathrm{c}}_1·{\mathrm{r}}_1({\mathrm{p}}_{id}^{(k)}-x_{id}^{(k)})+{\mathrm{c}}_2·{\mathrm{r}}_2({\mathrm{p}}_{gd}^{(k)}-x_{id}^{(k)}) \\ {x}_{id}^{(k+1)}=x_{id}^{(k)}+{\mathrm{v}}_{id}^{(k+1)} \end{gathered}$$

(2)

where v${}_{id}^{(k)}$ and x${}_{id}^{(k)}$ are the velocity and position of the d-th dimension of the i-th particle in the k-th iteration, respectively, ω is inertia weight, c₁ and c₂ are the acceleration coefficients. Parameters r₁ and r₂ are uniformly distributed random numbers within the range [0, 1]. The parameters p${}_{id}^{(k)}$ and p${}_{gd}^{(k)}$ indicate the best position that the i-th particle has found so far, and the best position that the global swarm has found so far. The procedure of PSO is shown in Algorithm 1.

Algorithm 1. The Procedure of Particle Swarm Optimization (PSO)

(1) Calculate the fitness function of each particle.

(2) Update p${}_{id}^{(k)}$ and p${}_{gd}^{(k)}$.

(3) Update the velocity and position of each particle.

(4) If the stop criterion is not met, go to (1), else p${}_g^{(k)}$ is the best position

3. The Proposed Algorithm

3.1. Three Kinds of Knowledge Based on Population

Situational knowledge is determined by all of the personal guide positions. Figure 2 gives the situational knowledge of a space with two objectives. The obtained local area is split into n parts identically in each dimension. Situational knowledge is used to generate the personal guide for each particle.

Normative knowledge is based on the minimum and maximum values of non-dominated particles. Normative knowledge of a space with two objectives is shown in Figure 3. The normative knowledge in Figure 3 is updated from L₁, U₁, L₂, U₂ to L₁’, U₁’, L₂’, U₂’.

L_i = min f_i(x), U_i = max f_i(x)

(3)

The normative space is later used to obtain the global guide of the MOPSO by establishing the topographical knowledge.

In order to construct the topographical knowledge, the normative knowledge is employed and then the space is divided into grids, in which each of the resultant hypercubes will then be characterized by all of the lower and upper limits of the hypercubes, as well as the number of non-dominated individuals that are located on that hypercube.

Figure 4 shows the topographical knowledge of a two-objective space. With each generation, the topographic knowledge will be updated. Topographical knowledge is used to find the global guide.

3.2. Personal Guide

As in cultural MOQPSO [1], one dimension of the position of one particle is substituted with a random number s generated in the band, decided by the personal guide of whole population. If the new one dominates the old one, then the new one is employed; otherwise, the old one continues to be used. The detailed procedure of a personal guide is shown in Algorithm 2.

Algorithm 2. The Detailed Procedure of A Personal Guide

For the d-th dimension of the i-th particle,

(1) bn+1 = max({pbestid|i = 1, 2, …, n});

b1 = min({pbestid|i = 1, 2, …, n});

bi = b1 + (bn+1 – b1)·(i – 1)/n, i = 2, 3, …, n

(2) si = rand(bi, bi+1);

(3) Replace x${}_i^d$ with si and get a new particle

(4) If the newly generated particle dominates xi, then it replaces xi; break;

if not, the newly generated particle is discarded

3.3. Global Guide

Making use of topographical knowledge, a hypercube is selected based on Roulette Wheel Selection, then one member of this hypercube is randomly selected to be the global guide. The probability of selecting one hypercube is inversely proportional to the exponential of the number of non-dominated particles in that hypercube, and the exponential parameter is used here to increase the probability of selecting the hypercube that has less non-dominated particles in it.

3.4. Dominate Criteria Based on the Hypercube

Inferred from [17], a new domination criterion is proposed. The location of each pair of members in the non-dominated archive is used together with the objective value to determine the dominance. The hypercubes that are dominated by hypercube 3, 2 are shown in gray in Figure 5.

3.5. Mutation Operator

A polynomial mutation is applied to each particle, and thereafter the non-dominated items produced are added into the repository. A polynomial mutation is applied to every member of the repository, the produced items are added into the repository, and then the repository based on dominates is updated. When the resulting offspring is non-dominated by the individuals in the repository, it is added into the repository. If any individual in the repository dominates the offspring, then the offspring is discarded.

The polynomial mutation is stated as [28]:

s_m = s + α·β_j

(4)

where α = $\{\begin{array}{c}{\left(2v\right)}^{1/\left(p+1\right)}-1,v<0.5\\ 1-{\left(2\left(1-v\right)\right)}^{1/\left(p+1\right)},\mathrm{otherwise}\end{array}$, β = u_j – l_j, p is a positive real number, v is a uniformly distributed random number in [0, 1], and β_j is the range of j-th dimension.

3.6. Updating the Repository

Inferred from [29], in which the archive is updated based on the convergence and diversity metric, the archive is updated based on crowding distance in this paper. When the number of non-dominated solutions exceeds the limit size of the repository, first the crowding distance of every item in the archive is computed, and then the one that has the minimum crowding distance is deleted. This is repeated until the number of non-dominated solutions is equal to the limit size of the repository. The update function of the dominance archive is shown in Algorithm 3.

Algorithm 3. The Update Function of the Dominance Archive

Suppose element f’ is a new element that is to be judged whether to enter the repository (rep).

(1) Input rep, f

(2) If ¬$\exists$f’$\in$ rep:hypercube(f)$\prec$hypercube(f’) then rep = rep∪{f}E

(3) If $\exists$f’:hypercube(f’) = hypercube(f)∩f$\prec$ε + f’ then rep = rep∪{f}f’

(4) If ¬$\exists$f’:hypercube(f’) = hypercube(f)∪hypercube(f’)$\prec$hypercube(f), then rep = rep∪{f}

(5) Output rep

4. Experiments

The proposed approach is compared to the following several multiobjective evolutionary algorithms to test its performance: MOPSO [4], cultural MOPSO [14], cultural quantum-behaved MOPSO [1], NSGA-II [2], and crowding distance-based MOPSO [9]. The implementation of some algorithms is available at http://delta.cs.cinvestav.mx/~ccoello/EMOO/EMOOsoftware.html.

ZDT and DTLZ test functions were taken to verify our approach [26]. The number of objectives in the ZDT test function is 2 and the number of objectives in the DTLZ test function is 3. The number of variables is 30 for ZDT1-3, 10 for ZDT4, 10 for ZDT6, 10 for DTLZ1-6 and 20 for DTLZ7.

The algorithm is run for a population size of 100, a repository size of 100 particles, and 30 divisions for the adaptive grid in each dimension of the objective space. The number of fitness function evaluations is 30,000. The results that obtained over 30 independent runs are reported in the following. The parameter of the PSO is time-varied as below [30,31]:

ω = 0.9 – (ita/max_ita) × (0.9 – 0.4)

(5)

c₁ = 2.5 – (ita/max_ita) × (2.5 – 0.5)

(6)

c₂ = 0.5 – (ita/max_ita) × (0.9 – 0.4)

(7)

4.1. Quantitative Performance Evaluations

Hypervolume is adopted as a quantitative measure in this article, as it is a measure of the performance in both convergence and diversity [32]. Hypervolume is defined mathematically as follows:

$$\mathrm{Hv}=^(\underset{p\in P{F}_i}{\cup }\{x|\mathrm{p}<x<\mathrm{r}\})$$

(8)

where ^ is the Lebesgue measure, and r ∊ Ri is a reference vector that is dominated by all of the valid candidate solutions in PFi.

Table 1. Hypervolume values for 30 runs of the test functions.

Test Functions		MOPSO	IC MOPSO	Cultural MOQPSO	Cultural MOPSO	NSGA-II	CD MOPSO
ZDT1	mean	5.56 × 10−1	8.73 × 10−1	8.42 × 10−1	8.66 × 10−1	8.11 × 10−1	8.74 × 10−1
	standard deviation (std)	9.03 × 10−2	5.38 × 10−3	8.01 × 10−2	7.00 × 10−3	8.46 × 10−3	4.76 × 10−3
	p-value	8.42 × 10−2		3.72 × 10−2	8.82 × 10−5	6.65 × 10−40	6.25 × 10−1
ZDT2	mean	2.88 × 10−1	5.38 × 10−1	4.72 × 10−1	4.20 × 10−1	1.20 × 10−1	5.37 × 10−1
	std	1.19 × 10−1	7.88 × 10−3	1.07 × 10−1	1.91 × 10−1	2.22 × 10−2	6.43 × 10−3
	p-value	1.69 × 10−16		1.40 × 10−3	1.31 × 10−5	6.53 × 10−66	8.25 × 10−1
ZDT3	mean	1.00	1.01	9.39 × 10−1	7.61 × 10−1	7.60 × 10−1	1.01
	std	5.17 × 10−3	4.01 × 10-3	1.43 × 10−1	4.00 × 10−1	1.31 × 10−2	4.37 × 10−3
	p-value	1.72 × 10−7		0.0110	1.35 × 10−3	2.35 × 10−66	8.76 × 10−1
ZDT4	mean	0.00	8.71 × 10−1	8.67 × 10-1	0.00	3.66 × 10−1	0.00
	std	0.00	6.08 × 10−3	5.18 × 10−3	0.00	1.77 × 10−1	0.00
	p-value	1.80 × 10−118		1.59 × 10−2	1.80 × 10−118	1.80 × 10−118	1.80 × 10−118
ZDT6	mean	4.94 × 10−1	5.04 × 10−1	4.65 × 10−1	5.03 × 10−1	0.00	5.05 × 10−1
	std	7.43 × 10−3	4.88 × 10−3	1.12 × 10−1	7.63 × 10−3	0.00	5.77 × 10−3
	p-value	5.97 × 10−08		6.31 × 10−2	8.61 × 10−1	3.40 × 10−110	8.76 × 10−1
DTLZ1	mean	0.00	1.30	1.30	0.00	0.00	0.00
	std	0.00	2.02 × 10−3	3.07 × 10−3	0.00	0.00	0.00
	p-value	2.90 × 10−156		6.44 × 10−1	2.90 × 10−156	2.90 × 10−156	2.90 × 10−156
DTLZ2	mean	6.74 × 10-1	7.15 × 10−1	4.27 × 10−1	2.46 × 10−1	6.93 × 10−1	7.20 × 10−1
	std	1.50 × 10−2	9.74 × 10−3	2.41 × 10−2	1.26 × 10−2	8.35 × 10−3	8.12 × 10−3
	p-value	4.76 × 10−18		3.83 × 10−54	1.33 × 10−78	3.17 × 10−13	3.85 × 10−2
DTLZ3	mean	0.00	7.16 × 10−1	7.15 × 10−1	0.00	0.00	0.00
	std	0.00	7.63 × 10−3	9.61 × 10−3	0.00	0.00	0.00
	p-value	8.20 × 10−108		9.47 × 10−1	8.20 × 10−108	8.20 × 10−108	8.20 × 10−108
DTLZ4	mean	6.82 × 10−1	7.29 × 10−1	5.46 × 10−1	1.67 × 10−1	6.65 × 10−1	7.22 × 10−1
	std	1.56 × 10−2	7.45 × 10−3	5.70 × 10−2	6.76 × 10−2	3.11 × 10−2	8.74 × 10−3
	p-value	1.55 × 10−21		9.53 × 10−25	6.11 × 10−47	9.22 × 10−16	8.66 × 10−4
DTLZ5	mean	4.34 × 10−1	4.39 × 10−1	4.35 × 10−1	2.43 × 10−1	4.38 × 10−1	4.38 × 10−1
	std	6.21 × 10−3	4.94×10−3	8.00 × 10−3	4.96 × 10−3	5.74 × 10−3	5.86 × 10−3
	p-value	1.65 × 10−3		3.72 × 10−2	2.13 × 10−77	5.64 × 10−1	3.14 × 10−1
DTLZ6	mean	4.29 × 10−1	4.40×10−1	4.21 × 10−1	2.36 × 10−1	0.00	4.39 × 10−1
	std	6.99 × 10−3	5.97×10−3	7.53 × 10−3	2.81 × 10−2	0.00	6.85 × 10−3
	p-value	1.43 × 10−8		5.64 × 10−15	3.58 × 10−43	9.89 × 10−102	6.28 × 10−1
DTLZ7	mean	1.27	1.43	1.21	5.44 × 10−1	1.47 × 10−1	1.38
	std	2.93 × 10−2	3.20 × 10−2	4.10 × 10−2	8.36 × 10−2	4.43 × 10−2	4.25 × 10−2
	p-value	1.19 × 10−28		2.91 × 10−31	1.98 × 10−51	5.57 × 10−73	5.54 × 10−6
Better/similar/worse		12/0/0		9/3/0	11/1/0	11/1/0	5/6/1

shows the mean and standard deviation values of hypervolume for over 30 independent runs obtained by all of the six algorithms. For further analysis, the T test is performed to evaluate the significant difference between two independent samples, and the results are illustrated in Table 1, where the value of α is set to 0.05. The best result of one test function with different methods is shown in bold. It can be seen that the proposed improved cultural MOPSO (ICMOPSO) has achieved the best performance in hypervolume for nine of the 12 test functions, and it is superior to MOPSO for all 12 of the test functions, superior to Cultural MOPSO except in ZDT6, superior to NSGA-II except in DTLZ5, superior to CD-MOPSO except in five functions, and finally inferior to CD-MOPSO only in DTLZ2.

4.2. Comparison of the Results Using Hypercube-Based $\epsilon$-Domination and Normal Domination Hypervolume

To demonstrate the effectiveness of the hypercube-based ε-domination criteria introduced in this paper, DTLZ1 and DTLZ7 are chosen as examples.

Figure 6 and Figure 7 show the Pareto fronts obtained using two kinds of dominate criterions for DTLZ1 and DTLZ7, respectively. It can be observed that the approximate Pareto fronts found by the improved Cultural MOPSO using hypercube-based ε-domination are more well-distributed.

The adjacent members of the repository may be deleted based on the hypercube-based domination criterion, because they may locate in the same grid which is not expected in scenes like ZDT4, ZDT6 and DTLZ6. So, for ZDT1, ZDT2, ZDT3, ZDT4, ZDT6, DTLZ5, and DTLZ6, normal domination criterion should be used.

4.3. Comparison of the Result with and without the Mutation Operator

The comparison of the Pareto fronts of ZDT1 and ZDT3 with and without the mutation operator is shown in Figure 8 and Figure 9. It can be seen that using the mutation operator can result in a more well-distributed approximation of the Pareto front.

5. Conclusions

A new approach has been proposed in this paper to improve the performance of MOPSO. It updates the personal and global guides based on three kinds of knowledge that are extracted from a cultural algorithm. An epsilon domination criterion has been utilized. What is more, a polynomial mutation operator is applied in order to raise the performance of the Pareto front. Experiments on a series of ZDT and DTLZ test functions has been conducted to compare the proposed method with several state-of-the-art MOPSO algorithms. The results indicate that the proposed approach is a viable alternative, since its average performance is highly competitive. The performance of the proposed algorithm in dynamic optimization problems merits further study in future works.