Optimizing Niche Center for Multimodal Optimization Problems

Many real-world optimization problems require searching for multiple optimal solutions simultaneously, which are called multimodal optimization problems (MMOPs). For MMOPs, the algorithm is required both to enlarge population diversity for locating more global optima and to enhance refine ability for increasing the accuracy of the obtained solutions. Thus, numerous niching techniques have been proposed to divide the population into different niches, and each niche is responsible for searching on one or more peaks. However, it is often a challenge to distinguish proper individuals as niche centers in existing niching approaches, which has become a key issue for efficiently solving MMOPs. In this article, the niche center distinguish (NCD) problem is treated as an optimization problem and an NCD-based differential evolution (NCD-DE) algorithm is proposed. In NCD-DE, the niches are formed by using an internal genetic algorithm (GA) to online solve the NCD optimization problem. In the internal GA, a fitness-entropy measurement objective function is designed to evaluate whether a group of niche centers (i.e., encoded by a chromosome in the internal GA) is promising. Moreover, to enhance the exploration and exploitation abilities of NCD-DE in solving the MMOPs, a niching and global cooperative mutation strategy that uses both niche and population information is proposed to generate new individuals. The proposed NCD-DE is compared with some state-of-the-art and recent well-performing algorithms. The experimental results show that NCD-DE achieves better or competitive performance on both the accuracy and completeness of the solutions than the compared algorithms.

To extend EC algorithms for solving MMOPs, the niching technique is commonly adopted to preserve the diversity of individuals in many studies [19]- [24]. For instance, the niching methods include crowding [22], fitness sharing [25], speciation [26]- [28], clustering [29], [30], hill-valley [31], and recursive middling [32]. The niching technique partitions the entire population into many niches, and each niche deals with one or a few global optima in a small region. A fundamental problem in niching techniques is how to create niches properly. Generally speaking, the difficulties of creating proper niches include how many niches should be partitioned and where the niche centers should be. First, the proper number of niches enables the algorithm to locate all the optima efficiently. Second, the proper location of the niche center can help the individuals belonging to this niche center to speed up converging to the corresponding optimum. However, the number of niches and the locations of the niche centers are usually different in different problems and should be accommodated with the distribution of the population during the different evolutionary states. This difference is attributed to the fact that niches depend on not only the landscape of the problem but also the current population distribution. Therefore, it is a challenge to design a method to distinguish the optimal niche centers (including the number of niches and the location of the niche centers) according to not only the MMOP's landscape but also the distribution of the current population. This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ To address this issue, we propose an interesting and effective idea by transforming the niche center distinguish (NCD) problem into an optimization problem. That is, to distinguish the optimal niche centers to partition the population into a set of promising niches, all the individuals are first regarded as the candidates to be the niche centers. The NCD problem is transformed to a 0/1 binary optimization problem to determine which individuals are selected as the niche centers (i.e., encoded as 1). Since GA can efficiently solve the 0/1 binary optimization problem, an internal GA is designed to optimize the NCD optimization problem. Note that in this transformed optimization problem, both the number of niches (i.e., how many individuals are selected) and the locations of the niches (i.e., which individuals are selected) are not fixed and are determined by the internal GA.
Therefore, this study proposes a ground-breaking and easyextended idea for dealing with MMOPs. Transforming the niche center selecting problem into an optimization problem is a novel method for niching. We do not need to design ad hoc niching approaches for different MMOPs but only design the objective function of the NCD optimization to determine the niche centers, that is, to find out the niche centers with good quality and diversity.
Based on the above idea and method, the NCD problem is solved to provide suitable niche centers and is integrated with the DE algorithm to efficiently solve MMOPs. Therefore, we propose a novel and efficient NCD-based DE (called NCD-DE) algorithm. The characteristics and advantages of NCD-DE are described as follows.
1) During the evolution, the NCD problem is treated as an optimization problem. To efficiently solve the NCD optimization problem, an internal GA is incorporated to select suitable individuals with both good quality and diversity as the niche centers automatically. In the internal GA, a novel fitness-entropy measurement (FEM) that both considers the quality and diversity is proposed to evaluate the selected niche centers. Since the population distribution is different in different problems and different evolutionary stages, dynamically selecting promising individuals as niche centers can satisfy the search requirements of different problems and different evolutionary stages. Note that the chromosomes in internal GA are evaluated via the FEM objective function rather than the objective function of the MMOP, thus the internal GA does not consume the fitness evaluations of the MMOP. 2) A niching and global cooperative mutation (NGCM) strategy that uses both niche-based and population-based information is designed to improve the abilities of exploration and exploitation. On the one hand, a niching mutation strategy with dual-scale local search is used in NGCM to enhance the capability of exploitation to speed up convergence. On the other hand, a global mutation strategy is used in NGCM to enable individuals to jump out from the current stagnate to another potential peak to better explore the search space. 3) During the evolution process, individuals may jump from one niche to another peak, leading to the extinction of some optimal solutions. Therefore, the archive strategy is applied to retain the elite individuals from extinction. The remainder of this article is organized as follows. In Section II, we introduce several relative works on MMOPs and the framework of basic DE are also presented. In Section III, we describe the NCD-DE and analyze the time complexity of NCD-DE. The experimental studies are provided in Section IV. Finally, we conclude and summarize this article in Section V.

