Region-based Memetic Algorithm with Archive for multimodal optimisation

Accepted Manuscript Region-based Memetic Algorithm with Archive for multimodal optimisation Benjamin Lacroix, Daniel Molina, Francisco Herrera PII: S0020-0255(16)30388-7 DOI: 10.1016/j.ins.2016.05.049 Reference: INS 12270 To appear in: Information Sciences Received date: 5 March 2015 Revised date: 13 April 2016 Accepted date: 29 May 2016 Please cite this article as: Benjamin Lacroix, Daniel Molina, Francisco Herrera, Region-based Memetic Algorithm with Archive for multimodal optimisation, Information Sciences (2016), doi: 10.1016/j.ins.2016.05.049 This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain. ACCEPTED MANUSCRIPT ACCEPTED MANUSCRIPT Region-based Memetic Algorithm with Archive for multimodal optimisation Benjamin Lacroixa, Daniel Molinab, Francisco Herrerac,d aIDEAS Research Institute, Robert Gordon University, United Kingdom bDepartment of Computer Science, University of C´adiz, 11003, C´adiz, Spain cDepartment of Computer Science and Artificial Intelligence, University of Granada, 18071, Granada, Spain dFaculty of Computing and Information Technology, King Abdulaziz University, 21589, Jeddah, Saudi Arabia Abstract In this paper we propose a specially designed memetic algorithm for multimodal optimisation problems. The proposal uses a niching strategy, called region-based niching strategy, that divides the search space in predefined and indexable hypercubes with decreasing size, called regions. This niching technique allows our proposal to keep high diversity in the population, and to keep the most promising regions in an external archive. The most promising solutions are improved with a local search method and also stored in the archive. The archive is used as an index to effiently prevent further exploration of these areas with the evolutionary algorithm. The resulting algorithm, called Region-based Memetic Algorithm with Archive, is tested on the benchmark proposed in the special session and competition on niching methods for multimodal function optimisation of the Congress on Evolutionary Computation in 2013. The results obtained show that the region-based niching strategy is more efficient than the classical niching strategy called clearing and that the use of the archive as restrictive index significantly improves the exploration efficiency of the algorithm. The proposal achieves better exploration and accuracy than other existing techniques. Keywords: Multimodal optimisation, memetic algorithm, niching strategy 1. Introduction1 Many real world problems offer various solutions considered as global op-2 tima. The identification of multiple solution has thus gained popularity in the3 research community. It is referred to as multimodal optimisation as the objec-4 tive is to retrieve more than one optima. While classical evolutionary algorithms5 Email addresses: [email protected] (Benjamin Lacroix), [email protected] (Daniel Molina), [email protected] (Francisco Herrera) Preprint submitted to Information Sciences June 1, 2016 ACCEPTED MANUSCRIPT ACCEPTED MANUSCRIPT (EA) were designed to identify a single optimum, some modifications have to6 be applied to identify multiple optima, preventing their premature convergence7 and maintaining the diversity in their population to ensure the exploration of8 distinct areas of the fitness landscape. Such techniques, known as niching strate-9 gies [6], are meant to stay in the population subgroups of individuals, or niches,10 in different parts of the search domain.11 Most existing techniques’ efficiency relies on two problem dependent param-12 eters, the niche radius and the population size [7, 16, 42]. The first one should13 be defined according to the distance between optima in the fitness landscape14 and the second one according to the number of optima to locate. Both data are15 however usually unknown in real world problems. Nowadays, research interest16 focuses on designing EA which are less dependent on those parameters.17 The main challenge when designing an EA for multimodal optimisation is to18 create an algorithm capable of approximating with the highest level of accuracy19 the different global optima.20 Memetic algorithms (MA) [35] are the hybridisation between EA and local21 search methods (LS) combining in one model the exploration power of the for-22 mer and the exploitation capacity of the latter. This hybridisation can achieve23 a good trade-off between the exploration of the domain search and the exploita-24 tion of found solutions, so it is important to obtain good results in EAs [59],25 and it also offers interesting properties when applying them to multimodal opti-26 misation problems from the multimodal optimisation point of view. Indeed, as27 we said before, niching techniques used with classical EA forms sub-populations28 destined to explore and optimise different areas of the search space with the29 same mechanism. MA separate these efforts, leaving the exploration task to the30 EA and the refinement of the most promising regions identified by the EA to31 the LS method.32 In a previous work [21], we designed a MA for global continuous optimisa-33 tion problems called region based memetic algorithm with local search chaining34 (RMA-LSCh). It proposed a novel niching strategy, the originality of which35 lies in the definition of a niche. While traditionally the niche surrounding a36 solution is defined by the radius around it, the proposed niching technique par-37 titions the search into equal hypercubes called regions. The dependency to the38 niche size (defined by the number of divisions of the search space) is reduced39 by increasing the number of divisions during the search. In this work we pro-40 pose a new algorithm specially designed for continuous multimodal optimisation,41 Region-based Memetic Algorithm with Archive (RMAwA). Although RMAwA42 maintains the same definition of a niche and alternatively applies the EA and43 the LS, the memetic scheme is modified and a novel archive is implemented to44 match the requirements of multimodal optimisation. First, while RMA-LSCh45 uses LS Chaining [32, 33] and thus limites the number of fitness evaluation per46 LS application, RMAwA applies the LS until it has reached a local or global47 optimum. Most importantly, regions intensively explored by LS are discarded48 by the proposal from further exploration. RMAwA contains an indexed archive49 with these regions to reduce the search domain in a very efficient way. Also,50 because the identified optima are stored into the archive and not into the pop-51 2 ACCEPTED MANUSCRIPT ACCEPTED MANUSCRIPT ulation, the number of optima that RMAwA can identify is not limited by the52 population size [12, 63, 64].53 RMAwA is tested using a specific benchmark for multimodal optimisation.54 The experiments carried out show that the use of the region based niching strat-55 egy coupled with an archive provides interesting improvements to the memetic56 framework, and that the RMAwA is a very competitive algorithm against ex-57 isting ones.58 This paper is organised as follows. In Section 2, we present a quick intro-59 duction on methods previously proposed to tackle multimodal problem optimi-60 sation. In Section 3, we present the RMAwA and detail each component. In61 Section 4, we explain the experimental framework used and the parameter set-62 ting of the algorithm. In Section 5, several comparisons are carried out to study63 the influence of the different components of the algorithm and our proposal is64 compared with other algorithms in the literature. Finally, In Section 6 some65 concluding remarks are pointed out.66 2. Background67 In order to identify multiple optima of a fitness landscape several techniques68 have been proposed. In this section, we give a brief overview of techniques that69 have been proposed to maintain the diversity in the population in order to pre-70 vent its convergence towards a single optimum. Such techniques are commonly71 called niching strategies and refer to the technique used for the discovery and72 preservation of distinct niches. This term is a reference to the ecological concept73 of niches referring to the formation of distinct species exploiting different niches74 (resources) in an ecosystem.75 The main challenge in multimodal optimisation is the unknown nature and76 characteristics of the objective function, specifically the number of global optima77 and their repartition on the search domain. The main goal of the proposals78 presented in this section is to tackle these issues. Alternatively, [55] proposes a79 preprocessing tool to estimate the number of basins of attraction in the fitness80 landscape.81 We have classified the methods proposed to tackle multimodal optimisation82 into two categories. The first one lists the classical niching strategies which83 mainly affect the replacement criterion of the EA they are applied to. The84 second one works with the idea of creating subgroups of solutions in different85 area of the search space by limiting the cooperation of each individual to its86 nearest neighbours. We refer to them as neighbourhood based techniques.87 In this section, we first describe the different elements composing those two88 categories by giving a general overview of the proposal making use of such89 techniques. In a third section, we briefly introduce proposals combining those90 techniques with MA which demonstrate that the use of a refinement method91 improves the performance of EAs for multimodal optimisation.92 3 ACCEPTED MANUSCRIPT ACCEPTED MANUSCRIPT 2.1. Classical niching techniques93 The first niching techniques consist in limiting the presence of multiple so-94 lutions within the same niche in order to keep the population highly diverse.95 When included in a classical EA, those mechanisms are mainly replacement96 strategies designed to remove solutions present in the same vicinity. We de-97 scribe here the four main methods to achieve this objective: crowding, clearing,98 fitness sharing, and speciation.99 2.1.1. Crowding100 Crowding is one of the first techniques proposed to tackle multimodal opti-101 misation problems [7]. After the generation of a new solution, a random sample102 of CF solutions is selected in the population. Each new solution competes with103 the closest solution of the sample to stay in the population. This technique’s104 main drawback is the definition of the crowding factor parameter (CF). A small105 value can lead to the replacement of a distant solution to the offspring and thus106 a loss of information, and a very large value has a high computational cost. The107 efficiency of this technique has proven to be limited [30] and advanced versions108 have been proposed:109 •Deterministic crowding proposed by [30] tries to limit the problem of110 replacement errors induced by the crowding technique by eliminating the111 need of defining the CF parameter. To do so, an offspring competes with112 its own parents to stay in the population.113 •Probabilistic crowding [31] on the other hand modifies the replacement114 strategy of the original technique. In this scheme, the offspring and its115 most similar individual in the crowding sample compete in a probabilis-116 tic tournament where the probabilities of winning for each individual X,117 p(X), is calculated according to their fitness:118 p(X) = f(X) f(X) + f(Y)(1) where f(X) is the fitness of the same solution Xand f(Y) is the fitness of119 the other solution. The idea is not to always show preference to solutions120 with higher fitnesses which may lead to the loss of niches.121 In [57], Thomsen proposed the popular crowding differential evolution (CDE)122 applying a classical crowding strategy on a differential evolution (DE) where a123 new solution is created by means of classical DE mutation and crossover scheme124 comparing with its closest solution in the whole population for replacement.125 CDE was then extended to multi-population crowding DE (MCDE) in [63]126 where multiple sub-population evolve in parallel using CDE. When all the sub-127 populations have converged, the optima identified by each of them are stored in128 an archive and the sub-populations are reinitialised.129 4 ACCEPTED MANUSCRIPT ACCEPTED MANUSCRIPT More recently, Qu et al. proposed the dynamic grouping of CDE (DGCDE)130 [45] with ensemble of parameters. The population is divided into three sub-131 population to which a set of control parameters is assigned.132 In [44], Qing et al. proposed a Crowding Clustering Genetic Algorithm133 (CCGA) using a clustering technique to eliminate the genetic drift introduced134 by the crowding strategy.135 2.1.2. Clearing136 Clearing techniques [42] lie in the principle of dedicating the limited re-137 sources of a niche to its best individuals. The population is sorted according138 to the individual fitness values. The solutions are then selected one after the139 other and the solutions with worse fitness falling within their niche radius σclear 140 are removed. Clearing has a low complexity and shows the best performances141 amongst the classical techniques but is highly sensitive to the niche radius [51].142 Variations have then been proposed to limit influence of the σclear parameter.143 For instance, in [47], similarly to the previously cited DGCDE, the authors144 propose an ensemble of clearing DE (ECLDE) in which the population was145 equally divided into 3 sub-populations each evolving in parallel using a clearing146 DE with different values of σclear.147 Some techniques use a redefinition of the niche in order to remove the use148 of the parameter σclear. In [11], the niches are defined through a hill-valley149 detection mechanism instead of using a niche radius. In [50], the niches are150 defined by fuzzy clustering of the solutions of the populations.151 2.1.3. Fitness sharing152 Contrarily to the clearing technique which consist in dedicating niche resources to a single solution, fitness sharing [16] consists in reducing the fitness of individuals present in densely populated regions. The fitness used of the ith individual, fshared(i), is calculated by: fshared(i) = foriginal(i) PNP j=1 sh(dij)(2) where foriginal is the original fitness function, NP is the population size, and sh function is calculated by: sh(dij) = (1−dij σshare α ,if dij < σshare 0,otherwise (3) where dij is the distance between individual iand j,σshare is the sharing153 radius and αis a constant called sharing level.154 In [57], Thomasen also proposed an DE using sharing where, after each gen-155 eration, the new shared fitnesses are calculated over the population individuals156 and the trial vectors, the best half being kept in the population.157 5 ACCEPTED MANUSCRIPT ACCEPTED MANUSCRIPT 2.1.4. Speciation158 Proposed in [22], speciation or species conservation introduces the notion of159 species by separating the population into several groups (species) according to160 their similarity. Those species are identified by a dominating individual called161 the species seed and a species distance σspecies defining the maximum distance162 between two individual of the same species. The set of species seed is build at163 each generation by iteratively adding individuals from the population that are164 further from any species seed than σspecies/2. The individuals are kept from one165 generation to another until a better solution is identified within their species166 while the classical recombination operators are applied.167 In [23], this concept is applied to a speciation-based PSO (SPSO). In SPSO,168 the particles are gathered into species to form sub-populations. This proposal169 was later extended to reduce its dependency to the species distance parameter170 by using population statistics [3] and a time-based convergence measure [49].171 2.2. Neighborhood based technique172 Another class of niching strategies can be referred to as neighbourhood-173 based. Contrarily to the previous section where the niching strategy could be174 seen as replacement strategy, these methods use the geographical information175 of the solutions in a population to modify the recombination scheme of a given176 EA. The main idea is to make solutions by only considering their neighbours in177 order to emphasize the speciation.178 Originally named spatially-structured EAs (SSEA) [58], these algorithms179 form sub-populations of individuals (called deme) based on their similarity and180 perform genetic operations within each deme.181 This idea has then been extended and two kinds of neighbourhoods can be182 identified in the literature:183 •Index-based neighbourhood [24] uses the indices in the population of a PSO184 to identify the neighbourhood of a solution. The velocity of a particle is185 thus influenced by the local best solution instead of the global best.186 •Distance-based neighbourhood uses the euclidean distance between individ-187 uals. In [26], the author proposed the FER-PSO algorithm where parti-188 cles are attracted towards the ”fittest-and-closest” neighbours. Similarly,189 the notion of neighbourhood is applied for DE in [13]. A new mutation190 strategy, DE/nrand/x is proposed. It uses as a base vector the nearest191 neighbour of each individual. This mutation strategy has then been used192 for more advanced models like in [12]. In [4] a neighbourhood mutation193 is proposed that considers normalized distance. Another option is to use194 the distance to create a clustering partition of the population to maintain195 diversity [15].196 Neighbourhood-based strategies have often been coupled with classical nich-197 ing strategies. For instance in [10], the authors propose including in a SSEA a198 fitness sharing and a clearing strategy.199 6 ACCEPTED MANUSCRIPT ACCEPTED MANUSCRIPT In [48], the authors use the DE/nrand/x operator with crowding, sharing200 and species-based niching strategies and obtain better results than the original201 algorithms.202 2.3. Memetic algorithms for multimodal optimisation203 As stated in the introduction, MA are the hybridisation of an EA and a204 LS method. This model is part of the more general Memetic computing (MC)205 family of algorithms which combine various optimisers (memes). The efficiency206 of these models have helped them gain popularity over the past decade [5, 37].207 The coordination of the memes is the main research topic in MC. Ong et al.208 [40] proposed a classification which was later updated by Neri et al. [37]:209 •Adaptive Hyper-heuristic [19]: the memes are coordinated by means of210 heuristic rules.211 •Meta-Lamarckian learning [39]: the probabilities of using the memes are212 based on their success, providing an online adaptability.213 •Self-Adaptive and Co-Evolutionary [20, 54]: the memes are encoded with214 the candidate solutions and evolve in parallel so the most appropriate can215 be selected.216 •Fitness Diversity-Adaptive [38]: the selection of the memes to be operated217 is based on the diversity measure of the population.218 MA are particularly adapted to multimodal optimisation problems as, when219 applied to different solutions, an LS method can offer a strong refinement of220 the promising solutions discovered by the EA, providing great accuracy for the221 identification of multiple optima. The use of such model has raised interest in222 the research community.223 For instance, the Sequential Niching Memetic Algorithm (SNMA) proposed224 by Vitela et al. in [60] and then extended in [61] is an MA which combines a ge-225 netic algorithm (GA) with a gradient-based LS method. Before each generation,226 the LS is applied to each solution of the population.227 In [46], Qu et al. included an LS method to various previously cited PSO for228 multimodal optimisation (FER-PSO, SPSO, rPSO). The LS method used con-229 sisted in generating at each iteration new solutions in the neighbourhood of the230 personal best of each particle to explore its surrounding. They demonstrated231 that the resulting memetic PSO obtained better results than the original algo-232 rithms. Similarly, Wang et al. proposed a memetic SPSO [62] which adaptively233 uses two different LS methods and came to the same conclusions.234 3. Region-based memetic algorithm with archive235 In this section we present the region-based MA with archive (RMAwA), an236 algorithm designed for multimodal optimisation which uses a niching technique237 to obtain as much optima as possible.238 7 ACCEPTED MANUSCRIPT ACCEPTED MANUSCRIPT RMAwA is a MA which alternatively applies an EA through a certain num-239 ber of evaluations and a LS method to the best solution in the population until240 stagnation. It then considers that an optimum has been reached, thus it stores241 that solution in an external archive and the EA is carried on.242 To maintain diversity during the search the algorithm divides each dimension243 in regions of same size, dividing the domain search in hypercubes. RMAwA uses244 these regions in two ways: First, only one solution is allowed in each region,245 thus when a solution generated by the EA falls in a region already occupied by246 a solution of the population the worst is removed. Second, regions in which one247 optimum has been found, by means of LS, are considered to be explored enough248 and discarded from the search space. The size of regions decreases during the249 run, by increasing the number of divisions per dimension.250 In order to efficiently discard regions from further exploration, this model251 maintains an index of the regions represented by a solution in the archive. Also,252 it stores all the found optima to recalculate the regions when its number changes.253 In the following subsections, we detail the algorithm. First, we briefly de-254 scribe the concept of the region-based niching strategy. Then, we explain the255 general scheme of the algorithm along with how the different components are256 integrated. Finally, we explain how the archive works in detail: its structure,257 which solutions are stored, and how it is used.258 3.1. Region-based niching strategy259 In [21], a novel niching strategy was proposed that redefines the notion260 of niche from the area surrounding each solution in the population to a fixed261 division of the search space. Each dimension of the search space is divided into a262 certain number of divisions, ND, creating a predefined grid of equal hypercubes263 representing the niches.264 In [56, 52], the authors use a similar partitioning of the search space to265 approximate the basin of attractions in multimodal fitness landscapes by means266 of clustered genetic search. In our algorithm, this fragmentation is used to define267 different niches in the search space. Ideally, regions contain a single basin of268 attraction but the unpredictability of the number of optima and their repartition269 in the search domain can not guarantee that. An illustration of the divisions of270 the search space can be seen in Figure 1. A solution sn∈RDis a real-parameter271 vector representing a solution to the problem at hand. It is associated with its272 region identified by its indices in each dimension, represented by a vector of273 integer values rn∈ND. The advantage of such definition is to allow faster274 retrieval of the existing niches by avoiding the computationally expensive cost275 of calculating the euclidean distance between solutions.276 In a region-based niching strategy, solutions generated in the evolutionary277 process compete with either the current solution present in the same region or278 the worst individual of the population. This technique can thus be assimilated279 to a clearing strategy in the sense that solutions compete to represent each niche280 in the population. The difference with classical niching strategy is the definition281 of the niche going from an euclidean distance-based representation to a region-282 based representation. In order to reduce the influence of the niche/region size, a283 8 ACCEPTED MANUSCRIPT ACCEPTED MANUSCRIPT Table 1: CEC’2013 benchmark problems Problem Function D Number of optima MaxF Es F1f11 2 5 ·104 F2f21 5 5 ·104 F3f31 1 5 ·104 F4f42 4 5 ·104 F5f52 2 5 ·104 F6f62 18 2 ·105 F7f72 36 2 ·105 F8f63 81 4 ·105 F9f73 216 4 ·105 F10 f82 12 2 ·105 F11 f92 6 2 ·105 F12 f10 2 8 2 ·105 F13 f11 2 6 2 ·105 F14 f11 3 6 4 ·105 F15 f12 3 8 4 ·105 F16 f11 5 6 4 ·105 F17 f12 5 8 4 ·105 F18 f11 10 6 4 ·105 F19 f12 10 6 4 ·105 F20 f12 20 8 4 ·105 over all global optima within the MaxF Es evaluations, and it is calculated by454 following Eq. 4:455 PR =PNR i=1 NPFi NKP ∗NR (4) where NPFiis the number of global optima found in the ith run, NKP is456 the number of known global optima and NR is the number of runs (for this457 benchmark, NR = 50). It is considered that an optimum optim is obtained if458 a solution sol is found where dist(sol, optim)≤, where dim is the Euclidean459 distance, and is a real value called accuracy level. The PR are calculated460 according to five different accuracy levels ={10−1, 10−2, 10−3, 10−4, 10−5}.461 Comparisons between algorithms have been performed for each accuracy462 level independently. For the comparison of two algorithms we considered non-463 parametric statistical tests [9]. More specifically, we used the Wilcoxon matched-464 pairs signed ranks tests for the direct comparison of two algorithms.465 4.3. Automatic configuration466 Setting the parameters of a new proposal can be a long and tedious task.467 Moreover, it does not ensure an optimal setting for these parameters. Consider-468 ing the novelty of certain components in this algorithm, it is more reliable to use469 an automatic configuration tool to assist in the design of the algorithm tuning470 the most critical parameters. To do so, we have used IRACE [29]. The IRACE471 package has already been extensively tested in several research projects, leading472 15 ACCEPTED MANUSCRIPT ACCEPTED MANUSCRIPT to successful improvement of the state-of-the-art, see for instance [28, 27]. The473 reader may refer to [41] for more information about IRACE and its parameters474 (we have used the recommended parameter values).475 We selected a set of parameters that we considered the most critical, and476 tuned them over the 20 problems of the CEC’2013 benchmark. For the non-477 tuned parameters we have selected commonly used values when not recom-478 mended values where given by from its authors. The list of tuned parameters479 can be seen in Table 2, showing for each parameter the explored range and the480 final value obtained by IRACE.481 Table 2: Tuned parameters and obtained values Parameters Descriptions Ranges Tuned iEA EA intensity, number of evaluations allocated to each EA application [100,1000] 550 iLS LS intensity, number of evaluations allocated to each LS application [100,1000] 150 ND0 Initial number of divisions, defines the size of the niches/regions [2,10] 2 uNumber of update to be performed [2,5] 4 muUpdate multiplier [1,5] 1.7 NP Population size of the EA [40,120] 70 αParameter for the BLX − αcrossover [0.1,0.9] 0.9 We can note that the EA intensity is almost four times the LS intensity.482 This is due to the fact that the LS is applied multiple times (until the improve-483 ments brought not significant enough) in each cycle. Concerning the number of484 division, we can see that the smallest number of divisions have been preferred485 (ND0= 2) along with a slow increase during the search by multiplying four486 times by 1.7: NDi+1 =ceil(1.7·NDi). The number of the divisions sequence487 is then [2,4,7,12,21]. Finally an important thing to note is the value of the α488 parameter for the BLX-α. Set to a high value (α= 0.9), it gives the EA a great489 exploration range.490 The other parameters listed in Table 3 were left to their default values taken491 from the corresponding papers. δmin LS defines the accuracy required for the search492 and is set to 10−6as the highest accuracy level required is 10−5. Concerning493 CMA-ES problems, we have set them to the default values as given in [17]. and494 the size of NAM selection method is taken from the previous work in [21].495 The parameters presented in Table 2 and 3 are the ones used in every ex-496 periment performed on every function and dimension of the benchmark.497 16 ACCEPTED MANUSCRIPT ACCEPTED MANUSCRIPT Table 3: Other parameters Parameters Descriptions Value λ Parameter to define the CMA-ES population size p= 4 + λln(D) 3 [17] µDefines the parent size for the CMA-ES p/µ 2 [17] NAMsize Size of the NAM selection method 3 [21] δmin LS Threshold for the LS stopping criterion 10−6 4.4. Possibility of finding all optima498 In this section, we discuss the ability of RMAwA to find all optima with an499 unlimited number evaluation. In other words, we wish to ensure that the search500 is not restricted to any subset of the whole search domain. For this model, we501 identify two phenomena that can cause such restriction and we discuss here if502 their occurrence is possible in the proposal.503 First, in population-based algorithm the risk of premature convergence of504 the population may lead to a genetic drift. The fact that RMAwA regularly505 generates new random solutions (when a solution is placed in the archive, it is506 replaced by a random solution) ensures sufficient diversity in the population to507 prevent premature convergence.508 The second risk that can be identified in this model is due to the restriction509 of the search to regions represented in the archive. Indeed, if a region repre-510 sented in the archive contains more that one optimum, some optima might be511 ignored. The probability of having more than one optimum present in the same512 region (noted M) is directly proportional to the hyper-volume of the regions Vr 513 calculated by Eq. 5:514 P(M) = a.Vr(5) where ais a variable that is dependent on the objective function fand the515 search domain. Basically, the smaller the region, the less probable that it will516 contain multiple optima. Thanks to the region size update, Vrkeeps decreasing517 during the search. In our algorithm, we make a limited number of reductions518 because the fitness evaluation number is also very limited. For an extremely519 large fitness evaluation number, the reductions would be applied repeatedly,520 reducing the hyper-volume of the regions each time. Thus, for an unlimited521 number of evaluations MaxF Es:522 lim MaxF Es→+∞Vr= 0 (6) 17 ACCEPTED MANUSCRIPT ACCEPTED MANUSCRIPT Hence:523 lim MaxF Es→+∞P(M) = 0 (7) Thus, there is no risk of limiting the search.524 5. Experimental results525 In this section, we are going to study the behaviour of the different compo-526 nents of our proposal, and we are going to compare our algorithm to previous527 algorithms in the literature. All the experiments are carried out following the528 experimental framework explained in previous section.529 The analysis of our proposal include the following experiments: First, we530 prove that using the region definition of a niche compared to the euclidean def-531 inition is more efficient in terms of computational time and exploration. Then,532 we demonstrate that using the solutions in the archive as excluding regions en-533 hance the performance of the model. We also analyse the influence of the region534 based niching strategy with the archive on the diversity of the population and in535 the exploration factor. Then, we analyse the memory and computational cost of536 the archive and the different components of the algorithm. Finally, we compare537 the proposed algorithm RMAwA with existing algorithms.538 5.1. Region niches versus classical niches539 Here, we assess the efficiency in terms of computation time and performance540 of the region definition of niches against the classical definition which implies541 calculating the euclidean distance between solutions. To do so, we consider the542 model presented without the use of the archive.543 The resulting algorithm here simply referred to as region based memetic544 algorithm (Region-MA) is opposed to an equivalent algorithm which uses the545 euclidean distance based definition of a niche as it is used in the classical clear-546 ing algorithm. This version is referred to as euclidean-distance based memetic547 algorithm (Euclidean-MA). On the generation of a new solution by the EA, the548 offspring created compete with the solutions falling within its niche radius σ,549 which is set to half the size of a region. In Region-MA, as it is explained in550 Section 3.1, new solutions created by the EA compete with the solutions already551 in the same regions.552 In order to simplify the display of the results, we will only focus on the553 highest level of accuracy (= 10−5). Indeed, the definition of a niche only554 affects the ability of the algorithm to explore the search space and not the555 precision of the solutions obtained.556 In Table 4, we show the PRs obtained by both versions along with the exe-557 cution time difference in percentage. We can see that the results of Region-MA558 are clearly better, and with significant differences (comparing with Wilcoxon’s559 test, the use of regions is statistically better with a p-value<0.001, see Table 5).560 Also, the execution time is much smaller, over the whole benchmark, using the561 region-based niches saves up to 17.4% of time.562 18 ACCEPTED MANUSCRIPT ACCEPTED MANUSCRIPT Table 4: PRs (for = 10−5) obtained by Region-MA and Euclidean-MA and execution time difference (in percentage) Problem F1F2F3F4F5 Region-MA 0.81 0.42 1 0.97 0.99 Euclidean-MA 0.77 0.56 1 0.36 0.87 Time difference (%) -35.88 -26.10 -28.36 -45.20 -43.57 Problem F6F7F8F9F10 Region-MA 0 0.7 0.06 0.22 0.94 Euclidean-MA 00.05 0.06 0.01 0.13 Time difference (%) -30.26 -39.05 -42.13 -38.96 -24.89 Problem F11 F12 F13 F14 F15 Region-MA 0.68 0.86 0.63 0.64 0.15 Euclidean-MA 0.27 0.14 0.2 0.18 0.14 Time difference (%) -19.42 -20.90 -28.38 -19.20 -21.11 Problem F16 F17 F18 F19 F20 Region-MA 0.36 0.16 0.17 0.13 0.13 Euclidean-MA 0.19 0.13 0.17 0.13 0.13 Time difference (%) -15.93 -1.56 -25.74 -21.19 -7.42 Table 5: Wilcoxon comparison of the P R obtained by Region-MA and Euclidean-MA (for = 10−5) R+ RRegion-MA Euclidean-MA p-value 189 21 0.0008 19 ACCEPTED MANUSCRIPT ACCEPTED MANUSCRIPT 5.2. Using the archive to reduce the search space563 The archive is used to store solutions considered as optima to allow the564 algorithm to remove them from the population without losing them. In our565 algorithm, it is used also to mark some regions as areas excluded for the search.566 In this section, we are interested in assessing how using the regions represented567 in the archive as excluded areas for the exploration of the EA improves the568 exploration of the search space and thus the discovery of more optima.569 In order to perform this comparison, we ran two versions of the algorithm.570 The first one is as presented in Section 3. The second one is the same algorithm571 without verifying that each solution created by the EA is present or not in572 the archive (steps 4-7 in Algorithm 2 are ignored). We thus compare here the573 proposed algorithm which uses an excluding archive (RMAwA) against one with574 a simple archive called RMA with Simple Archive (RMAwSA).575 As in the previous experiment, we will only focus on the highest level of576 accuracy (= 10−5). Indeed, the specific use of the archive mainly affects577 the algorithm’s ability to explore the search space and not the precision of the578 solutions obtained.579 Table 6: PRs of the RMA using an excluding archive (RMAwA) and a simple archive (RMAwSA) for = 10−5and computational time difference between the two versions. Problem F1F2F3F4F5 RMAwA 1.000 1.000 1.000 1.000 1.000 RMAwSA 1.000 0.312 1.000 1.000 1.000 Time difference (%) 22.6 23.3 7.7 15.5 3.1 Problem F6F7F8F9F10 RMAwA 0.000 0.917 0.824 0.513 1.000 RMAwSA 0.000 0.658 0.908 0.343 0.983 Time difference (%) 46.3 34.8 50.8 43.4 4.1 Problem F11 F12 F13 F14 F15 RMAwA 1.000 1.000 0.997 0.813 0.703 RMAwSA 0.667 0.930 0.667 0.667 0.648 Time difference (%) 5.8 1.5 2.2 21.8 15.3 Problem F16 F17 F18 F19 F20 RMAwA 0.670 0.660 0.233 0.128 0.125 RMAwSA 0.667 0.323 0.183 0.125 0.125 Time difference (%) 5.0 14.1 2.4 0.7 1.2 In Table 6, we show the PRs obtained by both versions of the algorithm and580 the time difference. Thanks to the excluding property of the archive, the per-581 formances of the algorithm are significantly improved (see Table 7 for Wilcoxon582 comparison). We also display in this table the CPU time increase caused by583 the use of the archive in the search. As we could have expected, this prop-584 erty implies more computational effort. However, the percentage increase in the585 20 ACCEPTED MANUSCRIPT ACCEPTED MANUSCRIPT computational time is reduced with the complexity and the dimensionality of586 the problem. This can be easily explained by the fact that in higher dimensions,587 the computational time of the evaluation increases while the time cost of the588 archive remains steady regardless the dimensionality. Also, considering the sum589 of the computational time for the whole benchmark, the runtime of RMAwA590 is 8.2% higher than RMAwSA’s (it cannot be calculated from table 6 because591 some functions take longer than others).592 Table 7: Wilcoxon comparison of the P R of the RMA with and without archive (for = 10−5) R+ RRMAwA RMAwSA p-value 186.5 23.5 0.00132 5.3. Diversity and Exploration593 In this section we analyse how RMAwA explores the search domain. First,594 we are going to study how the population diversity evolves along the search.595 Then, we visually analyse the exploration of the algorithm by plotting for several596 functions the solutions generated during the exploration phase.597 5.3.1. Population diversity: Influence of the Number of Divisions598 In this section we analyse the evolution of the population diversity during599 the search, and the influence of ND over the diversity. To do so, additional runs600 have been carried out and a diversity measure has been applied to the solutions601 into the population. The diversity measure applied is the following:602 DiversityP op =PN−1 i=1 PN j=i+1 dist(xi, xj) N·(N−1)/2(8) where Pop is the current population, Nis the population size, dist is the603 Euclidean distance, and xi, xjare solutions in the population.604 To study the influence of the current ND over the diversity, we are going605 to visualise and compare the diversity of the proposal (using the adaptive ND606 mechanism described in 3.2), with using a fixed ND.607 Figure 3 shows the evolution of the diversity for functions F7,F16 and F18.608 These functions have been selected for being representative of the different be-609 haviours detected in this benchmark. In axis xthere is the number of evalua-610 tions, and in axis ythe diversity measure. The vertical lines mark the updates611 of number of divisions (it only has influence over the adaptive ND version),612 dividing the axis xin five stages of the algorithm (each stage using a different613 ND). In the following, we are going to describe the main tendencies:614 1. In functions with a small dimension, like F7(where D=2), we can observe615 two phases. In the initial stages of the search (ND=2, ND=4), because616 21 ACCEPTED MANUSCRIPT ACCEPTED MANUSCRIPT (a) F7 (b) F16 (c) F18 Figure 3: Diversity of the RMAwA population using adaptive number of divisions and using different fixed number of divisions during one run 22 ACCEPTED MANUSCRIPT ACCEPTED MANUSCRIPT there are few regions, when a region is avoided the search space is reduced617 very quickly to a small portion of the whole space, thus the diversity618 decreases very quickly. Indeed, for these ND values the fixed ND version619 prematurely stops because all possible regions have a local optima. The620 subsequent updates in ND increase the number of regions, releasing space621 for the EA to explore and thus increasing the diversity. As compared with622 fixed ND, the diversity of adaptive ND is very similar in the first two and623 final stages, with a greater diversity in the stages inbetween.624 2. In functions with medium dimensionality, like F16 (D=5), the same phe-625 nomena is observed. However, after reaching a certain number of divisions626 per dimensions (third update) the diversity decreases, because the algo-627 rithm allows solutions more closer between them, reducing the diversity628 to enforce the exploitation of found solutions. Comparing adaptive ND629 with fixed ND, we can observe that diversity adaptive ND is actually very630 similar to ND in each stage.631 3. In functions with higher dimensionality, like F18 (D=10), we can see that632 the diversity constantly decreases at each increase of the number of divi-633 sions. In these functions, it seems that the niching model does not provide634 a good balance in the population diversity during the search. Comparing635 adaptive ND with fixed ND, we can observe that adaptive ND obtains636 very close results to obtained by the fixed ND in each stage.637 The previous section has shown the diversity differences comparing several638 fixed ND and the proposed dynamic ND. However, diversity itself is not our639 goal, thus we are going to compare the obtained PRs for each case. Table 8640 show the results, highlighting the results for those functions whose diversity has641 been analysed. We can observe that:642 •In functions with a small dimension, like F7, in which a higher ND implies643 a better diversity, the number of optima increases also with the ND. Better644 results are obtained with dynamic ND.645 •In functions with medium dimensionality, like F16, in which for certain646 ND values the diversity is reduced very quickly, the PR decreases when647 ND increases. Dynamic ND, on the contrary, obtains the best PR value.648 •Results obtained in functions with higher dimensionality, like F18, proves649 that there is noy a good balance in the diversity, and that it has bad650 consecuences for the obtained PR. In this case, dynamic ND obtains worse651 results than using ND=2 but better than the other values.652 In summary, Figure 3 shows that the number of regions and problem di-653 mensionality have strong influence over the diversity in the population and the654 number of found optima, and that diversity using an adaptive ND is very close655 to that obtained with a fixed ND in each stage, obtaining the most robust656 behaviour when finding the optima.657 23 ACCEPTED MANUSCRIPT ACCEPTED MANUSCRIPT Table 8: PRs of the RMAwA using different fixed numbers of divisions (ND) and with dynamic ND. Dynamic Function ND=2 ND=4 ND=7 ND=12 ND=21 ND F1 0.900 1.000 1.000 1.000 1.000 1.000 F2 1.000 1.000 1.000 1.000 1.000 1.000 F3 1.000 1.000 1.000 1.000 1.000 1.000 F4 0.750 1.000 1.000 1.000 1.000 1.000 F5 1.000 1.000 1.000 1.000 1.000 1.000 F6 0.000 0.000 0.000 0.000 0.000 0.000 F7 0.084 0.429 0.612 0.790 0.829 0.917 F8 0.023 0.290 0.610 0.458 0.853 0.824 F9 0.035 0.172 0.433 0.660 0.618 0.513 F10 0.923 1.000 1.000 1.000 1.000 1.000 F11 0.733 1.000 1.000 1.000 1.000 1.000 F12 0.470 0.840 0.875 0.955 1.000 1.000 F13 0.680 0.993 1.000 1.000 1.000 0.997 F14 0.760 0.813 0.940 0.727 0.647 0.813 F15 0.665 0.725 0.675 0.640 0.275 0.703 F16 0.667 0.667 0.533 0.300 0.273 0.670 F17 0.660 0.680 0.250 0.185 0.165 0.660 F18 0.473 0.167 0.167 0.167 0.167 0.233 F19 0.160 0.125 0.125 0.125 0.125 0.128 F20 0.125 0.125 0.125 0.125 0.125 0.125 Mean 0.555 0.651 0.667 0.657 0.654 0.729 24 ACCEPTED MANUSCRIPT ACCEPTED MANUSCRIPT As expected, we can observe a strong increase of the physical size used726 by the archive for the most complex problems. However, the memory used727 remains reasonable for today’s machines. In the most extreme problem, F20 728 where D= 20, the archive only uses 64.88 kB of memory. Even if it might appear729 irrelevant for such problems, the size of the archive can increase exponetially730 with the dimensionality and the number of evaluation. When tackling large731 scale problems, one may consider limiting the size of the archive.732 Table 9: Average number of elements in the archive’s list (|S|), the index (|R|) and total memory used by the archive (in kB) at the end of each run Problem D MaxF Es |S| |R|ArchiveSize F11 5.00 ·104135.92 4.58 1.08 F21 5.00 ·104130.24 9.96 1.06 F31 5.00 ·104129.32 10.52 1.05 F42 5.00 ·104106 22.76 1.83 F52 5.00 ·104112.76 14.5 1.88 F62 2.00 ·105425.52 112.64 7.53 F72 2.00 ·105448.28 100.18 7.79 F83 4.00 ·105681.84 398.62 20.65 F93 4.00 ·105811.64 389.08 23.58 F10 2 2.00 ·105431.28 100.68 7.53 F11 2 2.00 ·105372.72 106.42 6.66 F12 2 2.00 ·105326.04 104.42 5.91 F13 2 2.00 ·105349.52 121.84 6.41 F14 3 4.00 ·105583 283.48 16.99 F15 3 4.00 ·105581.6 278.68 16.90 F16 5 4.00 ·105524 259.42 25.54 F17 5 4.00 ·105516.64 270.26 25.46 F18 10 4.00 ·105446.84 187.36 42.23 F19 10 4.00 ·105338.52 168.28 33.02 F20 20 4.00 ·105343.8 142.92 64.88 5.4.2. Computational time of the different components of RMAwA733 In this section, we analyse the amount of time taken by the different com-734 ponents of RMAwA over a whole run, namely:735 •LS operations: the operations performed by CMA-ES during its search736 process.737 •EA operations: the operations performed by the SSGA to evolve the pop-738 ulation.739 •Niching: the time it takes for a new solution to go through the niching740 process (retrieval and comparison of the solutions present in the same741 31 ACCEPTED MANUSCRIPT ACCEPTED MANUSCRIPT region in the population).742 •Archive: the time implied by the excluding property of the archive (as-743 sessing the presence of the solution’s region in the archive’s index).744 First, to assess the computational time of each component, we use function745 f12. This function presents the advantage of being implemented in 4 dimensions,746 D={3,5,10,20}, allowing us to evaluate the scalability of the proposal. For747 those four problems, we calculate the CPU time used by each component to748 assess their scalability. The search effort is unequally divided between the LS749 and the EA (the number of evaluation at each EA application is fixed while the750 number of evaluation for each LS application is not limited). Thus, to perform751 a fair comparison, we only select the average time per evaluation. We plot the752 results in Figure 10.753 As far as we can see, the complexity of the niching strategy and the use of754 the archive are barely affected by an increase of the dimensionality. In the same755 way, the operations of the SSGA algorithms show interesting scalable properties.756 The main weakness lies in the use of CMA-ES as LS method. Although it offers757 a low complexity in the lowest dimensions, with more than ten variables, CMA-758 ES shows poor scalability in terms of complexity.759 In order to counterbalance the importance of this drawback, we show in ta-760 ble 10 the CPU time of each of the components along with the evaluation time.761 Here, we remind the reader of the notation used in this paper, we grouped the762 problems Fjby function fiin order to make for easier reading and see the rela-763 tions between the different dimensions of each function. From this table, when764 increasing the dimensionality, even if the proportion of the LS (i.e. CMA-ES)765 operations increases, the total CPU time is particularly affected by the compu-766 tational time of the evaluation which is independent of the algorithm. However,767 as the complexity of CMA-ES increases exponentially with the dimension, larger768 scale problems may require the use of another LS method.769 Sheet2 Page 14 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 0.0E+0 1.0E-2 2.0E-2 3.0E-2 4.0E-2 EA LS Archive Niching Dimension Time (in ms) Figure 10: CPU time (in ms) of each component per evaluations for problem f12 for different dimensions 32 ACCEPTED MANUSCRIPT ACCEPTED MANUSCRIPT Table 10: CPU time (in seconds) details of RMAwA for each problem Fj={fi, D}with the percentage in the whole optimisatioc process Problem F1={f1,1}F2={f2,1}F3={f3,1}F4={f4,2} Archive 0.150 (18.41%) 0.355 (40.84%) 0.148 (18.99%) 0.099 (14.96%) Niching 0.236 (28.91%) 0.221 (25.39%) 0.254 (32.43%) 0.180 (27.17%) EA 0.260 (31.86%) 0.257 (29.59%) 0.282 (36.08%) 0.207 (31.18%) LS 0.164 (20.08%) 0.027 (3.07%) 0.064 (8.23%) 0.172 (25.93%) Evaluations 0.006 (0.73%) 0.010 (1.11%) 0.033 (4.27%) 0.005 (0.76%) Total 0.816 0.869 0.782 0.664 Problem F5={f5,2}F6={f6,2}F8={f6,3}F7={f7,2} Archive 0.127 (9.96%) 1.085 (35.82%) 1.512 (29.74%) 1.354 (40.56%) Niching 0.204 (16.01%) 0.759 (25.05%) 1.188 (23.37%) 0.787 (23.58%) EA 0.233 (18.27%) 0.884 (29.18%) 1.513 (29.75%) 0.942 (28.23%) LS 0.696 (54.67%) 0.159 (5.24%) 0.494 (9.72%) 0.173 (5.18%) Evaluations 0.014 (1.10%) 0.143 (4.70%) 0.377 (7.41%) 0.081 (2.44%) Total 1.273 3.029 5.085 3.338 Problem F9={f7,3}F10 ={f8,2}F11 ={f9,2}F12 ={f10,2} Archive 2.210 (35.27%) 0.740 (28.09%) 0.615 (6.15%) 0.521 (5.38%) Niching 1.403 (22.38%) 0.759 (28.80%) 0.675 (6.74%) 0.586 (6.05%) EA 1.793 (28.61%) 0.934 (35.41%) 0.903 (9.02%) 0.802 (8.28%) LS 0.661 (10.55%) 0.153 (5.82%) 0.332 (3.32%) 0.207 (2.13%) Evaluations 0.200 (3.19%) 0.050 (1.88%) 7.479 (74.77%) 7.571 (78.16%) Total 6.267 2.636 10.003 9.687 Problem F13 ={f11,2}F14 ={f11 ,3}F16 ={f11 ,5}F18 ={f11,10} Archive 0.604 (6.15%) 0.910 (3.88%) 0.853 (2.60%) 1.065 (1.81%) Niching 0.630 (6.42%) 1.028 (4.39%) 0.951 (2.90%) 0.868 (1.48%) EA 0.852 (8.69%) 1.441 (6.15%) 1.300 (3.96%) 1.267 (2.16%) LS 0.210 (2.14%) 0.412 (1.76%) 0.671 (2.04%) 2.099 (3.58%) Evaluations 7.515 (76.59%) 19.641 (83.82%) 29.054 (88.50%) 53.400 (90.97%) Total 9.812 23.431 32.829 58.698 Problem F15 ={f12,3}F17 ={f12 ,5}F19 ={f12 ,10}F20 ={f12,20} Archive 0.982 (4.15%) 0.773 (2.33%) 0.723 (1.22%) 0.891 (0.71%) Niching 1.020 (4.31%) 0.932 (2.81%) 0.663 (1.12%) 0.714 (0.57%) EA 1.387 (5.86%) 1.282 (3.87%) 1.144 (1.93%) 1.435 (1.15%) LS 0.399 (1.69%) 0.580 (1.75%) 2.267 (3.81%) 12.974 (10.37%) Evaluations 19.867 (83.98%) 29.560 (89.23%) 54.636 (91.93%) 109.063 (87.20%) Total 23.655 33.127 59.434 125.077 33 ACCEPTED MANUSCRIPT ACCEPTED MANUSCRIPT 5.5. Comparison with existing algorithms770 In this section we compare the results obtained by our algorithm, RMAwA.771 We selected a number of algorithms from the literature along with algorithms772 presented for the CEC’2013 competition:773 •PNA-NSGAII [2] proposed for the competition, this algorithm is an im-774 provement of A-NSGAII [8]. These algorithms tackle the multimodal op-775 timisation problem by turning them into bi-objective problems. The first776 objective is the minimisation of the original function and the second one777 is the maximisation of the diversity brought by the evaluated individual.778 •dADE/nrand/1/bin [12] : a DE using a neighbourhood based mutation779 strategy and a dynamically updated archive.780 •DE/nrand/2 [13] : a DE using the neighbourhood based mutation strategy.781 •NVMO [34]: a Variable Mesh optimisation algorithm with niching strat-782 egy.783 •CMA-ES [18]: A version of CMA-ES that implements a simple archive.784 •NEA2 [43]: A version of CMA-ES that uses nearest-better clustering as785 niching strategy.786 These algorithms are the top six algorithms of the CEC’2013 competition.787 All the results used here were provided by the authors and used during the788 competition. The detailed results of each algorithm can be seen in the Appendix.789 We first analyse the overall performance of each algorithm on the benchmark790 and compare them with RMAwA. Then we study in detail their behaviour791 according to the problem’s characteristics.792 5.5.1. Accuracy level analysis793 We analyse here the general performance of these algorithms on the CEC’2013794 benchmark for each accuracy level. To support this analysis, we show in Ta-795 ble 11 the mean rankings of each algorithm according to the different accuracy796 levels and in Table 12 the Wilcoxon comparison of RMAwA with the other797 algorithms.798 First, when comparing with other algorithms using CMA-ES, we can see that799 RMAwA significantly outperforms the classical CMA-ES (with α= 0.1). This800 algorithm is not particularly designed for multimodal optimisation as it does801 not implement any niching mechanism. When comparing with NEA2, RMAwA802 offers similar performance.803 Then, we can see that RMAwA is third best for the smallest accuracy level804 (= 1E−1) behind NVMO and dADE although no statistical difference can be805 observed in Table 12.806 Excluding NEA2, RMAwA obtains better results than the other algorithms,807 and this superiority increases with the accuracy level, being specially remarkable808 34 ACCEPTED MANUSCRIPT ACCEPTED MANUSCRIPT Table 11: Mean rankings obtained by different algorithms over all functions CEC’2013 benchmark for each accuracy level Accuracy level 1E−1 1E−2 1E−3 1E−4 1E−5 PNA-NSGAII 4.53 5.18 5.28 5.45 5.43 DE/nrand/2 5.53 5.05 4.95 4.83 4.30 CMA-ES 4.58 4.00 4.08 3.93 3.90 NVMO 2.73 3.68 4.15 4.43 5.08 dADE 3.43 4.10 4.10 4.08 3.90 NEA2 3.70 2.75 2.55 2.68 2.70 RMAwA 3.53 3.25 2.90 2.63 2.70 for ={1E−4,1E−5}. Between NEA2 and RMAwA there is no statistical809 difference detected.810 This analysis highlights the difficulty of algorithms to properly balance the811 exploration and the exploitation. Indeed, when algorithms use the original812 CMA-ES, a very efficient method to obtain accurate solutions but also very813 costly, they generally perform better for higher accuracy levels. On the other814 hand, other algorithms (DE-based, NVMO, PNA-NSGAII) have better explo-815 ration efficiency but fail to identify accurate solutions.816 5.5.2. Problem specific performance analysis817 Let us now consider every problem individually. As it is the most challenging818 for this benchmark, we will consider here only the highest accuracy level (=819 1E−5). Table 13 lists the PRs obtained by each algorithm for this accuracy820 level.821 In this analysis we will focus on the problems offering the major differences822 between the results obtained by the compared algorithms. Concerning problems823 with highly multimodal fitness landscapes, F7to F9where the number of optima824 ranges from 36 to 216, RMAwA ranks amongst the best algorithm. It obtains825 the best results for problem F7and obtains the second best results of problem826 F8and F9after respectively dADE and NEA2.827 RMAwA also shows the best results in the problems with composition func-828 tions (F10 to F17). However the quality of the results decreases with the dimen-829 sionality (F18 to F20), being clearly worse than that obtained by NEA2.830 The improvable behaviour of RMAwA when the dimensionality increases is831 clear because for f11 and f12 results are very good with dimension 2, but not832 good with a higher dimension, like 10 (F18-F20). In Section 5.3.1, it can be833 observed that for these functions the results obtained are low for each possible834 ND, thus the results are not due to the ND adaptation mechanism. Another835 possible reason of the improvable results could be that the parameter values836 have been automatically tuned considering all functions, in which the majority837 has a very low dimension. Because of this, these parameter values could not be838 the more adequate for lower dimension problems. In order to reject or confirm839 35 ACCEPTED MANUSCRIPT ACCEPTED MANUSCRIPT Table 12: Wilcoxon comparison of the PRs of the RMAwA (R+) with other algorithms (R−) (for ={1E−1,1E−2,1E−3,1E−4,1E−5}) = 1E−1 RMAwA vs R+R−p-value PNA-NSGAII 128.5 65 2.27E-1 DE/nrand/2 180.5 29.5 3.40E-3 CMA-ES 168.5 41.5 1.62E-2 NVMO 62 132.5 1.96E-1 dADE 71.5 122 3.44E-1 NEA2 117.5 92.5 6.41E-1 = 1E−2 RMAwA vs R+R−p-value PNA-NSGAII 199.5 10.5 9.35E-5 DE/nrand/2 171.5 38.5 1.14E-2 CMA-ES 150.5 59.5 9.35E-2 NVMO 124.5 85.5 4.67E-1 dADE 134.5 75.5 2.71E-1 NEA2 82.5 127.5 4.01E-1 = 1E−3 RMAwA vs R+R−p-value PNA-NSGAII 185 7.5 8.39E-5 DE/nrand/2 171.5 38.5 1.14E-2 CMA-ES 139 53.5 9.98E-2 NVMO 166.5 43.5 2.04E-2 dADE 147.5 62.5 1.19E-1 NEA2 91.5 118.5 6.14E-1 = 1E−4 RMAwA vs R+R−p-value PNA-NSGAII 185 7.5 8.39E-5 DE/nrand/2 171.5 38.5 1.14E-2 CMA-ES 139 53.5 9.98E-2 NVMO 185.5 24.5 1.56E-3 dADE 165.5 44.5 2.27E-2 NEA2 95 97.5 1.00E+0 = 1E−5 RMAwA vs R+R−p-value PNA-NSGAII 199.5 10.5 9.35E-5 DE/nrand/2 151.5 42 3.23E-2 CMA-ES 138 54.5 1.09E-1 NVMO 189.5 20.5 7.79E-4 dADE 151.5 42 3.23E-2 NEA2 96.5 96 9.68E-1 36 ACCEPTED MANUSCRIPT ACCEPTED MANUSCRIPT Table 13: PRs obtained by each algorithm for = 1E−5 on the CEC’2013 benchmark. Values in the parenthesis represent the standard competition ranking of each algorithm for each problem Problem PNA-NSGAII DE/nrand/2 CMA-ES NVMO dADE NEA2 RMAwA F11.000 (1) 1.000 (1) 1.000 (1) 1.000 (1) 1.000 (1) 1.000 (1) 1.000 (1) F21.000 (1) 1.000 (1) 1.000 (1) 1.000 (1) 1.000 (1) 1.000 (1) 1.000 (1) F31.000 (1) 1.000 (1) 1.000 (1) 1.000 (1) 1.000 (1) 1.000 (1) 1.000 (1) F40.805 (7) 1.000 (1) 0.990 (5) 1.000 (1) 1.000 (1) 0.990 (5) 1.000 (1) F51.000 (1) 1.000 (1) 1.000 (1) 1.000 (1) 1.000 (1) 1.000 (1) 1.000 (1) F60.000 (1) 0.000 (1) 0.000 (1) 0.000 (1) 0.000 (1) 0.000 (1) 0.000 (1) F70.683 (5) 0.275 (7) 0.516 (6) 0.804 (3) 0.714 (4) 0.911 (2) 0.917 (1) F80.252 (4) 0.363 (3) 0.115 (6) 0.027 (7) 0.947 (1) 0.239 (5) 0.824 (2) F90.276 (4) 0.065 (7) 0.272 (5) 0.194 (6) 0.349 (3) 0.579 (1) 0.513 (2) F10 1.000 (1) 1.000 (1) 0.978 (6) 0.967 (7) 1.000 (1) 0.980 (5) 1.000 (1) F11 0.663 (7) 0.667 (5) 0.953 (3) 0.667 (4) 0.667 (5) 0.960 (2) 1.000 (1) F12 0.573 (7) 0.618 (5) 0.760 (3) 0.593 (6) 0.735 (4) 0.833 (2) 1.000 (1) F13 0.623 (7) 0.667 (4) 0.947 (2) 0.663 (6) 0.667 (4) 0.947 (2) 0.997 (1) F14 0.610 (7) 0.667 (4) 0.743 (3) 0.627 (6) 0.667 (4) 0.800 (2) 0.813 (1) F15 0.443 (5) 0.400 (6) 0.653 (3) 0.378 (7) 0.620 (4) 0.713 (1) 0.703 (2) F16 0.323 (7) 0.667 (3) 0.663 (5) 0.653 (6) 0.667 (3) 0.673 (1) 0.670 (2) F17 0.245 (7) 0.280 (6) 0.583 (3) 0.325 (5) 0.358 (4) 0.695 (1) 0.660 (2) F18 0.093 (7) 0.507 (3) 0.340 (4) 0.327 (5) 0.603 (2) 0.663 (1) 0.233 (6) F19 0.010 (6) 0.180 (3) 0.597 (2) 0.093 (5) 0.000 (7) 0.667 (1) 0.128 (4) F20 0.000 (5) 0.230 (3) 0.425 (1) 0.000 (5) 0.000 (5) 0.350 (2) 0.125 (4) 37 ACCEPTED MANUSCRIPT ACCEPTED MANUSCRIPT that hypothesis, we have carried out another automatic tuning considering only840 functions F18 −F20, but the results obtained were very similar. Thus, the841 improvable behaviour or RMAwA in higher dimensionality problems is kept as842 a open issue to be solved in the future.843 However, this previous behaviour is not unsurprising, because it has hap-844 pened to many others, as can be observed in Table 13. As was formulated by845 the No Free Lunch Theorem, designing an algorithm for an heterogeneous test846 bed of problems is very challenging. It is common for algorithms to perform847 well in problems with certain characteristics and poorly on others.848 In summary, the algorithm proposed of this paper, RMAwA, offers an overall849 performance significantly superior to the other algorithms by obtaining competi-850 tive if not better results in most problems (except in higher dimension problems)851 proposed in the CEC’2013 benchmark. Only NEA2, the winner of the CEC’2013852 competition obtains equivalent results.853 6. Conclusions854 In this paper, we present a novel model based on region-based MA to tackle855 multimodal optimisation problems. It uses a clearing strategy where niches are856 defined as regions. It implements an archive of solutions and indexed regions857 considered as explored and thus excluded from further exploration.858 In order to asses the efficiency of the model against existing ones, we have859 tested it on a MA which alternatively applies an EA (SSGA) to explore the860 search space and an LS (CMA-ES) to the best one until it does not improve for861 a certain number of evaluations.862 Various studies have been performed to study the performances of this863 model. First, we have demonstrated that the use of region-based niches was864 more efficient than that of the classical euclidean niches. We have shown that865 excluding regions explored by the LS allows the algorithm to reduce the search866 space leading to a more efficient exploration. Also, we have analysed the popu-867 lation diversity during the run and the degree of exploration in several functions.868 Finally, complexity testing show the good scalability of the proposal.869 We compared the resulting algorithm using the benchmark issued for the870 special session and competition on niching methods for multimodal function871 optimisation of the IEEE Congress on Evolutionary Computation in 2013. We872 noted that our algorithm was fairly independent to the different accuracy levels873 tested in this benchmark compared to the other algorithms obtaining signifi-874 cantly better results than most algorithms and similar performance to NEA2.875 This work opens the way of various potential future studies:876 •First, the behaviour of RMAwA with higher dimension problems should877 be studied more in detailed and improved it. Also, the memory cost of the878 archive may become consequent when tackling higher dimension problems879 and it could be interesting to study techniques which reduce or limit the880 size of the archive or remove similar solutions representing the same region.881 38 ACCEPTED MANUSCRIPT ACCEPTED MANUSCRIPT •As is often the case when the parameter defining the size of a niche, the882 number of divisions per dimensions is highly sensitive. Although the idea883 of a constant increasing during the search might offer interesting results,884 it may not be optimal in some cases. A more adapted strategy could be885 identified and researched. An other option would be to implement multi-886 population where, as it is done in Hierarchical Genetic Strategy [53], each887 population uses different numbers of divisions.888 •This model has proved to obtain interesting results when applied with889 both SSGA and CMAES as its EA and LS. Further experiments using890 different components or in a different memetic framework could lead to891 new efficient algorithms. If tested in higher dimensions, the necessity of892 changing CMAES to a more scalable LS method would be compulsory.893 Acknowledgements894 This Work was supported by the Spanish Ministry of Education and Science895 under Grants TIN2012-37930-C02-01, TIN2014-57251-P, TIN2013-47210-P and896 Research Regional Projects P10-TIC-6858, P12-TIC-2958.897 References898 [1] C. F. Agostinho, C. Fernandes, A. Rosa, A Study on Non-random Mating899 and Varying Population Size in Genetic Algorithms Using a Royal Road900 Function, in: Proceedings of the 2001 IEEE Congress on Evolutionary901 Computation, IEEE Press, 2001, pp. 60–66.902 [2] S. Bandaru, K. Deb, A parameterless-niching-assisted bi-objective ap-903 proach to multimodal optimization, in: Evolutionary Computation (CEC),904 2013 IEEE Congress on, 2013, pp. 95–102.905 [3] S. Bird, X. Li, Adaptively choosing niching parameters in a PSO, in: Pro-906 ceedings of the 8th annual conference on Genetic and evolutionary compu-907 tation, GECCO ’06, ACM, New York, NY, USA, 2006, pp. 3–10.908 [4] S. Biswas, S. Kundu, S. Das, An improved parent-centric mutation with909 normalized neighborhoods for inducing niching behavior in differential evo-910 lution, Cybernetics, IEEE Transactions on 44 (10) (2014) 1726–1737.911 [5] X. Chen, Y.-S. Ong, M.-H. Lim, K. C. Tan, A multi-facet survey on912 memetic computation, Evolutionary Computation, IEEE Transactions on913 15 (5) (2011) 591–607.914 [6] S. Das, S. Maity, B.-Y. Qu, P. Suganthan, Real-parameter evolutionary915 multimodal optimization a survey of the state-of-the-art, Swarm and Evo-916 lutionary Computation 1 (2) (2011) 71 – 88.917 39 ACCEPTED MANUSCRIPT ACCEPTED MANUSCRIPT [7] K. A. De Jong, An analysis of the behavior of a class of genetic adaptive918 systems., Ph.D. thesis, Ann Arbor, MI, USA (1975).919 [8] K. Deb, A. Saha, Multimodal optimization using a bi-objective evolution-920 ary algorithm, Evolutionary Computation 20 (1) (2012) 27–62.921 [9] J. Derrac, S. Garc´ıa, D. Molina, F. Herrera, A practical tutorial on the use922 of nonparametric statistical tests as a methodology for comparing evolu-923 tionary and swarm intelligence algorithms, Swarm and Evolutionary Com-924 putation 1 (1) (2011) 3–18.925 [10] G. Dick, P. Whigham, Weighted local sharing and local clearing for multi-926 modal optimisation, Soft Computing 15 (9) (2011) 1707–1721.927 [11] M. Ellabaan, Y.-S. Ong, Valley-adaptive clearing scheme for multimodal928 optimization evolutionary search, in: Intelligent Systems Design and Ap-929 plications, 2009. ISDA ’09. Ninth International Conference on, 2009, pp.930 1–6.931 [12] M. Epitropakis, X. Li, E. Burke, A dynamic archive niching differential932 evolution algorithm for multimodal optimization, in: Evolutionary Com-933 putation (CEC), 2013 IEEE Congress on, 2013, pp. 79–86.934 [13] M. Epitropakis, V. Plagianakos, M. Vrahatis, Finding multiple global op-935 tima exploiting differential evolution’s niching capability, in: Differential936 Evolution (SDE), 2011 IEEE Symposium on, 2011, pp. 1–8.937 [14] L. J. Eshelman, J. D. Schaffer, Real-coded Genetic Algorithms in Ge-938 netic Algorithms by Preventing Incest, Foundation of Genetic Algorithms939 2 (1993) 187–202.940 [15] W. Gao, G. G. Yen, S. Liu, A cluster-based differential evolution with self-941 adaptive strategy for multimodal optimization, Cybernetics, IEEE Trans-942 actions on 44 (8) (2014) 1314–1327.943 [16] D. E. Goldberg, J. Richardson, Genetic algorithms with sharing for multi-944 modal function optimization, in: Proceedings of the Second International945 Conference on Genetic Algorithms on Genetic algorithms and their appli-946 cation, L. Erlbaum Associates Inc., Hillsdale, NJ, USA, 1987, pp. 41–49.947 [17] N. Hansen, S. D. M¨uller, P. Koumoutsakos, Reducing the Time Complexity948 of the Derandomized Evolution Strategy with Covariance Matrix Adapta-949 tion (CMA-ES), Evolutionary Computation 1 (11) (2003) 1–18.950 [18] N. Hansen, A. Ostermeier, Completely derandomized self-adaptation in951 evolution strategies, Evolutionary Computation 9 (2) (2001) 159–195.952 [19] A. Kononova, D. Ingham, M. Pourkashanian, Simple scheduled memetic953 algorithm for inverse problems in higher dimensions: Application to chemi-954 cal kinetics, in: Evolutionary Computation, 2008. CEC 2008. (IEEE World955 40 ACCEPTED MANUSCRIPT ACCEPTED MANUSCRIPT Table 16: Results with DE/nrand/2 Pb Fun Dim Accuracy level 1E−1 1E−2 1E−3 1E−4 1E−5 F1f11 1.000 1.000 1.000 1.000 1.000 F2f21 1.000 1.000 1.000 1.000 1.000 F3f31 1.000 1.000 1.000 1.000 1.000 F4f42 1.000 1.000 1.000 1.000 1.000 F5f52 1.000 1.000 1.000 1.000 1.000 F6f62 0.669 0.669 0.669 0.669 0.000 F7f72 0.276 0.276 0.276 0.276 0.275 F8f63 0.365 0.365 0.365 0.365 0.363 F9f73 0.066 0.066 0.066 0.066 0.065 F10 f82 1.000 1.000 1.000 1.000 1.000 F11 f92 0.667 0.667 0.667 0.667 0.667 F12 f10 2 0.635 0.628 0.628 0.618 0.618 F13 f11 2 0.667 0.667 0.667 0.667 0.667 F14 f11 3 0.667 0.667 0.667 0.667 0.667 F15 f12 3 0.413 0.408 0.405 0.400 0.400 F16 f11 5 0.667 0.667 0.667 0.667 0.667 F17 f12 5 0.288 0.283 0.283 0.280 0.280 F18 f11 10 0.517 0.513 0.507 0.507 0.507 F19 f12 10 0.230 0.218 0.203 0.190 0.180 F20 f12 20 0.230 0.230 0.230 0.230 0.230 47 ACCEPTED MANUSCRIPT ACCEPTED MANUSCRIPT Table 17: Results with dADE/nrand/1/bin Pb Fun Dim Accuracy level 1E−1 1E−2 1E−3 1E−4 1E−5 F1f11 1.000 1.000 1.000 1.000 1.000 F2f21 1.000 1.000 1.000 1.000 1.000 F3f31 1.000 1.000 1.000 1.000 1.000 F4f42 1.000 1.000 1.000 1.000 1.000 F5f52 1.000 1.000 1.000 1.000 1.000 F6f62 1.000 1.000 1.000 0.988 0.000 F7f72 1.000 0.960 0.878 0.808 0.714 F8f63 0.990 0.991 0.985 0.958 0.947 F9f73 0.829 0.592 0.552 0.436 0.349 F10 f82 1.000 1.000 1.000 1.000 1.000 F11 f92 0.867 0.667 0.667 0.667 0.667 F12 f10 2 0.750 0.748 0.738 0.740 0.735 F13 f11 2 0.737 0.667 0.667 0.667 0.667 F14 f11 3 0.943 0.667 0.667 0.667 0.667 F15 f12 3 1.000 0.643 0.623 0.600 0.620 F16 f11 5 0.890 0.667 0.667 0.667 0.667 F17 f12 5 0.963 0.480 0.420 0.400 0.358 F18 f11 10 0.663 0.630 0.630 0.613 0.603 F19 f12 10 0.495 0.118 0.080 0.020 0.000 F20 f12 20 0.080 0.005 0.000 0.000 0.000 48 ACCEPTED MANUSCRIPT ACCEPTED MANUSCRIPT Table 18: Results with PNA-NSGA Pb Fun Dim Accuracy level 1E−1 1E−2 1E−3 1E−4 1E−5 F1f11 1.000 1.000 1.000 1.000 1.000 F2f21 1.000 1.000 1.000 1.000 1.000 F3f31 1.000 1.000 1.000 1.000 1.000 F4f42 1.000 1.000 0.995 0.985 0.805 F5f52 1.000 1.000 1.000 1.000 1.000 F6f62 0.562 0.536 0.523 0.473 0.000 F7f72 1.000 0.741 0.726 0.709 0.683 F8f63 0.352 0.330 0.310 0.275 0.252 F9f73 0.480 0.326 0.318 0.298 0.276 F10 f82 1.000 1.000 1.000 1.000 1.000 F11 f92 0.877 0.677 0.670 0.680 0.663 F12 f10 2 0.752 0.715 0.672 0.642 0.573 F13 f11 2 0.697 0.667 0.667 0.663 0.623 F14 f11 3 0.933 0.667 0.667 0.663 0.610 F15 f12 3 0.665 0.495 0.485 0.470 0.443 F16 f11 5 1.000 0.523 0.523 0.417 0.323 F17 f12 5 0.917 0.347 0.338 0.300 0.245 F18 f11 10 0.640 0.117 0.113 0.110 0.093 F19 f12 10 0.020 0.020 0.043 0.017 0.010 49 ACCEPTED MANUSCRIPT ACCEPTED MANUSCRIPT Table 19: Results with NVMO Pb Fun Dim Accuracy level 1E−1 1E−2 1E−3 1E−4 1E−5 F1f11 1.000 1.000 1.000 1.000 1.000 F2f21 1.000 1.000 1.000 1.000 1.000 F3f31 1.000 1.000 1.000 1.000 1.000 F4f42 1.000 1.000 1.000 1.000 1.000 F5f52 1.000 1.000 1.000 1.000 1.000 F6f62 1.000 0.996 0.944 0.681 0.000 F7f72 1.000 1.000 0.953 0.901 0.804 F8f63 0.411 0.300 0.276 0.198 0.027 F9f73 1.000 0.686 0.409 0.279 0.194 F10 f82 1.000 1.000 1.000 1.000 0.967 F11 f92 1.000 0.667 0.667 0.667 0.667 F12 f10 2 0.838 0.743 0.730 0.705 0.593 F13 f11 2 0.997 0.667 0.667 0.667 0.663 F14 f11 3 1.000 0.667 0.667 0.667 0.627 F15 f12 3 1.000 0.723 0.675 0.640 0.378 F16 f11 5 1.000 0.673 0.663 0.663 0.653 F17 f12 5 1.000 0.483 0.453 0.438 0.325 F18 f11 10 0.997 0.470 0.460 0.460 0.327 F19 f12 10 0.273 0.133 0.133 0.127 0.093 F20 f12 20 0.000 0.000 0.000 0.000 0.000 50 ACCEPTED MANUSCRIPT ACCEPTED MANUSCRIPT Table 20: Results with NEA2 Pb Fun Dim Accuracy level 1E−1 1E−2 1E−3 1E−4 1E−5 F1f11 1.000 1.000 1.000 1.000 1.000 F2f21 1.000 1.000 1.000 1.000 1.000 F3f31 1.000 1.000 1.000 1.000 1.000 F4f42 1.000 1.000 1.000 1.000 0.990 F5f52 1.000 1.000 1.000 1.000 1.000 F6f62 0.963 0.963 0.958 0.950 0.000 F7f72 0.946 0.925 0.918 0.914 0.911 F8f63 0.241 0.240 0.240 0.240 0.239 F9f73 0.622 0.595 0.584 0.581 0.579 F10 f82 1.000 1.000 1.000 0.988 0.980 F11 f92 0.980 0.967 0.967 0.960 0.960 F12 f10 2 0.853 0.850 0.843 0.840 0.833 F13 f11 2 0.977 0.970 0.960 0.957 0.947 F14 f11 3 0.830 0.817 0.810 0.807 0.800 F15 f12 3 0.743 0.723 0.720 0.718 0.713 F16 f11 5 0.673 0.673 0.673 0.673 0.673 F17 f12 5 0.695 0.695 0.695 0.695 0.695 F18 f11 10 0.667 0.667 0.667 0.667 0.663 F19 f12 10 0.667 0.667 0.667 0.667 0.667 F20 f12 20 0.363 0.360 0.360 0.360 0.350 51