Multi-Objective Particle Swarm Optimization Based on Gaussian Sampling

This paper proposes a multi-objective particle swarm optimization algorithm based on Gaussian sampling (GS-MOPSO) to locate multiple Pareto optimal solutions for solving multimodal multi-objective problems. In the proposed method, the Gaussian sampling mechanism is used to form multiple neighborhoods by learning from optimal information of particles. And particles search their own neighborhoods to obtain more optimal solutions in the decision space. Moreover, an external archive maintenance strategy is proposed which allows the algorithm to maintain an archive containing better distribution and diversity of solutions. Meanwhile, nine new multimodal multi-objective test problems are designed to evaluate the performance of algorithms. The performance of GS-MOPSO is compared with twelve state-of-the-art multi-objective optimization algorithms on forty test problems. The experimental results show that the proposed algorithm is able to handle the multimodal multi-objective problems in terms of finding more and well-distributed Pareto solutions. In addition, the effectiveness of the proposed algorithm is further demonstrated in a real-world problem.


I. INTRODUCTION
Multi-objective optimization problems (MOPs) involve the optimization of multiple objective functions. In real-world applications, many optimization problems always consist of a several conflicting objectives and a number of constraints. Without loss of generality, a MOP can mathematically be described as follows: . . , f m (x)) T subject to g i (x) ≥ 0, i = 1, . . . , q h j (x) = 0, j = 1, . . . , p x ∈ R n (1) where x = (x 1 , . . . , x n ) is an n-dimensional decision vector in the decision space R n , F(x) = (f 1 (x), . . . , f m (x)) T is an m-dimensional objective vector in the objective space R m , g i (x) ≥ 0 is the i-th inequality constraint, h j (x) = 0 is the j-th equality constraint. Since the objectives in MOPs are often The associate editor coordinating the review of this manuscript and approving it for publication was Md. Asaduzzaman .
contradicting with each other, the Pareto dominance relationship is commonly utilized to compare different solutions. For any two feasible solutions x and x , it can be said that x dominates another solution x , if ∀i ∈ {1, 2, . . . , n} , f i (x) ≤ f i (x ) and ∃j ∈ {1, 2, . . . , n} , f i (x) < f i (x ). If the solution x dominates any other solutions, x can be said as a nondominated solution. The set of all non-dominated solution in the decision space is called as Pareto-optimal Set (PS) [1]. The set of the vectors in the objective space that correspond to the PS is referred to as Pareto Front (PF) [2], [3]. The crux of solving the MOPs is to identify evenly spread solutions in the objective space. A number of multi-objective evolutionary algorithms (MOEAs) have been proposed to tackle MOPs over the past two decades [4]- [9]. Until recently, the MOPs which have multiple disjoint PSs corresponding to the same PF attract great widespread interest in the evolutionary computation research community [10]- [15]. This class of problems is referred to as multimodal multi-objective problems (MMOPs) by Liang et al. [10]. In Fig. 1, the two PSs in the decision space map to the same PF in the objective space, in which two solutions A 1 and A 2 in the decision space correspond to the same objective value A in the objective space. If the solutions A 1 and A 2 are obtained by the traditional multiobjective optimization algorithm concurrently, the solution A 2 will be deleted since the solutions A 1 and A 2 are too crowded in the objective space. However, if the solution A 1 becomes unavailable due to unexpected reasons, solution A 2 could provide an alternative option with the same quality for decision-maker. Therefore, both A 1 and A 2 should be retained simultaneously. Hence, identifying and reserving multiple PSs in decision space is an important task for addressing the multimodal multi-objective problems [16], [17].
There exist multimodal property for both single-objective and multi-objective optimization problems. For multimodal single-objective optimization, plenty of research works has been done to locate all the global and local optima instead of a single global optimum. Various niching methods have been developed at the same time, including crowding [18], fitness sharing [19], clearing [20], and speciation [21]. However, the involving research in multimodal multi-objective optimization has been little investigated [11]. Most researches focus exclusively on finding the Pareto front in the objective space but does not consider locate multiple PSs in the decision space. Even though a small number of algorithms have been proposed in the literature for multimodal multiobjective problems, further improvements are called for to alleviate observed shortcomings. First, the ability to find enough equivalent optimal solutions in the decision space is still relatively weak. Second, relatively simple test functions are adopted in the experiment. Finally, an algorithm's capacity to tackle real-world problems is less studied. Based on the above discussion, the multimodal multi-objective optimization is worthy to be further studied. Meanwhile, an effective strategy to locate the multiple PSs in the decision space needs to be further developed.
In this paper, we propose a multi-objective particle swarm optimization algorithm based on Gaussian sampling (GS-MOPSO). In GS-MOPSO, a Gaussian sampling mechanism is adopted to establish different neighborhoods. And the particles evolve within their own neighborhoods that can locate more solutions in the decision space. Moreover, an external archive maintenance strategy is employed to maintain the diversity of solutions. Then, nine multimodal multi-objective test problems are designed to assess the performance of the algorithms. In addition, a real-world problem is solved by the proposed algorithm to further prove its effectiveness.
The rest of this paper is organized as follows. Section II reviews the related works. Section III introduces the details of the proposed GS-MOPSO. Section IV describes the experimental settings. And Section V reports the experimental results and the relevant analysis. Section VI summarizes the conclusions and describes the future work.