A. Multimodal Optimization
To efficiently solve the MMOPs, multiple studies have been devoted in recent years. On the one hand, to locate as many peaks as possible, these proposed algorithms have to maintain good population diversity. On the other hand, since the number of function evaluations (FEs) is limited, convergence speed must also be considered to increase solutions' accuracy. Therefore, the niching method, in which the population is partitioned to form several niches and each niche evolves independently, has been widely incorporated with EC algorithms to search for multiple optima. The famous crowding [22] and speciation [26]- [28] techniques serve as examples. The crowding technique allows the offspring to compete with the most similar individual among the randomly selected C parental individuals (C is the crowd size). If the offspring is fitter, the parent is replaced; otherwise, the offspring is discarded. In the speciation approach, each subpopulation (specie) is formed based on a predefined radius r and each offspring is generated via evolutionary operations within its species.
Although the effectiveness and efficiency of the two niching techniques have been demonstrated in crowding-based DE (CDE) [22] and speciation-based DE (SDE) [26], two limitations remain. First, the results of the two techniques highly depend on the parameter settings. Second, experimental results show that they are not feasible for the problem whose search space is complex [33], [34], such as the highdimensional problem or the problem that contains many local optima. Therefore, some enhanced versions of these two niching techniques have been proposed. For instance, some studies have combined crowding and speciation techniques with the clustering technique to propose the neighborhoodbased CDE/SDE (NCDE and NSDE) [29], as well as the self-adaptive clustering-based CDE/SDE (Self_CCDE and Self_CSDE) [30]. However, in these algorithms, an additional parameter (i.e., clustering size) is also introduced. In addition, the clustering technique is too deterministic to distinguish promising niches in some situations [33]. Even if ignoring the disadvantages of the above algorithms, the effects of the "no free lunch" theorem still encourage researchers to design suitable algorithms for particular problems.
To enhance the performance, some niche parameterinsensitive approaches have drawn considerable attention. For example, an improved information sharing method is incorporated with niching DE (LoINDE) to induce efficient niching behavior in [35]. Yang et al. [33] took advantage of the estimation of distribution algorithm (EDA) and incorporated it with a clustering strategy, forming the LMCEDA and LMSEDA algorithms. Wang et al. [34] proposed an automatic niching DE (ANDE) to form the niches via affinity propagation clustering. Lin et al. [36] developed a DE-based algorithm that adopted nearest-better clustering to partition the population and species balance strategy to fine-tune the species' size. Chen et al. [37] proposed a distributed individuals DE (DIDE) algorithm. In DIDE, each individual is assigned with a virtual subpopulation to search for the optimal solution in a small region. Luo et al. [38] developed a nearest-better-neighbor clustering strategy and incorporated it with the PSO and the evolutionary strategy to better explore and exploit the search space. Huang et al. [39] designed a probabilistic niching framework for MMOPs, in which the historical information is stored in an enhanced binary space partition tree to help enhance the exploration and exploitation abilities.
In addition, some studies have focused on tackling MMOPs via multiobjective techniques. Yao et al. [40] and Deb and Saha [41] proposed to transform the MMOPs into multiobjective problems, in which the objectives are generated according to the fitness and the gradient. To maintain the population's elitism and diversity, Wang et al. [42] proposed a multiobjective optimization for MMOPs (MOMMOPs), in which two conflicting objectives were designed to address MMOPs. Cheng et al. [43] proposed to transform the MMOP into a biobjective optimization problem, whose objective functions are the fitness for MMOP and the grid-based diversity indicator.
Although many MMOP algorithms have been proposed to date, there are still some deficiencies, such as confronting difficulties in tackling high-dimensional problems or manyoptima problems. To overcome these deficiencies, enhancing the exploitation (to reach the elitism of the accuracy of the solutions) ability and improving the exploration (to locate all global optima regions) ability of MMOP algorithms are essential in solving MMOPs.

B. Differential Evolution
The evolution process in DE [12] is executed via the calculations of the difference between individuals. Similar to other EC algorithms, at the initialization stage, individuals are randomly generated in the search space as where x i,j indicates the jth element of the ith individual, and U j and L j are the upper bound and lower bound of the jth dimension, respectively. Here, rand(0, 1) is generated uniformly from interval [0, 1]. After initialization, the evolutionary operations in the main loop of DE, which are mutation, crossover, and selection operators, are executed iteratively before the termination. The mutation operator generates mutation vector v i according to the difference between randomly selected individuals. Some well-known mutation strategies are shown as follows.
In the crossover operator, the trial vector u i is generated via the recombination of the mutant vector v i and the parental individual x i , as follows: where CR indicates the crossover rate, and j rand is a randomly selected integer in [1, D] to avoid the trial vector u i being the same as the parental individual x i . After the crossover, the selection operator is applied to retain the optimal individual from the offspring and the parent. If the fitness of the offspring u i is better than that of the parental individual x i , the individual x i is replaced by u i . When dealing with a maximization optimization problem, the selection operator is shown as follows: where f (.) represents the fitness function.

