An improved differential evolution algorithm for multi-modal multi-objective optimization

Multi-modal multi-objective problems (MMOPs) have gained much attention during the last decade. These problems have two or more global or local Pareto optimal sets (PSs), some of which map to the same Pareto front (PF). This article presents a new affinity propagation clustering (APC) method based on the Multi-modal multi-objective differential evolution (MMODE) algorithm, called MMODE_AP, for the suit of CEC’2020 benchmark functions. First, two adaptive mutation strategies are adopted to balance exploration and exploitation and improve the diversity in the evolution process. Then, the affinity propagation clustering method is adopted to define the crowding degree in decision space (DS) and objective space (OS). Meanwhile, the non-dominated sorting scheme incorporates a particular crowding distance to truncate the population during the environmental selection process, which can obtain well-distributed solutions in both DS and OS. Moreover, the local PF membership of the solution is defined, and a predefined parameter is introduced to maintain of the local PSs and solutions around the global PS. Finally, the proposed algorithm is implemented on the suit of CEC’2020 benchmark functions for comparison with some MMODE algorithms. According to the experimental study results, the proposed MMODE_AP algorithm has about 20 better performance results on benchmark functions compared to its competitors in terms of reciprocal of Pareto sets proximity (rPSP), inverted generational distances (IGD) in the decision (IGDX) and objective (IGDF). The proposed algorithm can efficiently achieve the two goals, i.e., the convergence to the true local and global Pareto fronts along with better distributed Pareto solutions on the Pareto fronts.


INTRODUCTION
Most of optimal decision issues in research and application can be summarized as multiobjective problem (MOP),which is generally formulated as follows (Rao, 1991): min F (x) = (f 1 (x),f 2 (x),...,f m (x)) T  (1) where x = (x 1 ,x 2 ,...,x n ) T ∈ is the decision vector of dimension n and ∈ R n is the feasible space.The image set, S = {F (x)|x ∈ }, is called OS.The feasible solutions x 1 is dominated by x 2 on the condition that f i (x 1 ) ≤ f i (x 2 ),∀i=1 ,2,...,m and ∃j ∈ {1,2,...,m},f j (x 1 ) < f j (x 2 ).Furthermore, if a certain solution is not Pareto dominated by any other solutions, it is called a non-dominated solution.The set of all non-dominated solutions constitutes a Pareto optimal set (PS) and its image in objective space (OS) forms the Pareto front.
In recent years, the study on multi-objective evolutionary algorithms (MOEAs) (Deb, 2001) is in the development stage.Many efficient algorithms have been proposed, such as non-dominated sorting based genetic algorithm (NSGA-II) (Deb et al., 2002), multiobjective evolutionary algorithm based on decomposition (MOEA/D) (Zhang & Li, 2008), ε-dominance based multi-objective evolutionary algorithm (ε-MOEA) (Deb, Mohan & Mishra, 2005), improving the Strength Pareto Evolutionary Algorithm for Multi-objective Optimization (SPEA2) (Zitzler, Laumanns & Thiele, 2001).However, these MOEAs focused on finding the Pareto Front in OS, ignoring the convergence and diversity in decision space (DS).So, these algorithms are inefficient in solving MMOPs, which have multiple PSs in DS mapping to the same Pareto front in OS as shown in Figs. 1 and 2. Figure 1 depicts an MMOP example, illustrating that the distance between two points is large in the DS while mapping to the same objective vector.
Many problems in the real world have exhibited multi-modal characteristics, e.g., feature selection problems (Yue et al., 2019), the design of space missions (Schutze, Vasile & Coello, 2011) and imaging problems (Sebag et al., 2005).Obtaining all optimal solutions is necessary for decision-makers (DMs) to understand the problem entirely.Furthermore, it can be affected by some constraints that transfer to another solution.Therefore, decisionmakers hope to provide equivalent and diverse solutions for making decisions when dealing with MOPs.
To address MMOPs, many researchers have proposed different multi-modal multiobjective optimization algorithms (Qu, Suganthan & Liang, 2012;Li et al., 2017;Shi et al., 2019;Maity, Sengupta & Sha, 2019;Li, 2010) and have shown good performance on benchmark problems.The most representative is the omni-optimizer proposed by Deb & Tiwari (2005), which extended NSGAII with several techniques, including Latin hypercube sampling-based population, restricted mating selection and alternative crowding distance in both DS and OS.In 2016, Liang, Yue & Qu (2016) proposed a modified NSGAII called DN_NSGAII which adopted the crowding distance (CD) mechanism for forming the mating pool to enhance the diversity in the DS.Additionally, the CD method was used to maintain multiple PSs in the DS.However, the CD method is only considered in the DS without considering the OS.To address the limitations of existing multi-modal multi-objective particle swarm optimization techniques, Yue, Qu & Liang (2018) proposed a MO_Ring_PSO_SCD method, which utilized ring-topology to construct a stable niching to enhance the diversity.Moreover, a comprehensive special crowding distance (SCD) is designed to describe the crowding degrees considering both DS and OS.Then, Liu, Yen & Gong (2018) adopted two archive and recombination strategies in TriMOEA-TA&R to solve MMOP.Zhang et al. (2019) designed a new algorithm, using the clustering method to search for PSs and maintain the diversity.The proposed algorithm enhance the balance between exploration and exploitation.
Recently, numerous multi-objective multi-modal evolutionary algorithms (MMEAs) have been introduced for solving MMOPs.In 2021, Liang et al. (2021a) combined the clustering method and elite selection mechanism in MMODE.Clustering-based special crowding distance (CSCD) method is designed to obtain a well-distributed population in DS and OS.Kahraman et al. (2022) proposed an improvement of MOAGDE (Duman, Akbel & Kahraman, 2021) called DSC-MOAGDE.In this method, three spaces were created and strategies using different reference spaces were designed, which balances the exploitationexploration successfully and avoids the issue of premature convergence in both spaces.The DSC method has been proven to be an effective solution to overcome local solution traps.In 2022, Gao et al. (2022) proposed a decomposition-based multi-modal multi-objective evolutionary algorithm MOEA/D-SS, which adopts a density-based estimation strategy for estimating the number of PSs.MOEA/D-SS can faithfully locate all PSs more accurately.The environmental selection is developed to adjust sub-population size to maintain the diversity dynamically.To solve the problem of premature convergence, Keivanian & Chiong (2022) designed an enhanced ICA variant MOFAEICA to tackle MMOPs that