II. RELATED WORK A. PARTICLE SWARM OPTIMIZATION
Particle swarm optimizer (PSO), proposed by Kennedy and Eberhart in 1995 [22], mimics a flock of birds' foraging behavior to solve the optimization problems. Each particle can be viewed as a potential solution, which flies through the solution space [23]- [25]. The position and velocity of the particle are dynamically adjusted according to its historical personal best position (pbest) and the historical best position of its neighborhood (nbest). The velocity and position of the particle are updated according to (2) and (3), where v t (i), x t (i) denote the velocity and position of particle i of the tth generation, respectively. r 1 and r 2 are random numbers within the range [0,1]; w represents the inertia weight; C 1 and C 2 are the acceleration coefficients.
Owing to its simplicity, effectiveness, and reliability, PSO has been successfully used in many theoretic research and engineering applications. For practical purposes, various strategies are introduced into PSO to enhance the performance of the algorithm [26]- [32]. For instance, in order to avoid falling into local optimal solutions and improve the global searching ability, a stochastic inertia weight strategy is proposed [32]. In the proposed strategy, the inertia weight is adjusted by using the characteristics of random variables. The corresponding inertia weight update equation is as follows where σ is constant. µ min and µ max stand for the minimum and maximum value of the inertia weight, respectively. N (0,1) means random number of standard normal distribution.

B. PRIOR WORKS ON MULTIMODAL MULTI-OBJECTIVE OPTIMIZATION
Maintaining the distribution of solutions in the decision space has been done by a few researchers. Deb and Tiwari [33] introduced the concept of crowding distance into the decision space and proposed the Omni-optimizer to preserve the diversity of the solution in the decision space. Ulrich et al. [34] introduced decision space diversity into the hypervolume metric to increase the diversity in both decision space and objective space. Chan and Ray [35] adopted Lebesgue contribution and neighborhood count to enhance the diversity of population in the decision and objective space.
Rudolphet al. [36] employed a restart strategy to obtain a well-distributed Pareto set. Zhou et al. [37] proposed a probabilistic model based on MOEA to approximate the PS and the PF. Xia et al. [38] integrated crowding estimation method into MOEA, in which crowding distances were applied to guarantee diversity in the decision space. These works consider the distribution in the decision space, which produces positive effects on ameliorating the performance of the algorithm for handling multi-objective problems. However, the multiobjective problems where multiple PSs corresponding to the same PF is not well studied. Multimodal multi-objective optimization aims at locating multiple equivalent PSs. Therefore, niching techniques are incorporated into multi-objective algorithms by some researchers to promote the diversity of the population and improve the distribution of the solutions. Liang et al. [10] proposed a decision space based niching NSGAII (called DN-NSGAII), which employs the niching method to maintain the diversity in the decision space. Tanabe and Ishibuchi [39] proposed a decomposition-based evolutionary algorithm with addition and deletion operators to find more equivalent PSs. The results demonstrate that the proposed algorithm can maintain the population diversity in the decision space. Yue et al. [11] introduced a ring topology and special crowding distance into particle swarm optimization algorithm, in which a ring topology was deployed to form stable niches for maintaining multiple PSs, and a special crowding distance was utilized to balance the diversity in the both objective and decision space. Liang et al. [16] integrated self-organizing map network into particle swarm algorithm, in which the self-organizing map network is employed to establish the neighborhood of the individuals to identify a larger number of Pareto solutions. Liu et al. [15] introduced a multimodal multi-objective evolutionary algorithm using two-archive and recombination strategies to promote a good diversity solution in the decision space. Qu et al. [17] proposed a self-organized speciation based multi-objective particle swarm optimizer to locate multiple optimal solutions. The results show that the proposed method is competitive. Li et al. [40] combined reinforcement learning mechanism into differential evolution algorithm, in which the reinforcement learning mechanism is applied to increase the diversity of the solution set in the decision space.