C. Motivation
How to optimally divide the population into several niches is the key issue to improve the performance of the MMOP algorithms. As the related works reviewed above, some clustering approaches have been proposed to divide the population into several niches. Although the generated niches are suitable for effectively and efficiently solving the MMOP, there are still gaps between the generated niches and the optimal niches for the MMOP.
We note that finding the optimal niches is equivalent to searching for the optimal niches among multiple candidate options. Also, the aiming of solving an optimization problem is searching out an optimal solution among the search space. Considering this, we can build an optimization problem whose optimal solution is corresponding to the optimal niches. However, it is hard to construct an optimization problem that is directly relative to the niche-dividing problem. In MMOP, we can divide the population into niches by first selecting several individuals as niche centers and then assigning the other individuals to the nearest niche centers to generate niches. The problem of selecting optimal niche centers can be intuitively transformed into a binary-encoded optimization problem. Solving the binary-encoded optimization problem can provide a set of promising niche centers, and the niches can be generated according to these niche centers. Based on the above consideration, we can propose a novel and interesting idea that automatically divides the population into several niches via solving the optimization problem. The remained problems only include how to select the suitable optimization algorithm and the objective function to evaluate the niche centers, which are solved in the next section.

III. NCD-DE
This section details the proposed NCD-DE algorithm. In this section, the first part describes the NCD optimization. The NCD optimization problem is solved by a GA, namely, "internal GA," to distinguish the optimal niche centers and partition the population according to the niche centers. In internal GA, the FEM objective function is designed to evaluate the selected niche centers. The NGCM strategy, including the niching mutation strategy with dual-scale local search to improve the exploitation ability and global mutation strategy to enhance the exploration ability, is introduced later. Then, the description of the archive technique is given. Furthermore, the complete NCD-DE algorithm pseudocode and the advantages are provided. Finally, the time complexity of NCD-DE is analyzed.
A. NCD Optimization 1) Binary Encoding of Chromosomes: When using NCD-DE to solve the MMOPs, the NCD optimization is carried out in every generation to determine the optimal niche centers. The NCD problem is transformed into an optimization problem and is solved by an internal GA approach. To avoid confusion, a solution in GA is called a chromosome, whereas a solution in DE is called an individual.
The chromosome of the internal GA is encoded as a 0/1 binary vector with the length as NP, where NP is the number of individuals in DE. Each gene in the chromosome is corresponding to an individual of the DE population. If the value of the gene is 1, then the corresponding individual is distinguished as a niche center; otherwise, if the gene is 0, the corresponding individual is not selected as a niche center. For example, as illustrated in Fig. 1, the DE population with the size of 4 contains individuals A, B, C, and D, and a feasible set of niche centers are individuals B and C, so that the corresponding chromosome is [0, 1, 1, 0]. This way, each chromosome of the internal GA encodes a set of feasible niche centers and the optimal niche centers can be optimized by internal GA. Therefore, the number of niche centers is determined by number of 1s in the chromosome. Note that a niche center does not need to be the best individual in its niche, but the entire set of niche centers needs to make the corresponding chromosome in internal GA have a good quality to partition the population.
2) FEM Objective Function: The quality of a chromosome is determined by the objective function of internal GA. A set of promising niche centers should possess both good fitness and large diversity. First, the fitness of the individuals that are selected as the niche centers should be good enough since a niche center with better fitness means it is more likely to dominate a peak, enhancing the exploitation ability around this niche center. Second, the diversity of all the selected niche centers should be large. That is, because niche centers with more enormous diversity mean that they can cover a broader search space range, which helps discover more global peaks. Therefore, the objective function in internal GA should have two parts, which both focus on the fitness and the diversity of the selected niche centers. The fitness measurement of the selected individuals can be directly reflected by their fitness values on the MMOP, but the diversity measurement should be designed prudently.
Many MMOP algorithms [22], [23] use the Euclidean distance to measure the diversity among individuals. According to the chromosome in GA, the average Euclidean distance among selected niche centers can be calculated as follows: where Ch i,j is the jth gene on the ith chromosome. Since we only calculate the average Euclidean distance among the niche centers, the distance (j) is added when the value of Ch i,j equals 1. NN is the number of niche centers distinguished by this chromosome, which can be calculated as number of 1s in the chromosome. Here, distance(j) is the average distance from the jth individual to other niche centers. The formula of distance(j) is given as follows: where dis(x j , x k ) is the value of the Euclidean distance between individuals x j and x k . Intuitively, suppose the average Euclidean distance of a group of niche centers is relatively larger, niche centers can be distributed in a broad region, where most of the peaks might be covered. Therefore, the average Euclidean distance seems to be a feasible measurement for the niche centers' diversity. However, the average Euclidean distance may not be able to distinguish the most suitable niche centers. That is, because a relatively larger value of the average Euclidean distance only indicates the niche centers are distributed in a broad region, while the niche centers may not be uniformly distributed in  Fig. 2(b), the average distance is 1.307. In such a case, the niche centers in Fig. 2(a) are preferred because their average distance is relatively large. However, as we can see from these two figures and based on common sense, the distribution of niche centers in Fig. 2(b) is preferred because the niche centers are uniformly distributed and dominate a broad region. Therefore, the niche centers in Fig. 2(b) are more dispersed than those in Fig. 2(a), but the Euclidean distance diversity measurement may select niche centers with poor diversity.
To select the most suitable niche centers to maintain diversity, another diversity measurement, the entropy measurement, is employed in our proposed algorithm. The calculation of average entropy among niche centers is shown as where the entropy(j) is the entropy function of the jth individual. The value of entropy of the jth individual can be calculated as where p j,k can be viewed as the similarity between individuals x j and x k . exp(.) represents the exponential function, and ln(.) is the natural logarithm function. Note that the entropy function is calculated only when Ch j,k equals 1. A larger entropy indicates that the niche centers are with higher diversity. When encountered with the case noted in Fig. 2, the entropy measurement of (10) can obtain a better result than the average Euclidean distance measurement. The average entropy of niche centers in Fig. 2(a) is 0.196. In Fig. 2(b), the average entropy is 0.202. Therefore, the selection of niche centers in Fig. 2(b) is crossover rate P c ; mutation rate P m ; iteration epochs T GA ; Output: The best chromosome Begin 1: Randomly initialize and evaluate the population P of size m GA ; 2: For t = 1 to T GA