PRELIMINARIES OF STUDY
MODE is widely used for solving MMOPs, which are named as MMEAs.First, the definition of DE algorithm and APC is presented.The corresponding definitions, including local PS, global PS, local PF and global PF, can be referred to Li et al. (2022).In a typical MMOP, only one global PF maps to either a single PS or multiple global PSs, as shown in Fig. 1.However, in an MMOP with a local Pareto front (MMOPL), there exist local PFs and only one global PF, as depicted in Fig. 2.

Related studies
The literature shows that numerous researches have been carried out on (i) crowding distance method, (ii) mutation operator design, and (iii) environmental selection for better PF sets.
First, several crowding distance-based methods for selecting non-dominated solutions have been introduced in the literature.In NSGA-II, Deb et al. (2002) designed the crowding distance method to truncate the population and find PFs and PSs effectively.As mentioned in Kahraman et al. (2022), DS, OS and unified space are the main reference spaces in the crowding distance approach.Therefore, there are three cases to be discussed.
The first case only references the DS.For example, Liang, Yue & Qu (2016) proposed a modified NSGAII called DN_NSGAII, which employs the CD mechanism in DS.Wu a, Gong & Wang (2021) introduced a k-means and species-based algorithm (KSDE), which combines k-means clustering and species-based optimization techniques.In KSDE, clustering separates the population into sub-regions based on spatial positional relationships.The crowding factor in KSDE was equal to the population size.A few of these include evolutionary multi-objective optimization algorithm with a decomposition strategy in the decision space (EMO-DD) (Yang et al., 2021), a prediction strategy based on decision variable analysis (DVA) (Zheng et al., 2021), improved artificial electric field algorithm for multi-objective optimization (Petwal & Rani, 2020), and the evolutionary algorithm using a convergence penalized density method (CPDEA) (Liu et al., 2020).
The second case only references the OS.For example, Zadeh, Sayadi & Kosari (2019) presents the EMOPSO algorithm, an efficient meta-model-based collaborative optimization architecture for solving high fidelity multi-objective multidisciplinary design optimization problems.Star topology was employed in the Pareto-based crowding distance method.Aiming at the MOP with uncertain parameters, Huang et al. (2022) designed an outlier removal (OR) mechanism and a new crowding distance in the NSGA-II, the negative effects of parameter uncertainties are reduced.The new method, OR-NSGA-II, can adapt to the MOP with parameter uncertainties.A few of these include a reference-pointbased many-objective evolutionary algorithm (NSGA-III) (Deb & Himanshu, 2014), the multi-objective equilibrium optimizer with exploration exploitation dominance (MOEO-EED) strategy (Abdel-Basset et al., 2021), multi-constraint multi-objective evolutionary algorithm (MC/MOEA) (Liagkouras & Metaxiotis, 2021), and guided population archive whale optimization algorithm (GPAWOA) (Got, Moussaoui & Zouache, 2020).
The third case references both the OS and DS.For example, in the omni-optimizer algorithm (Deb & Tiwari, 2005), alternative crowding distances in both DS and OS are embedded to extend NSGAII.Yue, Qu & Liang (2018) modified the crowding distance in the omni-optimizer and proposed MO_Ring_PSO_SCD, which referenced both DS and OS adopting ring-topology to construct the stable niching for better explore ability.In this method, a comprehensive SCD is designed to describe the crowding degrees considering both DS and OS.Liang et al. (2021a) designed a CSCD method calculating the comprehensive crowding degree to achieve a evenly distributed population in both DS and OS.A few of these include a novel MMOEA with dual clustering in the decision and objective spaces (MMOEA/DC) (Lin et al., 2021), the Multi-modal Multi-objective Differential Evolution Algorithm Using Improved Crowding Distance (MMODE_ICD) (Yue et al., 2021), an evolutionary algorithm using a convergence-penalized density method (CPDEA) (Liu et al., 2020) and a novel multi-modal multi-objective differential evolution algorithm (MMO_CDE) (Zhou et al.,216).
Since the crowding-distance calculation is performed in the unified space, both the DS and OS simultaneously affect the density.Recently, Kahraman et al. (2022) proposed a crowding distance-based archive handling approach using unified space.
Second, studies are performed to balance exploitation and exploration capabilities by designing mutation operators.Yue et al. (2021) utilized a differential vector generation approach during the variation step to balance exploration and exploitation, thereby enhancing the diversity in both DS and OS.Javadi, Zille & Mostaghim (2021) employed weighted sum crowding distance and a neighborhood mutation operator to improve the convergence performance of the NSGA-II.Zhou et al. (216) proposed a neighborhoodbased dual mutation operator by developing two novel mutation operators, which provide a better balance between exploration and exploitation in locating multiple equivalent PSs.Liang et al. (2021b) designed a mutation strategy based on variable classification and elite individuals to improve the performance of offspring by using two different mutation strategies.To avoid premature convergence, Petwal & Rani (2020) combined bounded exponential crossover with polynomial mutation.Abdel-Basset et al. (2021) adopted the Gaussian method to improve mutation strategy.
Third, numerous recent works have been conducted on the environmental selection process.Wei et al. (2022) designed the environmental selection strategy by combining non-dominated sorting and hierarchical clustering to choose promising solutions for the next generation.Zhou et al. ( 216) developed a neighborhood-based dual mutation (NDM) strategy and a clustering-based environmental selection (CES) mechanism to maintain the convergence and diversity of Pareto optimal solutions.Gao et al. (2022) combined the estimation strategy and the greedy selection in the environmental selection phase to maintain the population diversity.Das et al. (2023) employed an improved environmental selection strategy to pick out high-quality solutions and balance convergence and diversity in solving MOPs.Li et al. (2021) designed an environmental selection strategy that designed an operator that removes inefficient solutions to maintain good convergence and diversity of PSs.Zhang et al. (2022) proposed a two-stage environmental selection strategy to obtain better convergence of the OS and the distribution of the DS.In addition, several successful studies on environmental selection mechanisms are a novel constrained many-objective optimization EA with enhanced mating and environmental selections (CMME) Ming et al. (2022), a dynamic multi-objective optimal evolutionary algorithm which is based on environmental selection and transfer learning (DMOEA-ESTL) He, Zheng & Peng (2022), an MOEA with a flexible reference point for clustering (MOEA/FC) Liu et al. (2021), and a Multi-modal Multi-Objective differential evolution algorithm using the clustering-based special crowding distance method and elite selection mechanism (MMODE_CSCD) (Liang et al., 2021a).
To summarize, to identify global and local PSs that meet user's preferences in resolving MMOPs, it is essential to carry out intensive research on the crowding-distance calculation, mutation operator design and environmental selection in MMEAs.However, compared to recent studies, there is no universally recognized approach for the above three procedures.We are confident that the approach introduced in this study can address the current gap effectively.