C. MOTIVATION
There is more than one PS for multimodal multi-objective problems. For example, for MMF2 [11], the function has two PS, namely PS 1 and PS 2 , as shown in Fig. 2. In Fig. 2, A and B are particles, and C and D are candidate nbest solutions.
To locate more Pareto solutions, observed two problems are considered to be solved. The first one is how to guide the particles close to one of the PS. Particle A is far away from PS 1 and PS 2 , and has relatively large rank according to nondominated sorting approach compared with other solutions (i.e., B, C and D). Particle A and solution C are more likely to be in the same niche as they are close to each other compared with solution D. If solution C is selected as nbest to guide the update of particle A, it is conducive to converge and locate the solution effectively. In this way, the population can also form multiple niches that can converge to different solutions. The second one is how to search more potential solutions after locating one of the PS. In Fig. 2, particle B is close to PS 2 , which indicated that it in a potential area. And it has a lower rank in the population in accordance with nondominated sorting approach. If a local search is performed on particle B, it is helpful to obtain multiple desired solutions. Based on the above discussion, for the larger rank particle, candidate solution which is close to it is used to guide its movement and thereby close to the true PS. In contrast, for the lower rank particle, the particle evolves within its own niche by local search using the proposed Gaussian sampling mechanism. And the detail of the proposed multi-objective particle swarm optimization algorithm based on Gaussian sampling (GS-MOPSO) is described in Section III.

III. PROPOSED ALGORITHM
In the section, we present the main framework of the multiobjective particle swarm optimization algorithm based on Gaussian sampling (GS-MOPSO), which is composed of two components: Gaussian sampling mechanism, and external archive maintenance strategy. The details are described in the following subsections.

A. MAIN FRAMEWORK OF THE GS-MOPSO
Algorithm 1 outlines the main framework of the proposed GS-MOPSO, where the notation POP t (i) represent the ith particle at the tth generation, and PBA{i} denotes the ith particle's personal historical best positions. PBA retains each particle's historical best positions, which is advantageous for each particle to improve its position for the next generation by learning its own historical experience.
The procedure of GS-MOPSO is as follows. First, the particle population (POP), the personal best archive (PBA) and the external archive (EXA) are initialized. Then, each particle in POP is assigned a rank value by employing the fast non-dominated sort approach. The purpose of this step is to obtain the maximum rank value of the whole population (MaxRank) and the rank value of each particle (Rank) for Gaussian sampling mechanism. Next, the first solution from Algorithm 1 Framework of GS-MOPSO 1 Input: maximum number of iterations MaxIter, population size N. 2 Output: the non-dominated solutions in external archive EXA. 3 Initialize a random population POP 0 , and Evaluation POP 0 . 4 //Initialize personal best archive PBA and external archive EXA. 5 for i = 1: N 6 PBA{i} = POP 0 (i). 7 end for 8 EXA= POP 0 . 9 for t = 1: MaxIter 10 Assign rank value for each particle in POP t according to fast non-dominated sort approach. 11 for i = 1:N 12 //Select pbest and nbest. 13 pbest = The first one in sorted PBA{i}.
the sorted PBA{i} is chosen as the pbest. The solution closest to the ith particle in EXA is chosen as the nbest. Afterward, POP t (i) is updated to POP t+1 (i) according to (4), (2),and (3). Then, the Gaussian sampling mechanism is employed. After evaluation, POP t+1 (i) is stored into PBA{i} and all solutions dominated by POP t+1 (i) are removed. Finally, EXA is updated according to external archive maintenance strategy. Repeat the above steps until termination conditions are satisfied.