3:
Initialize two empty sets Q and S;

4:
Evaluate the population P;

5:
For i = 1 to m GA

6:
Select two parents from P by binary tournament selection;

7:
If rand(0, 1) ≤ P c 8: Crossover the parents to generate two offspring, then randomly choose an offspring, mutate it and add it to Q; preferred. This finding indicates that the entropy measurement shows better performance than the Euclidean distance measurement in the above situation. Moreover, to demonstrate the general performance of the entropy measurement, comparison experiments between Euclidean distance and entropy are conducted on commonly used CEC2013 MMOP benchmark functions. The experimental results can be found in Section IV-C. Based on the above case studies, the FEM objective function, including the fitness part and the entropy part for internal GA, is proposed and shown as follows: fitness(j) × entropy(j) (13) where fitness(j) indicates the fitness of individual x j to accomplish the goal of making the fitness of niche centers better, and entropy(j) is the entropy given in (11) between individual x j and other niche centers. The evolution target of each chromosome in the internal GA is to make its FEM value optimal, which indicates that the corresponding niche centers not only obtain better fitness but also yield promising diversity.
3) Framework of Internal GA: To efficiently solve the 0/1 optimization problem above, the internal GA is carried out as a variant of the classical GA [8]- [10]. The internal GA is executed to distinguish a set of niche centers with maximum FEM value at the beginning of each DE generation. Similar to the basic GA, the internal GA maintains a population of chromosomes to search for the optimal niche centers via crossover, mutation, and selection operators.
The process of internal GA is described in Algorithm 1. First, m GA chromosomes that contain information about feasible niche centers are randomly initialized at the beginning, and each gene is a 0/1 binary number (line 1). After the initialization, the main loop (lines 2-18) of the evolution process With the crossover probability P c , the crossover operator is executed on the two parental chromosomes as follows: where Ch p1,j and Ch p2,j are the jth genes of the parental chromosomes, while Ch o1,j and Ch o2,j are the jth genes of the offspring chromosomes. j rand is the index of a randomly selected dimension, indicating the crossover position of the chromosomes. After that, one of the offspring is randomly selected and forwarded to the mutation operator. Otherwise, if the random number rand(0, 1) is larger than P c , the better one of the parental chromosomes is selected and forwarded to the mutation operator without crossover.
In the mutation operator, each gene of the offspring chromosomes can be changed to another value according to the mutation probability P m . Due to the 0/1 binary encoding of the chromosomes, each gene is either 0 or 1. Therefore, if the current gene is 0, its value is changed to 1 after mutation and vice versa, as follows: where P m is the value of the mutation rate. After crossover and mutation, in line 13, the best chromosome of the union set of the population P and the offspring Q is directly added in the next generation population S. Then, other m GA -1 chromosomes of the next generation are chosen via binary tournament selection from the union set of P and Q (lines [14][15][16]. Finally, the internal GA returns the chromosome with the largest FEM value.

4) Population Partition to Form Niches:
Until the completion of the internal GA, the population in DE is partitioned into several niches according to the optimal niche centers distinguished by the internal GA. Algorithm 2 is the pseudocode of the population partition method. First, the niche centers are selected according to the best chromosome Ch Best generated from the internal GA (lines 1-5). Then, each individual that is not a niche center is added to the niche whose niche center is nearest to this individual (lines 6-10).