Differential evolution
Differential evolution (DE) is a simple and powerful evolutionary algorithm introduced by Price & Storn (1997).DE uses a series of operations to generate offspring, such as mutation and crossover of current population individuals.DE algorithm also has lower spatial complexity.Because of its excellent characteristics, the DE algorithm is more conducive to the processing of large-scale, high-complexity optimization problems in many fields, including computer science, biology, and industry (Cappiello et al., ;Bano, Bashir & Younas, 2020;Ibrahim et al., 2015;Taha et al., 2023;Chen et al., 2022;Zhang, Yu & Wu, 2021).The classical DE algorithm (Storn & Price, 1997)

Affinity propagation
The clustering method (Marques & Wu, 2002;Frey & Dueck, ) divides objects into groups or clusters according to their similarities in some attributes.The general description of clustering analysis is as follows: Clustering is the process of dividing X into m non-empty subsets according to some criteria, which satisfies the following three conditions: The vectors contained in C i are more similar to each other.The C i are called clusters of X .The typical clustering process includes the following steps: feature selection, proximity measure, clustering criterion, clustering algorithm, and result validation.The similarity matrix is denoted by S N ×N , and the similarity is measured by the Euclidean distance between two data points.The elements on the diagonal line of the matrix S(i,i) are replaced by p(k),k = 1,2,...,N , which affects the final number of clusters.The larger value of p(k) is associated with the greater probability that this data point will eventually become a cluster.Therefore, P = [p(1),p(2),...,p(N )] is called the reference degree.In the affinity propagation algorithm, responsibility and availability are two significant factors to indicate information between points; the former is denoted as responsibility r(i,k), and the latter is availability a(i,k).The responsibility r(i,k) reflects the accumulated evidence for how suitable point x k is to serve as the exemplar for point x i , which is illustrated in Fig. 3A.The availability a(i,k) reflects the accumulated evidence for how appropriate it would be for point x i to choose point x k as its exemplar, as shown in Fig. 3B.In practice, the cyclic calculation process of AP algorithm is the process of exchanging massage of matrix A = [a(i,k)] N ×N and matrix R = [r(i,k)] N ×N .Initialize the factor a(i,k) = 0.Then, update the two factors with the following formula: (5) For the point x i , the value of k that maximizes a(i,k) + r(i,k) can identify point x i as an exemplar if k = i, otherwise point x i belongs to the cluster that k is the exemplar.The AP algorithm implements an iterative process searches for clusters until a high-quality set of exemplars and corresponding clusters are assembled.The procedure of the AP algorithm is depicted in the following Algorithm 2:

