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

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-point-based 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 neighborhood-based 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) contains four main steps during the optimization process. The algorithmic details of the four steps are outlined in Algorithm 1 .

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:

Let $X=\left\{x_1,x_2,\dots ,x_N\right\}\in R^n$, x_i is a vector of X. Clustering is the process of dividing X into m non-empty subsets according to some criteria, which satisfies the following three conditions: (2)${\displaystyle C_i\ne \varnothing ,i=1,2,..,m}$ (3)${\displaystyle {\cup }_{i=1}^m{C}_i=X}$ (4)${\displaystyle C_i\cap Cj=\varnothing ,i\ne j,i,j=1,2,\dots ,m.}$

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={\left[a\left(i,k\right)\right]}_{N\times N}$ and matrix $R={\left[r\left(i,k\right)\right]}_{N\times N}$. Initialize the factor a(i, k) = 0. Then, update the two factors with the following formula: (5)${\displaystyle r\left(i,k\right)=s\left(i,k\right)-ma{x}_{k\ne k^{\prime }}\left[a\left(i,k^{\prime }\right)+s\left(i,k^{\prime }\right)\right]}$ (6)${\displaystyle a\left(i,k\right)=\left\{\begin{array}{c}min\left\{0,r\left(k,k\right)+\sum _{i^{\prime }\ne i}\left[0,r\left(i^{\prime },k\right)\right]\right\},i\ne k\hfill \\ \sum _{i^{\prime }\ne i}max\left[0,r\left(i^{\prime },k\right)\right],i=k\hfill \end{array}\right.}$

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 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₁, 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₁, where $p_1=1-\frac{G_c-1}{maxgen}$, G_c refers to the current generation number while maxgen refers to the maximum generation number. These two parameters G_c and $\frac{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₁ 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₁ 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. Next, repair the candidate solution 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^{\prime }$. Here, we adopted one more convenient way to replace parameter values that violate boundary 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.

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 non-dominated 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: (7)${\displaystyle A_{dis}=\varepsilon \ast {\left(\prod _{i=1}^n\left(x_i^{max}-x_i^{min}\right)\right)}^{\frac{1}{n}}}$

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. Where n denotes the dimension of the DS. $x_i^{max}$ and $x_i^{min}$ 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: (8)${\displaystyle LocalM=\sum _{j=1}^{n_i}{B}_{i,j}}$ where B_i,j = 1 if 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 ( $\left|A_{t+1}\right|\le 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.

Experimental settings

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, 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 normal 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 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 distribution 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 $\left\{0.01,0.02,0.05,0.1,0.2,0.5,1,2,5,10,inf\right\}$. 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₁ 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₁ selected a linear decreasing function. Figure 6B denotes the parameter p₁ selected a constant 1, which presents the DE/rand/2 strategy is adopted only. The Fig. 6C denotes the parameter p₁ selected a constant 0.5. Figure 6D denotes the parameter p₁ 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.