B. GAUSSIAN SAMPLING MECHANISM
Particle swarm optimization is easy to plunge into the local optimal, and appear premature stagnation phenomenon. And it has a slow convergence to the optimal solution and can not locate more than one solution without the niching technique. According to Section II.C, more potential solutions can be obtained by searching in the neighborhood of lower rank particle, so as to improve the disadvantages of the algorithm. Therefore Gaussian sampling mechanism is proposed. The Gaussian sampling mechanism includes global Gaussian sampling and local Gaussian sampling, which aims to Algorithm 2 Gaussian Sampling Mechanism 1 Requirement: rank value of the ith particle Rank i , maximum rank value maxRank, maximum number of iterations MaxIter, current iteration number t, personal best position pbest, neighborhood best position nbest, population POP, personal best archive PBA, external archive EXA.
, std(sampleSet)+ε) 7 else 8 //Local Gaussian sampling. 9 Divide the set [POP; PBA; EXA] into N clusters by using K -means clustering method. 10 Identify the cluster which the particle POP t (i) belongs, and assign it to the sampleSet. 11 POP t+1 (i) = N (mean(sampleSet), std(sampleSet)+ε). 12 end if 13 end if achieve a trade-off between global exploration and local exploitation. In the early search stage, global Gaussian sampling is employed to explore the search space comprehensively and to seek global optimal solutions quickly for promoting the exploration capability of the algorithm. The main idea is to guide the particle to move toward the global optimum by learning the optimal information of particles. In the latter stage of the search process, local Gaussian sampling is adopted to exploit the neighborhood of the promising solution and to locate more optimal solutions in the decision space for enhancing the exploitation capability of the algorithm. The main idea is to induce multiple neighborhoods in decision space and to guide particles to evolve in their own neighborhood.
The procedure of Gaussian sampling mechanism is shown in Algorithm 2, while rand is a random number uniformly distributed between 0 and 1, Rank i is the rank of the ith particle, t is current iteration, MaxIter is maximum number of iterations, N () is a Gaussian random function, mean() is the arithmetic mean, std() refers to the standard deviation, and ε is a small number to avoid std() equals zero. The process starts from judging rand < Rank i /maxRank to decide whether to adopt Gaussian sampling mechanism. If rand < Rank i /maxRank, the Gaussian sampling mechanism is adopted. Then if rand> t /MaxIter, the global Gaussian sampling is performed, otherwise the local Gaussian sampling is employed.
For global Gaussian sampling, the best positions found by population and individual are used to evolve particle.
By combining the best position found so far, the particle can be guided toward the global optimal direction effectively. The steps are as follows. Frist, a sampling set (sampleSet) is defined, and contains pbest, nbest. Then, according to the sample set, a new position is generated by Gaussian distribution. And the update formula is expressed as follows: As a result, the new position will be around the pbest and nbest, which increase the chance of finding potential solutions.
For local Gaussian sampling, particle's neighborhood information is adopted to generate its next position. The procedures are as follows. First, according to K -means clustering method [41]- [43], the set [POP; PBA; EXA] is divided into N clusters (neighborhoods). Then, the neighborhood to which the particle belongs is identified, and assigns it to the variable sampleSet. Next, the position of particle is updated using Gaussian distribution. And the update formula is given below: In this way, particles fine search in their own neighborhood. Therefore, it has the ability to quickly find approximate solutions and the search accuracy can be improved at the same time.
To summarize, the new position generated by particles is randomly distributed around the optimal solution in a Gaussian distribution to enhance the population diversity and avoid falling into a local optimal solution. Moreover, by searching the neighborhood of the optimal solution, more solutions in the decision space are located. More specifically, in the early stage of iteration, pbest is far from nbest, and the standard deviation will be large, which facilitates to explore in a wide range for particles. Then by continuously learning the pbest and nbest, the particles gradually approach and converge to the optimal solution. In the later stage of iteration, particles exploit their neighborhood according to the neighborhood information. Each particle evolves in its own neighborhood. Therefore particles can find the optimal solution in its neighborhood, so the whole population can locate more solutions in different neighborhoods. Besides, by using the information of PBA and EXA, the particles move toward the optimal solution promptly and can find more potential solutions. Thus the search efficiency and accuracy of the algorithm are improved.
The process for the Gaussian sampling mechanism can be summarized as follows: Step 1: Check whether the condition rand < Rank i /maxRank is satisfied, if so, go to Step 2; otherwise, exit the program.
Step 2.1: A variable sampleSet is defined and set to an empty value. Step 2.2: Assign [pbest; nbest] to sampleSet.
Step 2.3: Calculate the mean and standard deviation of sampleSet. Use formula (A5) to update the position of the particle POP t+1 (i). Then go to Step 4.
Step 3.1: Set variable sampleSet to an empty value.
Step 3.3: Identify the cluster to which the ith particle belongs, and allocate it to sampleSet.
Step 3.4: Compute the mean and standard deviation of sampleSet. Update particle's position POP t+1 (i) according to equation (A6). Then go to Step 4.

C. EXTERNAL ARCHIVE MAINTENANCE STRATEGY
The external archive is utilized to store the non-dominant solutions found so far and to guide the particles toward the true Pareto set [44], [45]. In general, a maximum size is defined for the external archive. When the size of solutions in the external archive reaches its predefined maximum value, the external archive is maintained to determine which solution can be retained into the EXA. If the diversity of the non-dominant solutions in EXA is poor, the particles in the population will gather in a certain region. Suppose that the maximum size of the external archive is six, and six solutions need to be remained in the external archive. The distribution of non-dominant solutions in one iteration is shown in Fig.3 as a simple example. In Fig. 3(a), there are eleven solutions, and the solutions A, B, and C belong to the same cluster and are quite close. Compared to adding all three solutions to the external archive (see Fig. 3(b)), maintaining only one of them into the external archive (see Fig. 3(c)) can greatly increase diversity with the same archive size. Based on the above idea, an external archive maintenance strategy is proposed to obtain the even distribution of non-dominant solutions in the decision space. External archive maintenance strategy at the ith generation is described as following.
First, copy all individuals from archive EXA and population POP t to a combined population R t . The population R t is of size 2N. Then, the population is sorted according to the non-dominated sort method. In this process, individuals are assigned to several fronts F 1 , F 2 ,. . . , F l according to the level of non-domination, where l is the value of the last non-dominated front. If the size of F 1 is less than N , the first N better order individuals are added to the EXA directly. Otherwise archive maintenance strategy will be performed. The standard K -mean clustering algorithm divides F 1 into N clusters. Then randomly select a solution from each cluster and add it to the EXA. In this way, the external archive EXA is updated.
In theory, the K -mean clustering method divides similar solutions into the same cluster, where these solutions are dissimilar from solutions in other clusters. Therefore, only one of the solutions in the same cluster is stored in the external archive, which not only avoids over-exploitation of a certain region by the particle in the evolution process, but also provides well-distributed solutions for the external archive.