Framework of MMODE_AP
The procedure and flow chart of the MMODE_AP algorithm are depicted in Algorithm 3 and Fig. 4, respectively:

Solution generation
Since the differential evolution operator (Price, Storn & Lampinen, 2005) usually outperforms other mutation operators in single-objective optimization, the MMODE Algorithm 2 Outline of AP's main procedure Input: DATA: points to be clustered, G max : maximum iterative generations Output: cluster which each data point belongs Step 1. Initialization: generation counter t = 0; reference degree p(k); similarity matrix S; Step 2. main loop: 2.1 Calculate matrix A = [a(i,k)] N×N and R = [r(i,k)] N×N 2.2 Update the matrix A and R according to the formula: 5. Non-dominated sorting on P t to obtain P * 6.For each x i ∈ P t ,i = 1,2,...,n identify an exemplar 7. Generate a new solution offspring OP t = solution generation ( x i ) // Algorithm 4 8. Preserve the new solution P = OP t ∪ P * 9. Update the population P t +1 environmental selection ( P ) // Algorithm 5 10.Update the archive A t +1 // Algorithm 6 11.End while algorithm adopts it and uses the following two methods to improve its performance, which is shown in Algorithm 4.
Firstly, adopt the DE/rand/2 strategy or the DE/current-to-exemplar/1 strategy according to the random value.If the random value is less than p 1 , then select five random solutions from the whole population and utilize the DE/rand/2 strategy to produce a candidate solution.Otherwise, use the DE/current-to-exemplar/1 strategy with the probability 1−p 1 , where maxgen , G c refers to the current generation number while maxgen refers to the maximum generation number.These two parameters G c and G c −1 maxgen serve the same role of temperature and cooling schedule came from the method of simulated annealing (SA).The original DE does not have such a parameter and in many cases it may be considered an advantage and thus introduction of such a parameter is to be justified together with the formula for p 1 as other decreasing functions can be used here similar to SA cooling schedules.Hence, the exploration ability decreases from the initial process of evolution to the end as p 1 decreases from 1 to 0. In the proposed selection mechanism, individuals to generate difference vectors are selected in two methods adaptively.The first selection method enhances exploration ability and avoids iterations falling into local optima, the second method helps improve the convergence and diversity in DS and OS.Therefore, a probabilistic-based mutation strategy is used for the mutation.
Else Create a candidate solution use DE/current-to-exemplar/1 strategy: where a i and b i i = 1(,2,...,n) is the lower and upper boundaries.