B. NGCM Strategy
In many MMOP algorithms [23], [34], [36], the nichingbased DE mutation strategy generates a differential vector according to the parents chosen in the same niche. The niching mutation strategy can better exploit the local region to enhance the fitness value of the individuals. If the niching mechanism is not adopted in DE, the global mutation strategy [i.e., DE/rand/1 in (2)] generates differential vectors according to the parents globally selected from all individuals. This way, the search space is enlarged, individuals can migrate from one peak to another remote peak. Therefore, the global mutation strategy in DE can better explore more potential peaks. To combine the advantages of both the niching mutation and the global mutation, the NGCM strategy composed of the niching mutation strategy with dual-scale local search and the global mutation strategy is designed.
First, the niching mutation strategy with dual-scale local search is described, which corresponds to lines 7-22 in Algorithm 3. If the niche size is not smaller than three, the niching mutation strategy is adopted in this niche to generate offspring as follows: where nr1, nr2, and nr3 are three different indexes that are randomly selected in the niche of the individual x i . After niching mutation, the crossover operator in (6) is executed to recombine a new offspring. Once the offspring is created, it is forwarded into the crowding-based selection operator. The crowding-based selection strategy is suitable here for the following reasons. As the generated offspring u i may be far away from its parental individual x i , directly applying the selection operator in (7) may not be reasonable because u i and x i are more likely in the different niches. In such a situation, the crowding-based selection strategy that selects a better one within the offspring u i and its nearest individual x nn may be more suitable. When dealing with a maximization optimization problem, the crowding-based selection operator is shown as where x nn is the most similar parental individual of u i . Herein, the most similar refers to the nearest individual measured by the Euclidean distance. However, the niching mutation strategy cannot be carried out when the niche size is smaller than three. In such a situation, the dual-scale local search is performed to generate offspring. The niching mutation strategy is carried out by following (19) and (20), meaning the wide-scale and narrowscale, respectively. Specifically, with a half probability, the offspring is generated by a wide scale as follows: Copy the best individual in S i into the archive; Niching mutation strategy via (17);

16:
Create offspring by wide-scale local search in (19) where x nn is the most similar individual of x i . The normal (0, 1) is a random number sampled from the standard normal distribution, and dis(x i , x nn ) is the Euclidean distance between individuals x i and x nn .
With another half probability, the new offspring is generated by a narrow scale as follows: where nn is the index of x i 's nearest neighbor.
In the dual-scale local search strategy, as the offspring u i is sampled from a normal distribution centralized at its parental individual x i , the offspring u i is very similar to its parental individual x i . Therefore, the crossover operator between the new offspring u i and its parent x i is not necessary.
Similarly, as the new offspring u i is very similar to the parental individual x i , the crowding-based selection is also not necessary in this case. Therefore, the basic selection operator in (7) is carried out, where the comparison is performed between the offspring u i and its parental individual x i .
Second, the description of the global mutation strategy is given here, which corresponds to lines 25-29 in Algorithm 3. The global mutation strategy is executed periodically (i.e., by using the DE/rand/1 mutation strategy every λ generations) in the NGCM to help better discover more potential peaks. Specifically, after the interval of λ generations, each individual is updated via the global DE/rand/1 mutation strategy [i.e., (2)], and then the crossover operator in (6) and the crowding-based selection in (18) are carried out to generate offspring. Since each individual can mutate with any other individuals among the entire population, the individuals can vary in a broader region among the search space. Therefore, the global mutation strategy enables individuals to discover more potential peaks.

C. Archive Technique
Multiple algorithms designed for multiobjective optimization problems adopt an archive to preserve nondominated solutions [44]- [47]. Recently, this technique has also been used in a few MMOP algorithms to save historical information [48] or store stagnated individuals and reinitialize it to enhance diversity [49], [50].
In the proposed NCD-DE, it is observed that some niches may be extinct after the global mutation strategy because the individuals may jump to other niches or potential peaks. Even in the niching mutation strategy with dual-scale local search, the individuals may also jump to other nearby niches. Therefore, the archive technique is adopted in our algorithm to store optimal solutions. The idea of the archive technique is, after the population partition and before the NGCM strategy, the best individual in each niche is copied and stored into the archive to protect the found optimum. The effect of the archive technique is investigated in Section IV-C.

D. Complete NCD-DE
The complete NCD-DE algorithm is given in Algorithm 3. For a better understanding of the proposed algorithm, we illustrate the relationship among the NCD optimization, internal GA, and DE with NGCM strategy in Fig. 3. The advantages of the NCD-DE are summarized as follows.
1) Transforming the NCD problem into an optimization problem enables the NCD-DE to dynamically determine the number and locations of the niche centers depending on not only the optimization problem but also the distribution of the population, enhancing the performance on different problems and during different stages of the evolution. 2) Due to the FEM objective function in internal GA, the NCD optimization can select niche centers with both better quality and diversity. Compared with the average Euclidean measurement, the entropy measurement in the FEM objective function can better measure the diversity of the niche centers, which helps to distinguish the promising niche centers.
3) The NGCM strategy can enhance exploration and exploitation. The niche mutation strategy helps to exploit the nearby region, and the convergence speed is accelerated with the dual-scale local search strategy. The global mutation strategy can discover more potential optima to enhance exploration. 4) In the evolutionary process of NCD-DE, some individuals can jump into another peak, which may lead to the extinction of some niches. The archive technique can avoid the extinction of some niches by copying and storing the best individual of each niche in the archive.

E. Complexity Analysis
To For clarity, the time complexities of the recent wellperforming algorithms are listed in Table S.I in the supplementary material. Due to the space limitation, Table S.I is shown in the supplementary material. In Table S.I in the supplementary material, the time complexity of most algorithms is O(NP 2 ×D), the time complexity of MOMMOP algorithm is O(2NP 2 ×D 2 +4NP 2 ×D), and that of ANDE is O(NP 2 ×D+NP 3 ×NT), where NT is the iteration times of affinity propagation clustering in ANDE.