IV. EXPERIMENTAL DESIGN A. TEST FUNCTIONS
The research on multimodal multi-objective optimization is still in the emerging phase, related benchmark functions from the literature are relatively few. Furthermore, many researchers are becoming interested in this area [11], it is essential to design complicated test problems for comparing the effectiveness of the multimodal multi-objective algorithms systematically. Thus, according to the approach of designing test problem [10], [11], nine multimodal multiobjective test problems (i.e., F1-F9) are proposed in this paper. Table 1 lists the relevant features of these nine test problems in this work, where n is the number of variables, and m is the number of objectives.
As is shown in Table 1, the function F1 has four linear PSs, and its PF is composed of discontinuous pieces. F2 has two PSs, and its PSs and PF are linear. F3 and F4 are originated from unimodal multi-objective problem UF3 [46], which have four PSs. F5 developed from UF2 [46] has two PSs, where one of the PS is an irregular geometry which is shown in Appendix.A. The building block of F6 is from UF7 [46]. The PSs of F6 are nonlinear, and its PF is linear. F7 and F8 are constructed based on UF4 [46]. F7 has convex PF, and F8 has concave PF. F8 has eight PSs, which is more complex than F7 with two PSs. The PSs shapes of F9 include the nonlinear and linear geometry simultaneously, which increases the difficulty degree of locating multiple solutions. Further details of the nine test problems are in the attached Appendix. In addition, thirty-one other test functions, called MMF1-8 [11], [47], MMF10-MMF15 [47], MMF1_e [47], MMF14_a [47], MMF15_a [47], MMF10_l-MMF15_l [47], MMF15_a_l [47], MMF16_l1-MMF16_l3 [47], SYM-PART simple [36], SYM-PART rotated [36], SYM−PART rot.+ trans. [36], map-based test problem (denoted as MBP) [48], and the Omni-test function [33], are included as test problems in the experiments.

B. PERFORMANCE INDICATORS
Pareto Sets Proximity(PSP) [11] and Inverted Generational Distance (IGD) [49], [50] are used to evaluate the performance of algorithm. PSP reflects the similarity between the true PSs and the obtained solutions, and IGD indicator can measure both convergence and diversity of the obtained solutions in the objective space. PSP is used to compare the performance of the algorithm in the decision space, while IGD is utilized to assess the performance of the algorithm in the objective space. The larger value of PSP means the obtained solutions in decision space are well distributed. And an algorithm with a larger PSP is considered to be a better design for multimodal multi-objective optimization. The small value of IGD means the convergence and diversity of obtained solutions in objective space are good.

D. PARAMETER SETTINGS
For fair comparison, the population size is set to 800 and the maximal evaluation is set to 80,000 according to the original study [11]. For each problem, 25 independent runs are carried out. In GS-MOPSO, both C 1 and C 2 are set to 2.05 [11], and µ min , µ max and σ are respectively set to 0.4, 0.9, 0.15 [56], and the parameter ε is set to 10 −2 . Discussion about the  parameter ε is covered in Section 5.4. The other parameter settings of the compared algorithms are set as their respective references [10], [11], [15]- [17], [33], [40], [51]- [55].

E. RUNNING ENVIRONMENT
All experiments are conducted independently on the same computer. The hardware conditions of the computer are as follows: CPU is Intel Xeon E5-2640 with 2.0GHz and 8GB main memory; the software platform is Windows 8.1 operating system; the simulation software is MATLAB R2016a. All the algorithms are coded and run using MATLAB.

V. RESULTS AND DISCUSSIONS A. EXPERIMENTAL VERIFICATION OF THE EFFECTIVENESS OF THE PROPOSED ALGORITHM
To illustrate the effectiveness of Gaussian sampling mechanism and of external archive maintenance strategy, the proposed GS-MOPSO is compared with the basic MOPSO [11], MOPSO-I (MOPSO only with the Gaussian sampling mechanism), and MOPSO-II (MOPSO only with the external archive maintenance strategy). The mean PSP values on forty test functions are presented in Table 2. Moreover, the Wilcoxon's rank sum test is applied to determine the statistical significance of the advantage of GS-MOPSO. ''+'', ''-'', and ''='' in Table 2 denote the performance of GS-MOPSO is better than, worse than, and equal to that of the other algorithm. Furthermore, the PSs obtained by these four algorithms on MMF3 are presented in Fig. 4.
It is observed that the PSP value obtained by MOPSO-I and MOPSO-II are marginally larger than that of MOPSO, while the PSP value achieved by GS-MOPSO is much greater than that of MOPSO. In addition, the rank sum test results show that there are significant differences among GS-MOPSO and the other three algorithms. Moreover, GS-MOPSO locates more Pareto solutions in the decision space than other three algorithms as shown in Fig. 4. The reasons are discussed in the following paragraphs.
Gaussian sampling mechanism enables the algorithm to locate enough optimal solutions. This mechanism can guide particles to search for the potential region by using optimal information. Meanwhile, it establishes multiple niches, and particles evolve independently in their own niche  to fine search for more solutions. The introducing of external archive maintenance strategy helps obtain a good distribution of optimal solutions in the decision space. The similar solutions are removed preferentially. In this way, if solutions are not similar (crowded) in decision space, they are able to survive and maintain in the archive. Therefore, the distribution of obtained solution is improved.
In conclusion, both Gaussian sampling mechanism and external archive maintenance strategy make GG-MOPSO more effective in solving multimodal multi-objective problems.