4.
Crossover the candidate solution: Return y i .constraints with their upper and lower boundary values.Then, the crossover is used to mutate the candidate to improve its quality further.In this operator, C r denotes the control parameter.Sometimes, the newly produced candidate may violate the boundary constraints.In such cases, we simply replace them with the closest boundary value.Moreover, this approach does not require the construction of a new candidate.It is worth mentioning that the DE's mutation operation is only included in the MMODE_AP.That is, the selection operation is redefined in the next part.This design's purpose is that the DE mutation operator is invariant of any orthogonal coordinate rotation, which can solve multi-complicated Pareto sets (PSs) (Li & Zhang, 2009).The reason why these two mutation strategies are adopted in the solution generation process is that the exploration and exploitation abilities are balanced, and the diversities in DS and OS are improved.In the first strategy, the DE/rand/2 approach is used to enhance the exploration ability and avoid trapping into the local area.In the second strategy, select the one in the first non-dominated front as the exemplar.x exemplar refers to the learning model selected by non-dominated sorting methods, intending to achieve a uniform distribution of solutions on PS.Learning models represent individuals with excellent quality and low density, and we hope to generate more excellent solutions around these model individuals.'Comprehensive comparison' demonstrates the efficacy of the suggested adaptive individual selection approach.

Environmental selection
The environmental selection aims to truncate the population to obtain the best NP individuals.In MMODE_AP, the truncation, which consists of sorting the individuals with non-dominated sorting and evaluating the individuals with the crowding distance metric proposed in NSGA-II, is employed to create a promising population.Algorithm 5 shows the procedure of environmental selection.

Algorithm 5 Outline of environmental selection
Input: current population P Output: next generation population P t +1 1. Assign the solutions x ∈ P to fronts F 1 ,F 2 ,...,F l and calculate the CSCD using APC method// Algorithm 2 2. Add all individuals in P to the next parent population P t +1 is constructed from the member of the sets F 1 ,...,F l−1 , and from N −s members of the set F l according to the CSCD. 5. End if 6.Return P t +1 .

Archive update
After the environmental selection strategy, the diversity of the global optimal solution is maintained.However, the obtained solutions are not satisfactory regarding their convergence quality.Furthermore, all the local PFs are not completely reversed.The archive, including Archive_global and Archive_local, is applied to improve the convergence performance and control the quality of local PFs respectively.Notably, a predefined parameter is introduced to control the completeness of the local PSs and PFs.
The steps followed by updating the archive are explained in Algorithm 6.The primary task is to preserve two populations, Archive_global and Archive_local.First, execute nondominated sorting on the current Archive A t and pick up the non-dominated solutions to form Archive_global. Then add the remaining solutions into RemainPop and delete the solutions in RemainPop close to solutions in Archive_global.It is generally considered that if a point is close to one point in the DS, they are also likely to be near in the OS.Calculate the distance of all pairs of solutions in the RemainPop.Specifically, the Euclidean distance is adopted to measure.Moreover, a predefined parameter ε is introduced to calculate the Pareto domination relationship between solutions.Specifically, we have the following formula: Here, A dis represents the average distance of the DS.In the same layer, x j is considered the neighbor of x i if and only if the distance of x j and x i is smaller than A dis in the DS.denotes the dimension of the DS.x max i and x min i denote the maximum and minimum values of the i_th decision variable of current front, respectively.ε is a domain positive parameter controlling the average number of neighbors.Specifically, the optimal solutions of local PSs can not be dominated by any other solutions in the neighbor solutions.Moreover, the local PF membership is calculated within neighboring solutions.The neighbors of x i are denoted as N x .If there is no solution in N x dominates x i , then x i is a local Pareto optimal solution of current layer of RemainPop.This way, the solution in the local PS will not be compared with the solution in the global PSs.Therefore, the local PSs can be found.The local PF membership of x i is computed as follows: where

Benchmark problems
In order to validate the performance of the proposed algorithm, this chapter adopts the CEC'2020 multi-modal multi-objective optimization benchmark suite to test the algorithm.This benchmark suite includes 24 test functions with multiple features, as shown in Table 1.Details about the test suite can be found in the corresponding literature.Consequently, the experimental conditions for MMEAs were used (Liang et al., 2020).Experimental run settings are provided below: • All algorithms run 21 times independently on each benchmark problem.
• The population size is set as 200 * N_ops and N_ops means the number of local and global PS to be obtained for the problem.
• Maximal fitness evaluations (MaxFEs) was used as the termination criterion.The MaxFEs value was assigned as 10,000 *N_ops.
The details for each benchmark problem can be referred to (Liang et al., 2020).