A. Experiment Settings
The experiments are conducted in the well-known CEC2013 MMOP benchmark functions to validate the effectiveness and efficiency of the NCD-DE. CEC2013 MMOP benchmark contains 20 functions with multiple global optima [51]. The basic properties of the 20 CEC2013 benchmark functions are provided in Table S.II in the supplementary  material. Two commonly adopted performance metrics: 1) peak ratio (PR) and 2) success rate (SR), are adopted to evaluate the performance of NCD-DE and the compared algorithms. Under the settings of the accuracy level ε and the maximum FEs (MaxFEs), PR is calculated as the mean proportion of the found global optima over all the global optima in several runs of the algorithm, and SR is calculated as the percentage of the successful runs where all the global optima are found over multiple runs. The formulas used to compute PR and SR can be found in [51].
For the settings of accuracy level, three accuracy levels ε = 1.0E -03, ε = 1.0E -04, and ε = 1.0E -05 are adopted in the experiments to evaluate the performance of NCD-DE in different accuracy levels. The detailed settings of maximum FEs MaxFEs and population size NP on different problems are shown in Table S.III in the supplementary material [33], [34]. For the parameter settings in the NCD-DE, the population size m GA of the internal GA is 30, and the crossover rate P c and mutation rate P m of the internal GA are, respectively, specified as 0.9 and 0.1, which are two common parameter settings for GA. The iteration number T GA of the internal GA is set as 5. In the mutation and crossover operators of the NGCM strategy, the amplification factor F is 0.5 and the crossover rate CR is 0.9. The period λ of the global mutation strategy is set as 5. For the experimental results, each algorithm is independently executed 51 times to obtain the mean results.

B. Comparison With State-of-the-Art MMOP Algorithms
The experimental results for PR and SR at three different accuracy levels (i.e., ε = 1.0E -03, ε = 1.0E -04, and ε = 1.0E -05) obtained by the NCD-DE on F 1 -F 20 are provided in Table I. For clarity, the results where PR equals 1.0 are emphasized by boldface.
To further demonstrate the effectiveness and efficiency of NCD-DE, we compare the proposed NCD-DE with a group of state-of-the-art MMOP algorithms. These algorithms include nine DE-based algorithms, which are CDE [22], SDE [26], NCDE, NSDE [29], Self-CCDE, Self-CSDE [30], LoICDE, LoISDE [35], and ANDE [34]. Two EDA-based MMOP algorithms, including LMSEDA and LMCEDA [33], are also compared. Furthermore, a multiobjective-based MMOP algorithm MOMMOP [42] is included for comparison. Besides, NMMSO [52] that wins first place in the CEC2015 MMOP competition is also compared with NCD-DE. Note that the benchmark functions of the CEC2015 MMOP competition are the same as the CEC2013 MMOP benchmark functions [34]. The parameter settings for all the comparison algorithms are adopted from their corresponding papers. The NP and MaxFEs of all these algorithms are provided in Table S.III in the supplementary material. The results of the compared algorithms are derived from [34] and [50].
The experimental results of the NCD-DE and other algorithms at ε = 1.0E -04 are provided in Table II, and the rest  TABLE I EXPERIMENTAL RESULTS OF NCD-DE ON CEC2013 BENCHMARK FUNCTIONS F 1 − F 20 AT ε = 1.0E -03, ε = 1.0E -04, AND ε = 1.0E -05 of the results at ε = 1.0E -03 and ε = 1.0E -05 are provided in Tables S.IV and S.V in the supplementary material, respectively. To evaluate the performance of each algorithm in the statistic view, Wilcoxon's rank-sum test [53] at α = 0.05 is adopted on the values of PR between NCD-DE and other algorithms. The symbols "+ / ≈ / −" indicate the PR obtained by NCD-DE is "significantly superior/ equal/significantly inferior" to that obtained by the other algorithm, respectively. The analysis of the comparison results is given as follows.
1) For the first ten low-dimensional benchmark functions F 1 -F 10 , the experimental result of NCD-DE is significantly better than those obtained by most algorithms, while NMMSO dominates NCD-DE on F 7 and MOMMOP dominates NCD-DE on F 7 , F 8 , and F 9 . However, as the accuracy level increases, the results of MOMMOP rapidly decrease, while the results obtained by NCD-DE are more stable. 2) For the next five functions F 11 -F 15 , which are lowdimensional functions with many optima, NCD-DE also achieves promising performance. Although NCD-DE does not achieve the best PR on every function, the experimental result of NCD-DE is significantly better than those of most algorithms. Compared with the results on different accuracy levels, many algorithms' performance deteriorates as the accuracy level increases, such as CDE, SDE, NCDE, Self-CCDE, Self-CSDE, LoICDE, and MOMMOP. However, the performance of NCD-DE only slightly decreases with increasing accuracy levels. 3) On the last five high-dimensional functions F 16 -F 20 , NCD-DE also outperforms most of the compared algorithms. Although the experimental result of NCD-DE is dominated by that of LMSEDA on F 17 , NCD-DE achieves significantly superior performance than most algorithms. On these five high-dimensional functions, the results of most algorithms rapidly decrease as the accuracy level increases. For example, compared with the results on accuracy level ε = 1.0E -03 and ε = 1.0E -05, it is evident that the performance of many algorithms (such as Self-CCDE, Self-CSDE, and LMCEDA) deteriorates rapidly when the accuracy level increases. In contrast, NCD-DE's performance is stable. According to the results and analyses above, the effectiveness and efficiency of NCD-DE are demonstrated. NCD-DE outperforms most of the compared algorithms not only on the low-dimensional functions but also on the high-dimensional functions. Furthermore, the stable results indicate that the results obtained by NCD-DE do not deteriorate rapidly as the increasing of accuracy level ε, whereas the experimental results of other algorithms decrease quickly as ε increases. Therefore, it is concluded that the performance of NCD-DE is stable. Besides, the NCD-DE algorithm can find almost global optima and ensure the accuracy of solutions, which indicates that NCD-DE can achieve both great exploration and great exploitation abilities.