B. COMPARISON WITH OTHER ALGORITHMS
The mean values and standard deviations of the PSP metric for all the algorithms are presented in Table 3. The results show that the PSP values achieved by GS-MOPSO are highest on thirty-two test functions. In addition, according to the rank sum test, GS-MOPSO has better or similar performance as compared with the other compared algorithms on thirty-six test functions.
With respect to MMF3, the PSP value obtained by MMOPIO is a bit larger than GS-MOPSO. For MMF13_l, PEN-MOBA has a much larger PSP value than GS-MOPSO. For MMF15_a_l and F3, NMOHSA obtains the better PSP values than that of GS-MOPSO. However, the rank sum test results indicate that GS-MOPSO shows an equivalent performance with MMOPIO, PEN-MOBA, and NMOHSA. For MMF10, the SMPSO-MM is superior to GS-MOPSO. The reason may be that the SMPSO-MM employed the selforganizing map to extract neighboring relation information and help find the topological properties of the PSs, which makes the SMPSO-MM more suitable for solving MMF10. For MMF11_l and MMF12_l, GS-MOPSO is surpassed by ZS-MRPS. The key reason is that ZS-MRPS adopted the zoning search strategy to explore different subspaces and preserve the diversity of solutions for each subspace, which is more efficient for solving MMF11_l and MMF12_l. For F1, DE_RLRF obtains a higher PSP value. The reason is that DE_RLRF used the reinforcement learning to dynamically adjust the evolution direction of the population, which is more conducive to quickly find multiple PSs of F1. For the remaining thirty-two test functions, GS-MOPSO performs better than other algorithms. The main reason is that GS-MOPSO adopted the Gaussian sampling mechanism to learn from optimal information of the particles and establish the neighborhood of the particles, which is helpful in effectively finding more and better distributed solutions. In conclusion, GS-MOPSO obtains better performance for the PSP indicator comparing with twelve algorithms on the thirty-six test functions.
The IGD values are shown in Table 4 The IGD values of GS-MOPSO on the remaining twentythree test problems are the smallest among all the algorithms. In addition, the rank sum test results in Table 4 show that the GS-MOPSO performs statistically significantly better or equal when compared with other algorithms on majority of test functions. Therefore, it can be concluded that the GS-MOPSO obtains a better distribution in objective space.
A visual comparison is implemented to further demonstrate the performance of the algorithm. In this experiment, a complex function F8 is used, which has eight Pareto subsets in three-dimensional space. The obtained PS by different algorithms on F8 is shown in Fig. 5. Fig. 5 demonstrates that the proposed GS-MOPSO can find all the Pareto subsets. Meanwhile, GS-MOPSO achieves more uniform solutions than the other algorithms in each Pareto subset. In contrast, the Pareto solutions obtained by ZS-MRPS, DE_RLRF, SSMOPSO, MMOPIO, MMODE, PEN-MOBA, NMOHSA, SMPSO-MM, and MRPS are incomplete and uneven. Only parts of eight Pareto subsets are located by TriMOEATA&R, DN-NSGAII, and Omni-optimizer, meaning that the ability to identify more and well-distributed solutions in the decision space is relatively weak. The obtained PF on F8 is illustrated in Fig. 6. It can be observed that all of the algorithms obtain well-distributed PF on the F8 test problem. In conclusion, GS-MOPSO obtains the best distribution in decision space, which is an effective approach to deal with multimodal multiobjective optimization.