Performance metrics
Four performance indicators including rPSP (Yue et al., 2019), 1/HV (Liang et al., 2020), IGD (Zhou, Zhang & Jin, 2009) in DS (IGDX) and IGD in OS (IGDF) are employed to discuss the results of the competing algorithms by Liang et al. (2020).Among these four indicators, rPSP and IGDX measure the extent of spread attained in the obtained optimal solution in DS, while rHV and IGDF are adopted to examine the population distribution in OS.Furthermore, a satisfactory IGDX generally has a beneficial IGDF, while a reasonable IGDF does not mean a favorable IGDX (Liang et al., 2020).

Competitor algorithms
In order to validate the effectiveness of the proposed algorithm, the state-of-the-art MMEAs, DSC-MOAGDE (Kahraman et al., 2022), MO_PSO_MM (Liang et al., 2018), MO_Ring_PSO_SCD (Yue, Qu & Liang, 2018), DN-NSGAII (Liang, Yue & Qu, 2016), Omni_Opt (Deb & Tiwari, 2005), NSGA-II (Deb et al., 2002), and SPEA2 (Zitzler, Laumanns & Thiele, 2001) are chosen as the comparison algorithms.In order to get statistically sound conclusions in this comparison, the specific parameter settings in each algorithm are consistent with the original papers and the actual reference points for both PFs and PSs are provided in the original papers.Moreover, the Wilcoxon rank-sum test (Wilcoxon, 1945) with p < 0.05 and the Friedman test (Friedman, 1939) are employed to compare MMODE_AP with other competitor algorithms.If MMODE_AP performs better or worse than the competitor algorithm, '' + or ''− is labeled respectively.In addition, if the two algorithms have no obvious difference, ''≈ is labeled.

Comprehensive comparison
The performance of the proposed algorithm on test suites is analyzed in this section.The statistical metric values (i.e., mean and standard deviation values) are shown in Tables 2-5.
In comparison tables, the best index values are highlighted in bold.In addition, the average rankings calculated by the Friedman test are shown in Table 6.Due to article length limitations, the results of a typical approximated set obtained by all algorithms are provided in the Fig. S1.The superiority of MMODE_AP is depicted in the statistical data in Tables 2-6.Statistical comparison results show that MMODE_AP has the best solution performance for the test problems used in the experiment.As we know, the smaller the value of rPSP and 1/HV, the better the convergence and diversity of obtained solutions.Table 2 shows that NSGAII achieves the best average 1/HV values on ten functions.DSC-MOAGDE and DN-NSGAII rank second and third on 9 and 8 functions, respectively.MMODE_AP obtains good results, with a smaller 1/HV than other competitors, but occasionally worse than DN-NSGAII, DSC-MOAGDE and NSGA-II.For MMF7, MMF11, MMF12 and MMF13 test functions,