C. Effects of Components in NCD-DE
There are three independent components in NCD-DE: 1) the NCD optimization; 2) the NGCM strategy; and 3) the archive technique. Their properties and effects are analyzed in this section.

1) Effect of NCD Optimization:
The effect of NCD optimization is mainly influenced by the FEM objective function. Although the motivation of the proposed FEM objective function in maintaining the diversity of niche centers has been discussed in Fig. 2 and its effectiveness in helping NCD-DE solve MMOPs has been demonstrated in Section IV-B, the effect of the FEM in the NCD optimization is further studied here. To demonstrate the FEM is superior to the fitness-distance measurement (FDM), a comparison is conducted between the NCD-DE and a variant of NCD-DE with FEM objective function. The variant of NCD-DE is called NCD-DE-FDM (NCD-DE with FDM objective function). In NCD-DE-FDM, the average Euclidean distance between the niche centers is applied to replace the entropy in the FEM. The FDMs objective function FDM(.) can be calculated as follows: fitness(j) × distance(j) (21) where the formula of distance(j) is shown in (9). Another NCD-DE variant called NCD-DE-bq (NCD-DE with the best quarter of the population as niche centers) is designed to demonstrate the importance of the diversity of the selected niche centers. In NCD-DE-bq, the diversity of the niche centers is ignored. Instead, a quarter (i.e., 25%) of the total individuals with optimal fitness value are selected as niche centers. The comparison between NCD-DE and NCD-DE-bq reveals the necessity of the diversity of niche centers.  Table III Table  S.VII in the supplementary material. It is easy to be observed from the experimental results that the NCD-DE achieves more promising performance than NCD-DE-w/o-a. Specifically, the experimental results of NCD-DE are significantly better than those of NCD-DE-w/o-a in functions F 7 and F 9 , which contain many global optima (F 7 contains 36 global optima, and F 9 contains 216 global optima). In these functions, it is easy for individuals to jump from one peak to another, which leads to the extinction of several individuals. Therefore, the comparison result indicates that the archive technique helps to preserve individuals from extinction.

1) Sensitivity Analysis of T GA :
The control parameter T GA is an integer that determines the internal GA's iteration times in the NCD optimization. The setting of T GA should be carefully considered. On the one hand, if the value of T GA is relatively small, the internal GA cannot obtain promising niche centers with optimal FEM values in only a few iterations. On the other hand, if T GA is set as a relatively large value, the optimization for optimal niche centers can cost a relatively large computation complexity.
To investigate the influence of the parameter T GA , four NCD-DEs with different values of T GA are designed, which are NCD-DEs with T GA = 5 (original), 10, 15, and 20. The experimental results at ε = 1.0E -04 of all the NCD-DE variants with different values of T GA are provided in Table S.VIII in the supplementary material.
From Table S.VIII in the supplementary material, the experimental results of the NCD-DE variants are slightly improved as T GA increases in F 7 , F 9 , F 12 , F 15 , and F 20 . However, in F 8 , F 13 , and F 17 , the PRs perform undulate results as T GA increases. The performance of NCD-DE-T20 (i.e., NCD-DE with T GA = 20) for PR is generally better than that of NCD-DE, which can be concluded from the result of Wilcoxon's rank-sum test. However, NCD-DE outperforms NCD-DE-T10 and NCD-DE-T15.
In summary, the results of NCD-DE slightly increase as the T GA parameter increases. This finding is attributed to the fact that the internal GA chromosomes are refined during the evolutionary process iterations. If the iteration time increases, the final generation of chromosomes may obtain better FEM values. However, the performance increase is slight. Considering the balance between the time complexity and the algorithm's performance, the setting of the parameter T GA is recommended as 5.
2) Sensitivity Analysis of λ: In the NGCM strategy, the control parameter λ is an integer that determines the gap for the execution of the global mutation strategy. The value of parameter λ should be carefully selected. On the one hand, the global mutation strategy is executed frequently if λ is too small, so as to cost a large amount of FEs. On the other hand, if λ is relatively large, the global mutation strategy is seldomly executed, which can deteriorate the exploration ability of NCD-DE.
Therefore, the influence of the parameter λ should be investigated. Four NCD-DEs with different λ are designed, which are NCD-DEs with λ = 1, 5 (original), 10, and 20. These NCD-DEs are called NCD-DE (where λ = 5), NCD-DE-λ1, NCD-DE-λ10, and NCD-DE-λ20, respectively. The experiment results at ε = 1.0E -04 of all the NCD-DEs with different values of λ are given in Table S.IX in the supplementary material.
On functions F 1 -F 6 , all these NCD-DEs with different values of λ can find all the global optima. However, on F 7 -F 9 , the performance of NCD-DE-λ1 dominates that of other NCD-DEs. On F 12 -F 19 , the performance of NCD-DE dominates that of other NCD-DE variants. On F 20 , the best result is obtained by NCD-DE-λ20. From the experimental results, generally, the NCD-DE achieves the best performance. For example, the experimental results of NCD-DE with λ = 5 are significantly better than those of NCD-DE-λ1, NCD-DE-λ10, and NCD-DE-λ20 on seven functions, seven functions, and eight functions, respectively. We can find that the NCD-DE achieves the best performance when λ = 5. Therefore, the value of λ is recommended to be set as 5.