C. COMPARISON OF THE CONVERGENCE BEHAVIORS OF DIFFERENT ALGORITHMS
The convergence behavior of the algorithm is investigated on one typical test problem F6, which has two PSs in a three-dimensional decision space. The search space of F6 is divided into two sub-regions, namely Region 1{x 1 ∈ [−1, 0], There is one PS in each sub-region. The proportion of solutions in each sub-region for every generation is calculated to show the convergence behavior of the algorithm [11]. Ideally, if the algorithm performs well on F6, the proportion in each of the two regions should converge to 0.5. If the algorithm is trapped in the Region 1, the proportion in Region 1 will be greater than in Region 2.
The proportions for each region from the 1st to the 100th generation are shown in Fig. 7 (note that the ZS-MRPS uses the whole population to search one of the regions in one iteration, so the proportion of solutions in one region is 1 in an iteration, and the proportion of solutions in the remaining regions is 0. Hence it is excluded from this comparison study). Fig. 7 shows that GS-MOPSO maintains almost equal proportions in the two regions and the proportion in each  region is close to 0.5. This indicates that the GS-MOPSO algorithm converges to the two PSs of the F6 test problem. For DE_RLRF, the proportion of solutions in Region 1 is larger than that of Region 2 from the 40th to 100th generations. For SSMOPSO, the proportion of solutions in Region 2 is greater than that of Region 1 throughout generations. For MMOPIO, MMODE, SMPSO-MM, and MO_Ring_PSO_SCD, the proportion of Region 1 decreases to a small extent while the proportion in Region 2 increases slightly as the generation number reaches approximately 40. From the 40th to 100th generations, the proportion in Region 1 is smaller than those in Region 2, which indicates that most individuals converge to Region 2 in the evolution process. For TriMOEATA&R, DN-NSGAII, and Omni-optimizer, the proportion in Region 1 is much larger than 0.5 from the 20th to 100th generations, meaning that most individuals are more prefer to exploit Region 1. The convergence behavior of the PEN-MOBA is inferior to that of GS-MOPSO since from the 80th to 100th generation the proportion in Region 1 is greater than that in Region 2. For NMOHSA, the proportion in Region 2 is larger than that in Region 1 from the 60th to 100th generations. In addition, the proportion in each region fluctuates frequently throughout the iteration. In conclusion, GS-MOPSO has better convergence performance than other algorithms.

D. PARAMETER SENSITIVITY ANALYSIS
The parameter ε is designed to prevent the standard deviation from being a zero vector in Gaussian sampling mechanism. A large ε allows particles to explore a wide search space. However, this may lead to particles to deviate from VOLUME 8, 2020  the original sample area. A small ε would make particle to stagnate and even hardly escape from local optima, which is not conducive to locate more equivalent solutions. To investigate the sensitivity of the parameter ε, ε = {0, 10 −1 , 10 −2 , 10 −3 , . . . , 10 −15 } are tested on six test functions(MMF1-MMF6). Fig. 8 shows the distributions of the PSP values obtained by the proposed GS-MOPSO with different values of ε.
As reported in Fig.8, it can be observed that the parameter ε does affect the performance of algorithm. For MMF1, when the parameter ε is set to 10 −1 or 10 −2 , the algorithm achieves relatively high PSP value and exhibits a satisfying performance. When ε decreases from 10 −2 to 10 −4 , the PSP value gradually approaches the minimum PSP value. Then, the PSP value fluctuates slightly at the minimum PSP value as ε decreases from 10 −4 to 0. And a similar trend is also seen on MMF4, MMF5, and MMF6. For MMF2 and MMF3, the algorithm performance improves as ε decreases from 10 −1 to 10 −3 . When the value of ε to be 10 −3 GS-MOPSO obtains the highest PSP value. Then the PSP value fluctuates around the maximal PSP value as ε decreases from 10 −3 to 0. It is concluded that GS-MOPSO with ε ∈ {10 −1 , 10 −2 } has better performance on MMF1, MMF4, MMF5, and MMF6, while GS-MOPSO performs better on MMF2 and MMF3 when

E. APPLICATION
To further verify the algorithm's performance, this section considers its application to a real-world multimodal multiobjective problem. One of the well-known problems is mapbased test problem (MPB), which has three disconnected Pareto optimal subset [48]. This problem is solved with the proposed GS-MOPSO and the other twelve algorithms. And the results are provided in Table 3. From Table 3, the proposed GS-MOPSO obtains the highest PSP value than that of the other algorithms on the MBP. It can be concluded that GS-MOPSO shows a superior performance relative to the other algorithms. In addition, the Pareto solution sets obtained by each algorithm are displayed in Fig.9. It can be observed that the GS-MOPSO is able to cover all the three Pareto optimal regions. In contrast, for ZS-MRPS, MMO-PIO, MMODE, PEN-MOBA, NMOHSA, SMPSO-MM, and MO_Ring_PSO_SCD, most parts of the Pareto optimal sets are obtained. For DE_RLRF, SSMOPSO, TriMOEATA&R, DN-NSGAII, and Omni-optimizer, only a few solutions are found. In conclusion, the GS-MOPSO algorithm is effective in solving the practical multimodal multi-objective problem.

F. THE ROLE OF GAUSSIAN SAMPLING
The role of the Gaussian sampling mechanism is to quickly identify the global optimal solution according to the best information of particles in the early stage of iteration, and to locate more potential solutions by exploiting the neighborhood of the optimal solution in the later stage of iteration. In order to verify its effectiveness, the MOPSO with the Gaussian sampling mechanism (denoted as MOPSO-I) is compared with the basic MOPSO. The results are provided in Table 2. It can be seen that the performance of MOPSO-I is better than that of MOPSO on all test functions, which indicated that the Gaussian sampling mechanism improves the performance of the algorithm.
The Gaussian sampling mechanism consists of two parts: global Gaussian sampling and local Gaussian sampling. These two parts play different roles in the searching process, which are analyzed as follows.
The global Gaussian sampling is designed to increase the exploration capability of the algorithm, while local Gaussian sampling is designed to boost the exploitation capability of the algorithm. Therefore, the role of global Gaussian sampling is considered to explore the whole search space and to find the global optimum as quickly as possible. The role of local Gaussian sampling is considered to exploit the neighborhood of the optimal solution and to obtain more promising solutions as accurately as possible.
To further analyze the role of Gaussian sampling, MOPSO-I compares with MOPSO-I-v1 (MOPSO-I without the global Gaussian sampling) and MOPSO-I-v2 (MOPSO-I without the local Gaussian sampling) on forty test functions. Table 5 presents the mean and standard deviation values of PSP obtained by MOPSO-I, MOPSO-I-v1, and MOPSO-I-v2. It can be observed that the performance of MOPSO-I is significantly better than or equal to that of MOPSO-I-v1 on all the test functions, which indicates that global Gaussian sampling can help this algorithm to obtain diverse solutions and improve its exploration ability. MOPSO-I is significantly superior to or equal to  MOPSO-I-v2 on all the test functions, which reveals that the local Gaussian sampling can help this algorithm to converge on promising solutions and enhance its exploitation power.
To intuitively illustrate the role of Gaussian sampling, the SYM-PART simple test function as an example is chosen. This test function has nine Pareto sets. The final population distribution of different algorithms on SYM_PART simple is shown in Fig. 10. It can be observed from Fig. 10(a) that MOPSO-I is capable to locate all optima solutions. In Fig.10(b), the whole population is distributed around the optimal solutions, which denotes that the local Gaussian sampling makes the MOPSO-I-v2 has good convergence performance. In Fig.10 (c), the population distributed about the whole search space, which hints that the global Gaussian sampling makes the MOPSO-I-v2 has good diversity performance.
Overall, the comparison results illustrate that Gaussian sampling can help increase the performance of the algorithm.

G. LIMITATIONS AND FUTURE WORK
Though the proposed algorithm has good optimization performance on the majority of test functions, there are limitations of our research. On the one hand, the proposed algorithm does not perform well on test functions with discontinuous PS (e.g. MMF12_l, and F1) in Table 3. The reason is that the algorithm does not design relevant strategies according to the geometrical characteristics of PS. In fact, few strategies have been proposed by analyzing the characteristics of PS. Due to the PS geometry may be connected/disconnected, linear/nonlinear and other complex shapes, locating more solutions on complex PS is a very challenging task. One possible way to alleviate this difficulty is to adopt multiple strategies to guide the evolution of the population according to the different characteristics of PS. On the other hand, the parameter ε is set to a fixed value in the proposed algorithm. According to the analysis in Section V.D, the parameter ε affects the performance of the algorithm. In fact, the comparison algorithms (e.g. ZS-MRPS, SSMOPSO, and Tri-MOEATA&R) have similar situations. The reason is that it is difficult to determine the relationship between the parameter and the characteristics of PS. Therefore, future work should look deeply into this matter.

VI. CONCLUSION
In this paper, a multi-objective particle swarm optimization algorithm based on Gaussian sampling (GS-MOPSO) is proposed to solve multimodal multi-objective optimization problems. The Gaussian sampling mechanism adopts historical information of particles to divide multiple neighborhoods, and constantly search in the particle's neighborhood to locate more potential solutions in the decision space. In addition, an external archive maintenance strategy is utilized to improve the quality of the obtained solutions in the decision space. Moreover, a set of multimodal multi-objective test problems are proposed in this paper, i.e., F1−F9.
The proposed GS-MOPSO algorithm is compared with twelve multimodal multi-objective algorithms. All the algorithms are tested on forty multimodal multi-objective test functions. The experimental results show that GS-MOPSO has better performance than the compared algorithms in terms of PSP and is able to locate more and better distributed Pareto-optimal solutions. Finally, the proposed GS-MOPSO is applied to a real world problem, the map-based test problem. Results indicate that GS-MOPSO is feasible and effective.
In a future study, some of the real-world problems will be modeled to generate new multimodal multi-objective test functions with different characteristics. In addition, future work should further extend the problem set with multimodal many-objective or higher dimensional optimization. (A1) Its true PS is where0 ≤ x 1 ≤ 10.
Its true PS and PF are illustrated in Fig. 13.
Its true PF is where 0 ≤ f 1 ≤ 1. Its true PS and PF are illustrated in Fig. 14.
Its true PS is Its true PF is where 0 ≤ f 1 ≤ 1. Its true PS and PF are illustrated in Fig. 16 where −1 ≤ x 1 ≤ 1.
Its true PS and PF are illustrated in Fig. 18.