Notes.
Bold numbers indicate the best result in the row of data.
MMODE_AP yield equal mean 1/HV values compared with the winner algorithms.
Based on archive update and environmental selection mechanism, MMODE_AP achieves an excellent approximation to the true PF on most of these MMOPs.However, it is difficult to obtain a better approximation to the true PS and true PF simultaneously.That is, the number of the best mean 1/HV values obtained by MMODE_AP is less than other algorithms.The reason might be that MMODE_AP is focused on solving MMOP with the local Pareto front.Table 3 shows that MMODE_AP achieves competitive performance compared to its competitors regarding rPSP.Especially, its performance on functions with multiple PFs is outstanding.Accordingly, the results show that the performance of MMODE_AP whose framework embedding the AP approach is excellent in terms of convergence.Compared with DN-NSGAII and NSGA-II, MMODE_AP has a narrow advantage in some instances.Table 4 shows that MMODE_AP obtains better results, with the smaller IGDX compared with some state-of-the-art algorithms on the test problems.In particular, MMODE_AP wins 12 instances over 24 test suites.From the tables, MMODE_AP and MO_PSO_MM algorithms have demonstrated significant superiority over other existing algorithms, whose highlighted mean value is the most.MMODE_AP obtained one order of magnitude better than the traditional migration algorithm in terms of IGDF in solving most MMOPLs (MMOP with a local Pareto front).DN-NSGAII and Omni-optimizer performed poorly on any of the 24 benchmark problems, although they are specially designed for MMOPs.Table 5 lists the mean and standard deviation IGDF results of the competitor algorithms.Intuitively, MMODE_AP wins 18 instances over 24 test problems.MMODE_AP has a great superior performance over competitors on the benchmarks.Note that for some instances, such as MMF5, MMF8, MMF10, MMF15 and MMF15_a, there is no substantial variation between the MMODE_AP and winner algorithms.Similar to that of IGDX, MMODE_AP obtained one order of magnitude better than the traditional migration algorithm in terms of IGDF in solving most MMOPLs.Tables 3-5 shows that MMODE_AP performs similarly to competitor algorithms for  6, which shows that MMODE_AP obtains the best result on rpsp, IGDX and IGDF indicators.In terms of 1/HV, NSGAII achieved a better ranking than MO_Ring_PSO_SCD and SPEA2.This finding makes the 1/HV metric controversial as studies have reported MO_Ring_PSO_SCD and SPEA2 to be more comparative than NSGAII (Shang et al., 2020;While et al., 2006).Moreover, the 1/PSP, IGDX, and IGDF metrics give more consistent results than the 1/HV metric.In summary, comprehensive results reveal that the best comprehensive result is achieved by MO_PSO_MM, followed by MMODE_AP and DSC-MOAGDE.This indicated that the top three algorithms performed competitively.Now, we discuss the performance of MMODE_AP.The most outstanding difference occurs with the results for MMF15_l on three objectives.Figs.S2 and S3 shows the  of the obtained solutions in the OS and DS on the MMF15_l benchmark, which has a global PF corresponding to one PS and a local PF corresponding to one PS.Notably, DSC-MOAGDE is considered the latest effective algorithm that can solve MMOP.As the figures illustrate, MMODE_AP shows significantly better performance compared with DSC-MOAGDE.Meanwhile, SPEA2 loses one PS and can find only a single non-uniform distribution PS.Table 7 shows the four indicators obtained by MMODE_AP on MMF15_l instances and the best values are highlighted in bold.From this table, we can observe that MMODE_AP performs best regarding rPSP, IGDX and IGDF indicators.However, other algorithms have worse rHV results than DN_NSGAII.It is reasonable to conclude that MMODE_AP locates all PSs with good convergence and overlap ratio, unlike most competitor algorithms that only detect one PS.
Due to article length limitations, the final approximated set for two (three) objectives of the test suite of CEC'2020 are not presented.For more information, please refer to Fig. S1.The visual comparisons show that the optimal solutions yielded by MMODE_AP successfully approximate the whole PFs.Especially for the top ten benchmarks, we can observe that the objective points generated by MMODE_AP distribute alone the PFs well on two-or three-dimensional space from Fig. S1.To summarize, MMODE_AP shows the superiority of balancing the convergence and diversity of solutions.

Parameter analysis
In the Archive update strategy, ε is an important parameter that may affect the results of MMODE_AP.Improper parameter settings may lead to significant performance deterioration.In this section, the effect of the parameter on the performance of MMODE_AP is analyzed by performing relevant experiments on three representative functions (MMF5, MMF1_e, and MMF15_l).First, ε is empirically set as 0.01,0.02,0.05,0.1,0.2,0.5,1,2,5,10,inf .Figure 5 presents the average values calculated over the 3 test problems for each indicator (i.e., IGDX and IGDF).From Fig. 5A, we can observe that a good average IGDX is obtained when the parameter value is in [0.01,0.2].For the MMF1_e test problem, the results are clearly bad when ε is too large.In Fig. 5B, the IGDF value remains constant for the MMF5 and MMF15_l test problems as ε increases.For the MMF1_e test problem, the IGDF values are fluctuating.Based on the above discussion, the performance of MMODE_AP is favorable when ε ∈ [0.01,0.1].However, limitations also exist in this experiment.For example, the parameter range for obtaining satisfactory results is wide.Therefore, it increases the cost of time for decision-makers.

Analysis of the effectiveness of the proposed solution generation
To verify the effectiveness of the proposed selection mechanism, MMODE_AP with different p 1 are compared.The results on MMF3 are taken as an example and they are similar to the other test functions.The obtained PSs with median 1/PSP in the 21 times of these parameters are shown in Fig. 6.The Fig. 6A denotes the parameter p 1 selected a linear decreasing function.Figure 6B denotes the parameter p 1 selected a constant 1, which presents the DE/rand/2 strategy is adopted only.The Fig. 6C denotes the parameter p 1 selected a constant 0.5.Figure 6D denotes the parameter p 1 is a constant 0, which presents the DE/current-to-exemplar/1 strategy is adopted only.As shown in Figs.6B-6D, the population doesn't converge to the true PSs well.The obtained PSs of MMF3 are complete as shown in Fig. 6A.Overall, the the proposed selection mechanism is effective.