E. Real-World Application
To verify the effectiveness and efficiency of NCD-DE, the experiments on a real-world MMOP, that is, multiple competitive facilities location and design (MCFLD) problem [54], are conducted in this section. In the MCFLD problem, a chain wants to construct several new facilities in a given area. However, with the condition that many facilities offering the same service already exist or might be soon constructed in the given area, the new facilities of the chain have to compete with the existing facilities in the market. The goal of the MCFLD is to determine the optimal location and quality of the new facilities to maximize the total profit. Since the MCFLD problem is a nonlinear problem with multiple optima, the chain may want to find out several available solutions in the real-world application.
In specific, according to the model of MCFLD problem given in [54], the total profit that constructs p new facilities is formulated as M nf 1 , nf 2 , . . . , where n, m, and k are the number of demand points, the number of existing facilities, and the number of facilities that belong to the chain, respectively. nf j stands for the properties of the jth new facility, which includes three variables: 1) quality α j ; 2) x-coordinate x j ; and 3) y-coordinate y j . Thus, the MCFLD problem includes 3p decision variables. The detailed description and settings of the parameters can be found in [54].
Since the location and the fitness of the real optima of the MCFLD problem are unknown, we first estimate the optimal fitness and then determine which solution is located at the peak. To estimate the real optimal fitness, the entire search space is divided into 10 6 grids, and a position is sampled in each grid. The fitness value of the best position in the 10 6 sampled positions is the estimated optimal fitness. A solution is determined to be peak if its fitness is larger than the estimated optimal fitness and there is no solution determined as a peak within 0.1 distance.
The results of mean peak number (PN) and the corresponding rank obtained by NCD-DE are compared to those obtained by CDE, NCDE, and ANDE. Each algorithm is independently executed 20 times to obtain the mean results. The experimental results on six MCFLD problems with settings p = 2 and (n, m, k) = (21, 5, 2), (50, 2, 1), (50,10,4), (71, 5, 2), (100, 2, 1), and (100, 10, 4) are given in Table IV. From the results, it is obvious that the NCD-DE can find more global optima than the compared algorithms. In each problem, the NCD-DE can find more than 60 optima, while the number of optima found by CDE and NCDE is less than 5. Although the ANDE can find out nearly 40 optima when (n, m, k) = (21, 5, 2), its performance in the other cases is not promising enough.

V. CONCLUSION
In this article, we proposed an inspiring and practical mechanism that transforms the NCD problem into an optimization problem, and applies an extra EC algorithm to solve it. Through the transformation, we do not need to design ad hoc niching approaches for different MMOPs but only design the objective function of the optimization problem to fulfill our demands of the niche centers.
Based on the above mechanism, a novel NCD-DE algorithm composed of the NCD optimization, NGCM strategy, and archive technique has been proposed for solving MMOPs. In the NCD optimization, the NCD problem is transformed into an optimization problem and solved by the internal GA. The NCD optimization calculates the FEM of the selected niche centers to ensure that the niche centers with the optimal fitness value and great diversity are selected. The niche centers optimally selected by the NCD optimization are distributed widely and are with good fitness values on the MMOPs, which helps to locate more peaks and to increase the solution accuracy. To further increase the capability of exploration and exploitation, the NGCM strategy was proposed. On the one hand, the niching mutation strategy with dual-scale local search helps to search the circumambient region of each niche. On the other hand, the global mutation strategy enables each individual to jump from one niche to another niche or some undiscovered region. In addition, the archive technique was proposed to copy each niche's best individual and store it into an archive, which preserves the optimal individuals from extinction.
Based on the three components described, the proposed NCD-DE achieved better performance than multiple state-ofthe-art and recent well-performing algorithms in dealing with MMOPs, as verified on the CEC2013 benchmark functions.
In future work, we will comprise the NCD optimization with other EC algorithms, such as PSO and EDA, to propose more alternative algorithms to deal with MMOPs. Moreover, we will make efforts to extend the NCD-DE to solve realworld MMOP applications. Besides, extensive research into updating the proposed NCD-DE as a distributed [13], [55] or parallel [56] algorithm, or using matrix-based framework [57], scale-adaptive fitness evaluation [58], or region-based encoding scheme [59] to reduce the running time are also worthy studied.