CONCLUSION
In order to enhance the exploration and exploitation of MMODE, this article proposes an improved MMODE named MMODE_AP to solve MMOPs.In MMODE_AP, the affinity propagation clustering method groups individuals in non-dominated layers into different classes, allowing individuals and their neighbors to be grouped in the same Pareto set.
Based on this method, the convergence of this algorithm can be accelerated in the process of optimization evolution.Then, two mutation strategies are adopted in the solution generation process.Based on the strategy, not only the exploration and exploitation performance are balanced but also the diversities in DS and OS are improved.Moreover, the CSCD is employed to create a promising generation in environmental selection.The comprehensive CD can reflect the actual crowding degree accurately.Additionally, an archive updating approach is designed to balance the convergence and diversity quality of solutions, in which a predefined parameter is designed to control the completeness of the local and local PFs.To demonstrate the potential benefit of the novel algorithm, a series of comprehensive experiments on CEC'2020 multi-modal multi-objective benchmark instances are conducted.Here is a summary of the pros and cons of the suggested approach: The proposed method has several advantages.Thanks to the implementation of the APC technique, the local and global PS (PF) can be well maintained.However, some limitations are also found in the experiment.For example, the performance of MMODE_AP in OS needs further improvements as it loses to competitors in terms of HV metric.Except for benchmark suites, developing potential multi-objective real-world test problems that contain multi-modal challenges for researchers is also essential.In future work, the implementation of the algorithm is to be discussed, including the potential for parallelization and optimization.Furthermore, more efficient maintenance mechanisms of local PSs need to be studied.

Algorithm 1 Outline of DE's main procedure Step
contains four main steps during the optimization process.The algorithmic details of the four steps are outlined in Algorithm 1. 1. Initialize population P randomly and evaluate the individuals of P.
Determine the center of the points 2.4 Distribution of each point 2.5 t = t + 1 Step 3. Stop iteration: if t > G max , print out the result , otherwise turn to Step 2.
Algorithm 3 Main steps of MMODE_AP Input: MOPs: Multi-objective optimization problems; n: population size ; G max : Maximum Generation; Output: If r < p 1 Select five different random solutions x r 1 ,x r 2 ,x r 3 ,x r 4 ,x r 5 from the whole population.Create a candidate solution use DE/rand/2 strategy: Kozlov, Samsonova & Samsonova (2016)to avoid violating the boundary constraints.There are other approaches for handling boundary constraints in DE.Lampinen (2002)replaced parameter values that violate boundary constraints with random values generated within the feasible range.Kozlov, Samsonova & Samsonova (2016)introduce a new parameter to transform y i .Here, we adopted one more convenient way to replace parameter values that violate boundary Algorithm 4 Outline of solution generation Input: parent solution x; adaptive method parameter p 1 ; Output: candidate solution y 1. Generate a rand value r 2.
and only if x j ≺ x i and B i,j = 0; otherwise.n i denotes the number of neighbors of x i in the set N x .Based on the definition of local PF, the one whose LocalM equals 0 can be considered a local Pareto optimal solution.An appropriate value of ε can balance the quality of the local PF (if it exists) and enhance the efficiency of finding all local Pareto solutions.It is worth noting that the setting of ε is dispensable on the prior information about the problem.Actually, ε can serve as an indicator of the difference between the global and local PFs.If the problem possesses several global and local PS, we suggest setting ε ∈ [0.01,0.1]according to the experiment results.Otherwise, ε = inf and MMODE_AP transform into a regular MMEA focusing only on global PS in this situation.By setting different values of parameter ε, the local PSs and the solutions surrounding the global PS can be maintained.Notably, if there are acceptable solutions in A t +1 ( |A t +1 | ≤ A size ), then we return A t +1 .Otherwise, if we prune Archive_local by non-domination sorting approach until the population size of A t +1 is less than the threshold value A size .The prune method is consistent with the environmental selection operation.

Table 5 Comparison of IGDF mean and standard deviation obtained by all algorithms for CEC 2020 benchmark problems.
MMOPs but significantly better for MMOPLs.In summary, MMODE_AP is competitive for MMOPLs and competent for normal MMOPs.Moreover, the Friedman test results are shown in Table normal