Adam Gaier

Alexander Asteroth

Jean-Baptiste Mouret

Aerodynamic Design Exploration through Surrogate-Assisted Illumination

A new method for design space exploration and optimization, Surrogate-Assisted Illumination (SAIL), is presented. Inspired by robotics techniques designed to produce diverse repertoires of behaviors for use in damage recovery, SAIL produces diverse designs that vary according to features specified by the designer. By producing high-performing designs with varied combinations of user-defined features a map of the design space is created. This map illuminates the relationship between the chosen features and performance, and can aid designers in identifying promising design concepts. SAIL is designed for use with computationally expensive design problems, such as fluid or structural dynamics, and integrates approximative models and intelligent sampling of the objective function to minimize the number of function evaluations required. On a 2D airfoil optimization problem SAIL is shown to produce hundreds of diverse designs which perform competitively with those found by state-of-the-art black box optimization. Its capabilities are further illustrated in a more expensive 3D aerodynamic optimization task. Nomenclature m() Mean function k() Covariance function GP () Gaussian Process model K Kernel matrix k Covariance vector D Observation set (input/output pairs) µ Mean σ 2 Variance C D Coefficient of Drag C L Coefficient of Lift X up Position of highest point on along upper surface of airfoil Z up Height of highest point on upper surface of airfoil

I. Introduction

Design optimization techniques are often thought of by their creators as the final step in the design process. Algorithms are developed to squeeze every last bit of performance out of a solution, and emphasis is placed on the refinement of designs to their most optimal form. [START_REF] Hornby | Automated Antenna Design with Evolutionary Algorithms[END_REF] If, however, the goal is to support designers then this focus on optimality may very well be misplaced. In an interview study conducted by Autodesk 2 to better understand how professional designers, engineers, and architects use optimization tools they found that optimization was most commonly used not at the end of the design process, but at the beginning. Rather than using optimization algorithms to solve design problems, they were instead used to explore the range of possible solutions. By first generating a range of candidate solutions, designers familiarize themselves with the design space to better understand design concepts and trade-offs.

When faced with a high dimensional design space with many interacting parameters it can be difficult to form an intuition about which combinations of features will lead to high performance. In these cases optimization tools can be employed to produce a batch of design alternatives which best satisfy the objectives and constraints of the problem. These design alternatives act as concrete way points in high dimensional parameter space, helping designers understand various approaches to satisfy their objectives. These designs illustrate the trade offs in the design space, as well as the assumptions and consequences inherent to the problem definition and constraints. As constraints and objectives are reconsidered and adjusted, their effect is made clear by the designs which result from the optimization process. It is within the design regions found by this collaborative human-machine exploration that the designer works, refining designs to match precise requirements along with intangibles such as aesthetics.

The tool most commonly cited in these interviews to produce the variety of designs needed for design exploration was multi-objective optimization. [START_REF] Bradner | Parameters tell the design story: Ideation and abstraction in design optimization[END_REF] If the objectives being considered are in opposition to each other the result of the optimization process is a Pareto front of solutions, each representing a different tradeoff. [START_REF] Deb | Unveiling innovative design principles by means of multiple conflicting objectives[END_REF] It may, however, be the case that during this explorative process the designer's interest is not only in the maximization of performance, but in the effect of different design features on performance, along with the interaction of the features themselves. A better understanding of the problem space can then be gained through experimenting with the effect of different features, and this process can lead to more inspired designs than maximization of objectives alone.

To assist designers in this generative exploration process a new class of algorithms known as illumination algorithms, or quality-diversity (QD) algorithms, 4 could be applied. These algorithms use evolutionary approaches to produce a variety of high-performing solutions, rather than converging on a single optimum. One such algorithm, Multi-dimensional Archive of Phenotypic Elites (MAP-Elites) [START_REF] Cully | Robots that can adapt like animals[END_REF][START_REF] Mouret | Illuminating search spaces by mapping elites[END_REF] explicitly explores the relationship between user-defined features and performance. Designers select a few features deemed interesting or important, such as weight or structural strength, and MAP-Elites produces high-performing solutions which span the possible variations of those features. Known as illumination, this process reveals the performance potential of the different regions of the feature space.

While MAP-Elites is effective at finding a variety of high-performing solutions, this feat is only accomplished at great computational cost. In its best known application MAP-Elites was used to build a repertoire of thousands of different hexapod robot controllers, which then allowed the walking robot to quickly recover from damage to its limbs and motors [START_REF] Cully | Robots that can adapt like animals[END_REF] . While this damage recovery could be achieved online by the robot after only trying a few alternate behaviors, the repertoire itself was produced offline and required twenty million walking simulations. In fluid dynamics applications, where a single simulation can take hours, application of such a computationally intensive algorithm is simply unrealistic.

It is common when faced with computationally expensive optimization problems to make use of approximate models of the objective function, or surrogate models. These models predict the performance of solutions that have not been evaluated, based on the known performance of the solutions which have. [START_REF] Jin | A comprehensive survey of fitness approximation in evolutionary computation[END_REF][START_REF] Forrester | Recent advances in surrogate-based optimization[END_REF][START_REF] Shahriari | Taking the human out of the loop: A review of bayesian optimization[END_REF] Though it is unrealistic to model the entire design space, it is possible to model the high performing regions. Modeling is done by the careful sampling of solutions, coordinated by the optimization algorithm itself. By optimizing an acquisition function which balances exploring new areas of the design space and refining the model in regions already known to perform well, accuracy in high performing regions can be improved and optima found. These computationally efficient models can be used in place of the objective function during optimization, greatly accelerating the process. By incorporating these surrogate-assistance techniques into the evaluation-heavy illumination process, it is potentially possible to use MAP-Elites even in computationally expensive design problems.

In this paper we introduce the Surrogate-Assisted Illumination (SAIL) algorithm for design space exploration. We integrate surrogate-assistance techniques with the MAP-Elites algorithm, allowing it to be used in computationally expensive design problems while maintaining its original capabilities. Before demonstrating its use in a 3D aerodynamics problem we analyze the performance of the algorithm in a simpler 2D airfoil design case, showing that SAIL is:

• Divergent -Produces a diversity of solutions which vary across a user-defined continuum;

• Accurate -Predicts behavior of the objective function in high-performing regions;

• High-Performing -Produces near optimal solutions;

• Efficient -Performs under computational constraints.

Figure 1 describes the operation of the SAIL algorithm: 1) an initial set of observations is produced by sampling the parameter space and evaluating the resulting designs, 2) a surrogate model is produced from these parameter/performance pairs, 3) the MAP-Elites algorithm is used to maximize the acquisition function in every region of the feature space, producing an acquisition map of solutions that are most likely to improve the model in high performing regions, 4) a sample of solutions is drawn from the acquisition map, evaluated, and added to the set of observations. Steps 2-4 are repeated until the the computational budget is reached. The performance predictions of the model can then be used to produce a prediction map of the estimated optimal designs in each region of the feature space, illuminating the performance potential of each region of the feature space and the relationship between features.

1) Sample design space 2) Construct model 3) Maximize acquisition function 4) Sample acquisition map 5) Maximize performance estimation

II. Related Work

Quality Diversity and MAP-Elites

Quality diversity (QD) algorithms 4 use evolutionary methods to produce a set of diverse, high quality solutions within a single run. Rather than seeking a single global optimum, QD algorithms discover as many types of solutions to a problem as possible, and produce a best possible example of each type. For this reason they are also referred to as illumination algorithms, as they illuminate the performance potential of different regions of the solution space.

One such illumination algorithm, the MAP-Elites algorithm [START_REF] Cully | Robots that can adapt like animals[END_REF][START_REF] Mouret | Illuminating search spaces by mapping elites[END_REF] is designed to produce high-performing solutions across a continuum of n user-defined feature dimensions. MAP-Elites first divides the feature space into a grid, or map, of n-dimensional bins. This map holds the set of solutions, with each bin holding a single solution. When the map is visualized, with each bin colored according to the performance of the solution it contains, it provides an intuitive overview of the performance potential of each region of the feature space.

To initialize MAP-Elites a set of random solutions are first evaluated and assigned to bins. The bin to which a solution is assigned is based on its features. If, for example, the feature space is 2D with one dimension for weight and another for cost, a low cost and low weight solution would be placed in the low cost, low weight bin location of the map. If the bin is empty, the solution is placed inside. If the bin is already occupied the two solutions are compared, and the one which performs better is placed in the bin while the other is discarded. Through this competitive process each bin contains the best solution found so far for each combination of features. These solutions are known as elites.

To produce new solutions, elites are chosen randomly, their parameters perturbed by adding Gaussian noise, and then evaluated and assigned a bin based on their features. These new solutions have two ways of being added to the map: discovering an unoccupied bin, or out-competing an existing solution for its bin. Repetition of this process produces an increasingly explored feature space and an increasingly optimal collection of solutions, illuminating the performance potential of the entire feature space. The MAP-Elites algorithm is summarized below in Algorithm 1. X ← X initial , P ← objective_f unction(X initial )

Algorithm 1 MAP-Elites

5:

for iter = 1 → I do 6:

x ← random_selection(X )

7:

x ← random_variation(x)

8: b ← f eature_descriptor(x ) 9:
p ← objective_f unction(x ) end for 14:

return (X , P) Return illuminated map 15: end function MAP-Elites has been shown to be effective in exploration and optimization in a variety of domains including the design of walking soft robot morphologies, [START_REF] Mouret | Illuminating search spaces by mapping elites[END_REF] the generation of images that fool deep neural networks, [START_REF] Nguyen | Deep neural networks are easily fooled: High confidence predictions for unrecognizable images[END_REF] and the evolution of robot controllers capable of adapting to damage. [START_REF] Cully | Robots that can adapt like animals[END_REF]

Surrogate-Assisted Optimization

Evolutionary optimization approaches typically require a large number of evaluations before acceptable solutions are found. In many applications these performance calculations are far from trivial, and the computational cost of repeated evaluations becomes prohibitively expensive. In these cases approximate models of the objective function, or surrogate models, are used in their place. Surrogate-assisted optimization has been a particularly useful approach in the computationally demanding context of fluid dynamics. [START_REF] Giannakoglou | Aerodynamic shape design using evolutionary algorithms and new gradient-assisted metamodels[END_REF][START_REF] Dumas | CFD-based optimization for automotive aerodynamics[END_REF][START_REF] Obayashi | Visualization and data mining of Pareto solutions using self-organizing map[END_REF][START_REF] Jeong | Data Mining for Aerodynamic Design Space[END_REF] Due to the sheer number of evaluations required by MAP-Elites, even when evaluations are inexpensive surrogate-assistance has the potential to accelerate the illumination process dramatically.

Surrogate-assisted optimization often takes place within the wider framework of Bayesian optimization (BO). [START_REF] Cully | Robots that can adapt like animals[END_REF][START_REF] Shahriari | Taking the human out of the loop: A review of bayesian optimization[END_REF][START_REF] Brochu | A tutorial on Bayesian optimization of expensive cost functions, with application to active user modeling and hierarchical reinforcement learning[END_REF] BO approaches the problem of optimization not just as one of finding the most optimal solution, but of modeling the underlying objective function in high-performing regions. To estimate the objective function probabilistic models are used, giving each sample a predicted objective value (mean) and a confidence in that prediction (variance). New samples are chosen where the model predicts a high objective value and where prediction uncertainty is high. How the emphasis of predicted performance and prediction uncertainty is balanced is determined by an acquisition function. In the active learning process of BO, the sample which maximizes the acquisition function is chosen as the next observation, evaluated, the result included in the set of observations, the model updated, and the process repeated.

A variety of data-driven machine learning techniques such as polynomial regression, support vector machines, and artificial neural networks can be used to construct surrogate models. [START_REF] Jin | A comprehensive survey of fitness approximation in evolutionary computation[END_REF][START_REF] Forrester | Recent advances in surrogate-based optimization[END_REF] However, for model-based methods that require probabilistic results, such as Kriging [START_REF] Jeong | Efficient optimization design method using kriging model[END_REF] and BO, [START_REF] Brochu | A tutorial on Bayesian optimization of expensive cost functions, with application to active user modeling and hierarchical reinforcement learning[END_REF] Gaussian process (GP) models 17 are typically used.

Gaussian Process Models

In the active learning context of surrogate-assisted optimization a measure of model uncertainty is essential, as this allows for the balancing of exploring unknown regions of the search space and exploiting regions known to contain high-performing solutions. GP 17 models are often used as for surrogate-assisted optimization as they are effective even when only a small number of samples are available and, crucially, their predictions include a measure of certainty.

Gaussian process models are a generalization of the Gaussian distribution: where a Gaussian distribution describes random variables, defined by mean and variance, a Gaussian process describes a random distribution of functions, defined by a mean function m, and covariance function k.

f (x) ∼ GP (m(x), k(x, x )) (1)
In much the same way as a neural network model can be thought of as a function that returns a scalar value given an arbitrary input vector x, a GP model can be thought of as a function that, for a given x, returns the mean and variance of a normal distribution, with the variance indicating the certainty of the prediction.

Gaussian process models calculate predictions of unknown samples based on their distance to known samples in input space, a relationship defined by a covariance function. A common choice is the squared exponential function: the closer the points are in input space the more closely correlated they are in the output space:

k(x i , x j ) = exp - 1 2 x i -x j 2 (2)
Given an observation set D = (x 1:t , f 1:t ) where f 1:t = f (x 1:t ), a matrix of covariances is built. In the simple noise-free case this results in the kernel matrix:

K =     k(x 1 , x 1 ) • • • k(x 1 , x t ) . . . . . . . . . k(x t , x 1 ) • • • k(x t , x t )     (3)
When considering a new point (x t+1 ) the value (f t+1 = f (x t+1 )) can be derived using the normal distribution built from this kernel matrix a :

f 1:t f t+1 ∼ N 0, K k k T k(x t+1 , x t+1 ) (4)
where k = [k(x t+1 , x 1 ), k(x t+1 , x 2 ), . . . , k(x t+1 , x t )] T allowing us to compute the GP as:

P (f t+1 |D 1:t , x t+1 ) = N µ t (x t+1 ), σ 2 t (x t+1 ) (5)
a for simplicity we assume the zero mean function: m(x) = 0

where:

µ t (x t+1 ) = k T K -1 f 1:t (6) σ 2 t (x t+1 ) = k(x t+1 , x t+1 ) -k T K -1 k (7)
gives us the predicted mean and variance for a normal distribution at the new point x t+1 . If we were then to evaluate the objective function at this point, we would add it to our set of observations D, reducing the variance at x t+1 and at other points near to x t+1 . In this pure generalized form, our GP model weighs variations in every dimension equally, applying the same squared exponential relationship regardless of input dimension. For higher dimensional problems each dimension's effect on the output can also be weighted via a technique known as automatic relevance detection (ARD). The hyperparameters which weigh each dimension are set by maximizing the likelihood of the model given the data. [START_REF] Rasmussen | Gaussian Process for Machine Learning[END_REF] This increases model accuracy, and the weighting provides an human-readable estimation of the relative importance of each input dimension.

Visualization of Surrogate Models

Effective surrogate models are precise enough to be used by optimization algorithms, but as black-box models of high dimensional problem spaces, the knowledge they contain is inherently difficult to grasp. Despite the wealth of information encoded within these models, few attempts in the literature have been made to make this knowledge accessible.

Obayashi [START_REF] Obayashi | Visualization and data mining of Pareto solutions using self-organizing map[END_REF] and Jeong [START_REF] Jeong | Data Mining for Aerodynamic Design Space[END_REF][START_REF] Jeong | Efficient optimization design method using kriging model[END_REF] use self-organizing maps (SOMs) to project the results of single surrogateassisted optimization runs from a high dimensional parameter space to a 2D space based on similarity of the samples in the observation set in parameter space. This allowed different design regions to be identified, as well as the design trade-offs they represent. In work by Jouhaud 18 an airfoil design optimization problem was solved using a surrogate model, and the resulting model visualized. As only two parameters were optimized and used by the model the search space could be easily sampled and presented as a 2D plane with the fitness, drag, and lift of different parameter regions shown by isolines.

Neither of these approaches can be used to intuitively visualize a high-dimensional surrogate model. In the approach of Obayashi and Jeong only precisely evaluated solutions are visualized, making the SOM not a description of the design space as interpreted by the current surrogate model, it is instead a history of the surrogate model's construction. While Jouhaud's approach provides the encoded knowledge of the surrogate model in an easily understood form, it is only a viable technique because of the low dimensionality of the problem, and cannot be scaled to higher dimensional design problems. By making use of a surrogate model to perform illumination of the feature space, SAIL allows intuitive visualization of the surrogate model's knowledge through the lens of chosen low-dimensional features.

III. Surrogate-Assisted Illumination

SAIL operates by modeling the underlying objective function in different high-performing regions of the search space. Knowledge of the objective function can then be applied to the production of high performing designs in different regions of the feature space. These high performing designs can in turn be examined as a whole in order to improve our understanding of the relationships between features, performance, and parameters. Deliberate sampling of the objective function to better understand its behavior and locate its optima is also the goal of Bayesian optimization (BO), [START_REF] Shahriari | Taking the human out of the loop: A review of bayesian optimization[END_REF][START_REF] Brochu | A tutorial on Bayesian optimization of expensive cost functions, with application to active user modeling and hierarchical reinforcement learning[END_REF] and we borrow from its methods and vocabulary.

BO requires a probabilistic model of the objective function, such as a GP, and an acquisition function, which describes the utility of sampling a given point. The model of the objective function is produced from previously observed samples, the point with the greatest utility is found, evaluated, and added to the set of observations. The model is then reconstructed with these additional samples, and the process repeated.

It is infeasible to model the entirety of the search space when evaluations are computationally expensive, so sampling is directed towards the interesting regions of the search space. While in BO this interesting region is the region around the highest performing point, in SAIL we would like to find optima in all regions of the feature space, and so must model the region around a high performing slice.

How we define utility in the acquisition function determines how the search space is sampled. Balance should be maintained between exploration, sampling points with high uncertainty, and exploitation, sampling points which are likely to have high performance. In SAIL, the upper confidence bound (UCB) [START_REF] Srinivas | Gaussian Process Optimization in the Bandit Setting: No Regret and Experimental Design[END_REF] is used as an acquisition function. UCB judges potential observations optimistically, favoring uncertainty under the assumption that higher uncertainty hides a potentially higher reward. A high mean (µ(x)) and large uncertainty (σ(x)) are both favored, with relative emphasis tuned by the parameter κ.

U CB(x) = µ(x) + κσ(x) (8)
UCB performs competitively with more complex acquisition functions such as Expected Improvement (EI) and Probability of Improvement (PI). [START_REF] Brochu | A tutorial on Bayesian optimization of expensive cost functions, with application to active user modeling and hierarchical reinforcement learning[END_REF][START_REF] Calandra | An experimental comparison of Bayesian optimization for bipedal locomotion[END_REF] These acquisition functions rely on comparisons to the current optimum, while UCB is based only on the underlying model. This is more appropriate for SAIL, as solutions in lower performing regions of the feature space would all have vanishingly small probabilities of improving on the optima found in the highest performing regions. A more local approach, such as calculating EI only within bins, is also fraught with difficulties. In many cases there will be more bins than observations, and many bins will not contain previous observations for comparison, leaving local performance measures incalculable.

To explore a diversity of designs SAIL models the objective function around high-performing solutions that cover the entire feature space. In order to predict performance in this entire slice of the search space, we first produce designs with every combination of features. This is done by dividing the feature space into bins, and using MAP-Elites (see Section II: Quality Diversity and MAP-Elites) to produce designs which maximize the acquisition function in each bin. The resulting map, which contains designs with maximal predicted performance and uncertainty, we call the acquisition map.

It is from this acquisition map that we select new samples to be evaluated. In order to draw new observations from the feature space in a uniform fashion, we employ a Sobol sequence. [START_REF] Niederreiter | Low-discrepancy and low-dispersion sequences[END_REF] This space-filling sequence allows us to select bins in an increasingly dense, but evenly spaced manner. Once evaluated the performance of the samples can be added to our set of observations and the GP model reconstructed. We can then produce a new acquisition map with the updated GP model, and repeat the process until our evaluation budget is exhausted. GP ← Gaussian_process_model(D)

Algorithm 2 Surrogate-Assisted Illumination

8:
acquisition() ← (U CB(), GP())

9:

(X acq , P acq ) = MAP-Elites(acquisition(), X ) Select solutions from acquisition map for PE 10:

x ← X acq (Sobol iter )

11:

X ← X ∪ x, P ← P ∪ P E(x) 12: end for 13: 14: 2) Produce Prediction Map 15: D ← (X , P)

Observation Set: Parameters, Performance 16: GP ← Gaussian_process_model(D) 17: prediction() ← mean(GP(x)) 18: (X pred , P pred ) = MAP-Elites(prediction(), X )

The SAIL algorithm is more precisely defined in Algorithm 2. An initial set of designs is first produced. For maximum coverage of the parameter space we again use a Sobol sequence, [START_REF] Niederreiter | Low-discrepancy and low-dispersion sequences[END_REF] this time to ensure even coverage of the parameter space. These designs and their performance form a set of observations D, which is used to construct a GP model. An empty acquisition map is then created and filled with the designs from D, along with their utility as judged by the acquisition function. These designs are taken as the starting population for MAP-Elites (Algorithm 1) which then illuminates the map as described in Section II: an elite is selected and its parameters perturbed to produce a new solution, it is assigned a bin based on its features, and it finally competes for the bin if it is not occupied. This illumination process repeats for a number of iterations, and results in an acquisition map of elite solutions who each maximize the acquisition function in their bin.

From the acquisition map we select the next samples for evaluation. To ensure even coverage of the feature space, we again employ a Sobol sequence to direct the sampling, this time producing coordinates in feature space rather than parameter values. These coordinates indicate the bin to be sampled, and the individual stored is precisely evaluated on the true objective function. Once evaluated these new designs and fitness pairs are added to our observation set D and the process can be repeated.

The mean prediction of the resulting GP model can then used as the objective function of MAP-Elites, and a prediction map produced. This map is an estimate of the relationship between features and performance, including a predicted optimal design for each bin. As only the surrogate model of the objective is used for evaluation, this prediction map can be produced with minimal computation.

Not only does this prediction map contain a variety of high-performing designs it also presents the effects of features on performance in an intuitive visual format, allowing optimal or interesting regions to be easily identified. Parameter values of the found optimal designs can also be viewed through this feature space lens, allowing different design regions to be distinguished and providing insight into the solution representation. If, for example, the optimal parameter values are all at the extreme of the allowed range, its bounds are likely too restrictive. If optimal parameter values across the feature space appear noisy and randomly distributed, it can be deduced that the parameter value has little effect on performance, either because it is genuinely unimportant or the parameter range is so restrictive that variation within that range has little effect.

IV. Experimental Setup

Objectives and Constraints

To test the ability of SAIL to create a variety of high performing designs we evaluate its performance on a 2D airfoil optimization problem. Our objective is the reduction of the C D of an airfoil, while matching the area and lift of the high-performing RAE2822 airfoil. Foils are evaluated at an angle of attack of 2.7 • , at Mach 0.5 and Reynolds number of 10 6 . To ensure that area and lift constraints are closely followed we include quadratically increasing fitness penalties for violations. As the C D values are quite small we take the negative logarithm of C D as our drag, with higher values being more desirable. The evaluation criteria are defined as:

fitness(x) = drag(x) × penalty lift (x) × penalty area (x) (9)
where drag(x) = -log(C D (x))

penalty lift (x) =    C L (x) lift RAE2822 2 , if C L (x) < lif t RAE2822 1, otherwise (10)
penalty area (x) = 1 - |area -area RAE2822 | area RAE2822 7 (11)
During optimization both C D and C L values are approximated with Gaussian Process models. The UCB of the drag prediction is taken as the drag component:

drag(x) = µ drag (x) + κσ drag (x) (12)
We treat the lift prediction as a classification problem. Based on the mean and variance of the predicted lift, we penalize foils based on the probability that their lift is less than that of the base RAE2822 foil. The GP prediction based lift penalty is defined as:

penalty lift (x) = 1 -P (C L (x) < lif t RAE2822 ) (13)

Representation

Airfoils are encoded using the airfoil-specific PARSEC parameterization. [START_REF] Sobieczky | Parametric airfoils and wings[END_REF] Polynomial expressions are used to encode design features, such as the height of the airfoil and the angle of the trailing edge. The PARSEC representation allows a wide variety of designs to be expressed with a small number of design parameters. The sharpness and end point of the trailing edge is fixed to that of the base RAE2822 foil, and the remaining 10 parameters are shown below in Figure 2. Illumination algorithms allow us to explore the design space along a continuum in feature space. In this example we choose two of the PARSEC parameters: the height of the highest point on the top side of the foil (Z up ), and the location along the length of the wing of this highest point (X up ). These two parameters were chosen because of the large effect on the drag of the foil they were found to have in early tests. We divide the possible range of Z up and X up into 25 partitions, resulting in a 25×25 grid, for a total of 625 bins.

Zxx lo rLE up rLE lo Z up Z lo α TE β TE X lo

X up

Zxx up

In practice we would typically define feature dimensions that do not correspond to parameter values. This allows the exploration of higher level features that are the result of parameter combinations, enabling the projection of a high dimensional representation onto a lower feature dimensional space of the designer's choosing (see discussion in Section VI: Representation). These higher level features should be chosen based on the designer's experience, or to explore a particular relationship. Here we select parameter values as features so that performance to enable comparison to other, more traditional, optimization algorithms.

Baseline and Hyperparameters

As one of the primary objectives of the SAIL algorithm is to achieve results in computationally expensive domains, the unit of comparison used is the number of precise evaluations (PE) required to achieve a certain result. Model construction and other algorithmic costs are ignored and assumed to be dominated in truly expensive optimization problems. The optimality of the designs in the prediction maps produced by SAIL are compared to those produced by 1) the original MAP-Elites algorithm, 2) the black-box optimization algorithm Covariance Matrix Adaptation Evolution Strategy (CMA-ES), [START_REF] Hansen | Completely derandomized self-adaptation in evolution strategies[END_REF] and 3) a surrogate-assisted variant of CMA-ES (SA-CMA-ES).

SAIL is given a total computational budget of 1000PE. 50PE is used to evaluate a set of designs produced through parameter sampling which form the basis of the initial GP model. At each iteration (Algorithm 2 lines 4-12) 10 additional individuals are evaluated and incorporated into the GP model until the computation budget is exhausted. This result was compared with the MAP-Elites algorithm given a computational budget of 10 5 PE.

To better understand the true performance of our algorithm in finding optimal solutions, we compare it to the black-box optimization algorithm CMA-ES, [START_REF] Hansen | Completely derandomized self-adaptation in evolution strategies[END_REF] which is designed to produce a single solution. As the dimensions of variation correspond to parameter values, by restricting the allowed parameter ranges of the representation we can confine the search within a single cell. Each bin is treated as an independent search problem which CMA-ES is given a budget of 1000PE to solve. In addition we employ a surrogate-assisted variant of CMA-ES, SA-CMA-ES, to solve each of these subproblems. Within a single bin 25 initial solutions are sampled from the parameter space and used to produce an initial GP model. CMA-ES is then used to maximize the UCB based acquisition function described above. The single found optimum is evaluated and incorporated into the GP model. This process is repeated until a budget of 100PE is exhausted.

Each approach was replicated 20 times with different random seeds. b Unless otherwise mentioned, all values are medians across these replicates. As valid initial designs whose highest point was found at the leading edge of the wing (high Z up and low X up ) could not be found due to geometric constraints inherent to the PARSEC representation, [START_REF] Padulo | Airfoil Design Under Uncertainty with Robust Geometric Parameterization[END_REF] only the remaining 577 bins were considered in comparisons. Where possible, standard implementations of algorithms were used, including CMA-ES, [START_REF] Hansen | CMAES Version 3.61.beta[END_REF] Gaussian Processes, [START_REF] Rasmussen | Gaussian Process Regression and Classification Toolbox[END_REF] and airfoil simulation. [START_REF] Drela | XFOIL airfoil simulator[END_REF]

Accuracy

At the conclusion of the SAIL algorithm a prediction map is produced by illuminating the feature space without rewarding uncertainty. This gives us our model's 'best guess' of the optimal designs in the feature space. In this two-dimensional airfoil case it was possible to precisely evaluate every design in the prediction map and so judge its accuracy in different regions of feature space. The median prediction, true values, and their absolute percentage error are shown in Figure 4 below.

Over the majority of the feature space, especially for C D , the predictions made by the model are reliably accurate. 90% of drag (log(C D )) predictions and more than 80% of lift (C L ) predictions are within 5% of their true value. Drag errors are clustered in one region of the feature space, where the flow simulator was less likely to converge and produce valid results. In these experiments when the simulator did not converge the result was simply discarded and the next sample taken, leaving this region with fewer samples on which to base a prediction.

To allow for effective optimization our models are designed to accurately predict performance in optimal regions of the design space across a slice of the feature space. This is done by directing the sampling process toward high performing solutions in different regions of the feature space. To measure the accuracy of our models in this high fitness slice we test their ability to accurately predict the performance of the optimal solutions found by CMA-ES. We compare the models built through selection guided by SAIL to models built by evenly sampling the parameter space with a Sobol sequence. To further tease apart the factors which contribute to SAIL's modeling of the feature space we also examine the effect of the choice of acquisition function. As discussed in Section III, we guide selection through the upper confidence bound (UCB), a combination of the mean and variance of the prediction. This should guide the search process to promising unexplored regions. The MAP-Elites illumination process is also designed to find diverse high quality solutions. To differentiate the contribution of each to directing the search towards the high fitness slice we test SAIL using an acquisition function of only the mean, and only the variance.

The accuracy of each resulting model's drag prediction on the optimal design set is measured at various stages of the sampling process. Figure 5 shows that using the SAIL algorithm to select samples from across the feature space improves accuracy over evenly sampling the parameter space, regardless of whether uncertain or high-performing solutions are favored. When both the uncertainty and performance are considered, as it is with the UCB, SAIL produces models which are an order of magnitude more accurate on the high fitness slice than uniform sampling. Mean squared error (log scale) of drag prediction on optimal designs. Models constructed using designs sampled from parameter space using a Sobol sequence or selected from acquisition maps produced with the mean, variance, or the UCB of the prediction.

Accuracy on High Fitness Slice

Number of Samples

Optimality and Efficiency

It is not our goal to compete with algorithms designed to find one very high-performing solution. If we are to accurately portray the performance of feature space, however, it is required that the solutions which we produce represent the best designs in the region. We again rely on the designs produced through repeated applications of CMA-ES as our benchmark of optimality. The median designs found by SAIL in each bin after 1000PE were compared to the best found by CMA-ES in each bin after 20 runs for 1000PE in each bin (≈ 11.5 million PE in total). In Figure 6 we see the median predicted fitness in each bin, the true fitness of the designs found after evaluating them, and the 'optimal' fitness in each bin found by CMA-ES. Though the solutions found by CMA-ES have higher performance, as we would expect from a state-ofthe-art black box optimization algorithm using four orders of magnitude more evaluations, SAIL provides a near-optimal approximation of the feature space. The fitness potential of each bin is well illuminated: found designs in nearly half the bins are within 5% of the optimum, and the general relationship between features is accurately portrayed.

We are unaware of such feature space exploration algorithms beyond MAP-Elites, and so judge the data efficiency of our algorithm in comparison to CMA-ES. CMA-ES was not designed for use across a multitude of subproblems, and so the number of PEs required to arrive at an optimized feature map is highly dependent on the number of bins that are defined in the map. For a more informative comparison we therefore compare SAIL to the performance of CMA-ES within a single bin as well as its performance in optimizing the map as a whole. We include a surrogate-assisted variant of CMA-ES (SA-CMA-ES) (see Baseline and Hyperparameters, previous section) to control for the ease of modeling the problem, and also compare its performance in a single bin and across the entire map.

Figure 7 shows these comparisons of performance per evaluation. As optimization progress may vary according to the bin, we take the single bin performance as the median performance over all bins. Performance itself is measured as the percent of the optimum found by CMA-ES in each bin. As this optimum is the highest performance over all runs, the median CMA-ES run does not typically reach this value. Performance of the individuals produced to construct the initial surrogate models is set to 0%, with the first valid performance indicators at 25PE/bin for SA-CMAES and 50PE for SAIL.

Given the budget required by CMA-ES to find a near optimal solution in a single bin, SAIL finds designs of comparable performance in every bin. Comparisons between CMA-ES and MAP-Elites and their surrogateassisted variants SA-CMA-ES and SAIL reveal that the gains from surrogate-assisted optimization are even greater for MAP-Elites than the traditional optimizer. Even when MAP-Elites is given an evaluation budget two orders of magnitude greater than SAIL it cannot produce solutions of similar performance.

98.9% (Bin

VI. Case Study: Velomobile Design

Having established the capabilities of the SAIL algorithm in a 2D case, we now present design exploration results in the more computationally demanding context of 3D aerodynamics. We apply the SAIL algorithm to the problem of creating aerodynamic shells for fully faired recumbent bicycles. Such streamlined vehicles hold human powered top speed (144.17 km/h c ) and distance (680km in 12 hours, 1219km in 24 hours d ) records due to their highly tuned aerodynamics. Due to constraints such as rider movement and comfort, three-wheeled designs built for distance races, known as velomobiles, often have non-intuitive shapes (see Figure 8). These high-performing, but odd, designs suggest that the domain may be rich in interesting design concepts. By producing a set of varied designs, we illustrate the advantages of SAIL over approaches focused purely on optimization or exploration, and within the limited computational budget enforced by the use of 3D CFD simulation.

Representation

In this work we focus on a simplified base body of a velomobile only, without more complex structures such as wheel housings or a cockpit for the rider's head. We describe this three dimensional shape with a series of airfoil-like curves defined with the PARSEC 22 representation (see Section IV: Representation). Viewed from the top, the vehicle is defined as a symmetrical airfoil (Top curve). The side profile is defined by three curves: one fixed curve for the bottom and two parameterized curves for the top, one in the center of the vehicle (Mid curve) and one which forms a 'ridge' for the knees (Ridge curve). A final curve connects this ridge to a flat bottom to form the view from the front (Rib curve). This Rib is defined along a unit vector, and is scaled to maintain the curves defined by the Top and Ridge curves in 32 sections. From the top view ridges follow the same curve as the Top curve of the body. All curves are connected to each other by splines. An example of these curves and the 3D shape they produce is illustrated below in Figure 9. Each PARSEC airfoil can be described by 10 or more parameters, but as we are only producing one curve rather than an entire foil, only the parameters which describe one side are necessary. After initial tests with SAIL we examined the automatic relevance detection coefficients in our Gaussian Process models (see Section II: Gaussian Process Models). From the values of these coefficients we could identify parameters which had negligible influence on performance across the entire features space. By fixing these less important parameters we were able to reduce the representation to the 16 values in Table 1. Reducing the number of parameters allows both optimization and modeling of the design space to be performed more effectively.

Objective and Dimensions of Variation

The produced velomobile shells are judged on the drag force they produce when traveling at 20 m/s (72 km/h). Flow simulation is performed with the OpenFOAM e CFD toolbox. We explore two dimensions of variation: volume and curvature. Volume is computed as the convex hull of the design, curvature is the average curvature of all ribs. Volume is necessary for both the rider and the equipment inside of the shell. It is obvious that larger volume shapes will tend to have a larger amount of drag, but it is useful for a designer to be able to select a high-performing design, and if it is too small, select a slightly larger one. Curvature of the design can allow for a smaller frontal area while still maintaining volume, for example to accommodate the rider's knees. It also has the effect of adding structural integrity to the light-weight shell. In areas where the shell is flat it is also weak, enough to buckle and alter shape at high speeds or in high winds, as well as being easily damaged in accidents.

VII. Results

SAIL was run with a budget of 1000 PE, with 200 PE used for the initial sampling, and 10 new samples drawn after illumination. The volume and curvature feature space was divided into 20 sections, producing a 20 × 20 feature grid. To analyze the true performance of the algorithm we evaluated all 400 designs in the prediction map at the end of the run. As SAIL relies on stochastic processes for optimization this was repeated 15 times: unless otherwise mentioned all values refer to the median performance of these replicates. The prediction map produced along with the true performance of the designs in the prediction map are shown in Figure 10. While the model is optimistic in its predictions, the error of the prediction is quite low, with the majority of designs performing only .15N less than predicted, and clearly captures the unsurprising relationship between increased volume and greater drag force. In Figure 11 a selection of designs produced in a single run from different regions of the map are shown. We can see the variety of designs in the two specified axes of variation. A clear continuum can be observed, with each design similar to its neighbors, and those at opposite ends of the feature space having little in common. This variety allows a designer to quickly scan through a collection of high-performing designs according to characteristics and requirements of their own choosing.

For example, the design at the bottom right, while performing very well has such a low volume that a rider may not fit inside. The designs in the middle row, while large enough for the rider, may require some adjustment if the proper gearing is to fit inside. In the case that adjustments to the shell or creative configuration of the internal mechanisms is not sufficient, a larger high-performing design such as the shell in the middle of the top row can be taken as a starting point.

More Explorative than Pure Optimization

When producing a single design a traditional optimizer will perform better than a design space exploration algorithm, but is seriously limited in its ability to explore according to the interest of the designer. This is despite the fact that an optimizer will not always converge on identical solutions. CMA-ES was given the same 1000PE budget to find optimal solutions for our problem and, depending on the initial starting solutions, converged on different designs. Two such designs are shown in Figure 12 These designs are in many ways dissimilar, one with a single gentle curve on the side, and the other favoring high curvature. The progress of the optimization was tracked and the parameter values of the best found individuals is shown at the bottom of Figure 12. Each run shows convergence on different optima, each in different areas of the parameter space.

While this seems to make the case for exploration using traditional algorithms, if optimization progress is instead tracked in feature space rather than parameter space (Figure 13) an altogether different picture emerges. As volume is closely linked with drag force the search for optimality quickly converges on the smallest designs allowable by the representation. The search is confined to a very small area of the feature space, leaving us with only a narrow range of designs. This may be desirable behavior in an optimization context, but provides little assistance to a designer trying to understand the problem space more broadly.

To a designer the optimal volume may depend on other constraints not easily inserted into the algorithm, such as how other mechanisms and equipment could be arranged within the shell. By producing a selection of designs that vary according to their needs the designer can find promising designs by exploring a menu of options.

More Optimal than Pure Exploration

To produce representative designs across the feature space, sampling of the parameter space alone will typically not be sufficient. In the earlier airfoil experiments the feature space was a subset of the parameter space. In that case evenly sampling the parameter space produced an even sampling of the feature space. In this case, however, the volume and curvature features are the result of a combination of parameters. This causes evenly spaced sampling of the parameter space to result in an uneven sampling of the feature space.

To illustrate this, 1000 samples were produced by sampling the parameter space with a Sobol sequence, recording the resulting designs position in feature space. This process was repeated 10 times and the number of times each bin was discovered recorded. The frequency at which each bin was sampled over each run is shown in the left of Figure 14, revealing a clear bias in parameter space toward low curvature shapes, in addition to the extremes of the feature space being left undiscovered in every attempt.

Even when the feature space is well sampled the designs are not necessarily 'good' examples of designs in this feature region. As volume is tightly coupled with performance their general relationship is clearly visible from parameter sampling, but as examples of high-performing designs they are insufficient. There is no guarantee that the designs found by sampling are any better than average for their bin, as no optimization of any kind has taken place. The median fitnesses in Figure 14 display this volume/drag force relationship so smoothly only because they are median results of repeated sampling. In the case of more limited sampling the results can easily be much noisier. Furthermore as the designs which fill each bin were produced pseudo-randomly without any method of improvement, searching through these solutions for common design characteristics makes little sense.

SAIL simultaneously explores the feature space and optimizes within each bin. This allows SAIL to not only more thoroughly explore the feature space, but produce high-performing examples of designs for each feature combination. The improvement of designs produced by SAIL over those found by parameter sampling is shown at the bottom of Figure 14.

VIII. Conclusion and Discussion

A design space exploration algorithm, Surrogate-Assisted Illumination (SAIL), was introduced. By combining surrogate-assisted optimization techniques with quality-diversity methods SAIL produces a large variety of high-performing solutions with limited computational effort. Produced designs vary according to features described by the user, giving the designer the flexibility to explore design characteristics that interest them, allowing their own domain knowledge and intuition to lead them to interesting design regions. By producing accurate models of the design space and designs which those models predict will perform well, an explanatory map of the relationship between features and performance is produced. By accelerating the illumination process, exploration of the design space can be done rapidly and in an any-time fashion. Once new samples are available the surrogate model can be reconstructed and the estimated performance potential of the entire feature space updated. This 'real-time' view of the search space differentiates the SAIL algorithm from techniques which visualize the precisely evaluated solutions used to produce a model. By visualizing the current best predictions rather than the set of evaluated solutions we see the search space as the model currently understands it; the evaluated solutions are only a history of the optimization process.

The same model can be used to quickly produce prediction maps of different resolution and different features combinations. This allows more detailed exploration of interesting areas, and visualization of the design space through various feature lenses. Data mining techniques designed to extract common design principles from a set of high-performing designs can still be used in conjunction with SAIL. That designs are produced along a continuum defined by the designer provides the added ability to choose elements of variation in high-performing designs to examine more closely.

The capability to rapidly understand the performance potential of the design space through concrete highperforming examples is a potential boon to designers. SAIL not only accelerates the generative design cycle, but allows the effect of user-defined features to be examined, adding new flexibility to cooperative humanmachine design exploration. Generative design tools such as SAIL, which consider more than objectives, can help designers understand what is possible, beyond what is optimal.

Figure 1 . 2 ) 3 )

 123 Figure 1. Surrogate-Assisted Illumination (SAIL) 1) Sample design space to produce initial solutions. 2) Construct model of the objective function based on samples. 3) Maximize the acquisition function, which balances exploitation and exploration, in every region of the feature space, producing an acquisition map. 4) Draw the next samples to test on the objective function from the acquisition map. Repeat steps 2-4. 5) Maximize fitness as predicted by the resulting model, producing a prediction map of high-performing designs.

1 :

 1 function MAP-Elites(objective_f unction(), X initial )

10 :

 10 if P(b ) = ∅ or P(b ) < p then 11: P(b ) ← p , X (b ) ← x

1 :

 1 Initialize with G solutions drawn from Sobol sequence 2: X ← Sobol 1:G , P ← P E(X ) PE = precise evaluation 3: 4: 1) Produce Acquisition Map 5: for iter = 1 → precise_evaluation_budget do 6: D ← (X , P) Observation Set: Parameters, Performance 7:

Figure 2 .

 2 Figure 2. The 10 parameters used to define an airfoil. Dimensions of variation (Xup and Zup) in gold.

21

 21

Figure 4 .

 4 Figure 4. Drag and Lift Predictions and True Value in Each Bin Predicted and true values of drag and lift for designs in each bin after 1000PE. Median values over 20 replicates.

Figure 5 .

 5 Figure 5. Accuracy of Sampling StrategiesMean squared error (log scale) of drag prediction on optimal designs. Models constructed using designs sampled from parameter space using a Sobol sequence or selected from acquisition maps produced with the mean, variance, or the UCB of the prediction.

FitnessFigure 6 .

 6 Figure 6. Performance of Designs Found By SAIL Top: Median predicted and true performance of designs found by SAIL with a budget of 1000PE compared with performance of optimal designs found by 20 runs of CMA-ES per bin (≈11.5 million PE) Bottom : Optimality of SAIL designs per bin.

Figure 8 .

 8 Figure 8. Record-setting velomobiles have non-intuitive shapes

Figure 9 .

 9 Figure 9. A three-dimensional velomobile shape produced from a set of curves

Figure 10 .

 10 Figure 10. Predicted performance of designs found by SAIL, the true performance of the designs found through simulation, and the absolute error in each bin. Median results after 15 replicates each with a budget of 1000 PE.

Figure 11 .

 11 Figure 11. A selection of designs produced in a single run of SAIL. Left: Lower Curvature, Up: Increased Volume Parenthesis after category indicate the bin where each design is located

Figure 12 .

 12 Figure 12. CMA-ES Shapes found in two separate runs using 1000PE with different starting positions. Top: Designs and their component curves Bottom: Parameter values of best performing design found at each CMA-ES iteration in the separate runs, warmer colors indicate later iterations.

Figure 13 .

 13 Figure 13. Feature space exploration in two sample runs of CMA-ES. First iteration in which a design is found in each bin. White bins indicate that no solution was ever ever found with this combination of features.

Figure 14 .

 14 Figure 14. Left: Probability of sampling a bin with 1000 samples drawn from a Sobol sequence. Middle: Median drag force values found by sampling the parameter space. Right: Median drag force values found in a median run with SAIL. Bottom: Improvement per bin of SAIL over parameter sampling with a Sobol sequence. Only bins which were sampled by parameter sampling 3 or more times in the 10 trials are shown.

  Optimization Efficiency in a Single Bin and Over the Entire Design Space Computational efficiency of CMA-ES, SA-CMA-ES, MAP-Elites, and SAIL measured in in precise evaluations. Bin : median progress towards optimum in every bin. Map: performance of CMA-ES and SA-CMA-ES is median bin performance multiplied by number of valid bins. Performance of individuals produced to construct initial models is set to 0%. Bounds indicate one standard deviation over 20 replicates. PEs and performance in log scale.

			Optimization Performance Per Precise Evaluation
	10 1 99% 90% Percentage of Optimum 0%	10 2	) SA-CMA-ES	98.5% (Bin) CMA-ES 97.6% (Map) SAIL 10 3 Precise Evaluations	98.9% (Map) SA-CMA-ES 10 4	91.9% (Map) MAP-Elites 98.5% (Map) CMA-ES 10 5
		CMA-ES	SA-CMA-ES	SAIL	MAP-Elites
	Figure 7.

Table 1 .

 1 The 16 parameters used to determine the shape of a velomobile shell. The tallest point of the ridge is based on the height of the middle point.

	Curve		PARSEC Parameter
	Top Leading Edge Radius	Position of Tallest Point Tallest Point	Curvature
	Mid Leading Edge Radius	Position of Tallest Point Tallest Point	Curvature
	Ridge Leading Edge Radius	Position of Tallest Point Relative Height
	Rib	Leading Edge Radius Angle of Trailing Edge	Position of Tallest Point Tallest Point	Curvature

of 20 American Institute of Aeronautics and Astronautics

b One replicate, including data gathering, with 8 cores of a Intel Xeon 2.6GHz processor required approximately: SA-CMA-ES:32h, CMA-ES:80h, SAIL:12h, MAP-Elites:14h

c http://www.recumbents.com/wisil/whpsc2016/results.htm d http://www.whpva.org/land.html

e http://www.openfoam.com/

Acknowledgments

This work received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (grant agreement number 637972, project "ResiBots") and the German Federal Ministry of Education and Research (BMBF) under the Forschung an Fachhochschulen mit Unternehmen programme (grant agreement number 03FH012PX5 project "Aeromat").

V. Results

X up

In Figure 3 we see SAIL's prediction of the best performing designs in each region of the feature space. The values and designs shown are the median predicted performance over 20 runs of SAIL. The designs produced in these median runs are similar to those optimal designs found by CMA-ES, a selection of which are shown in black around the border. These were the best designs produced over all 20 runs of CMA-ES in each bin, using nearly 11.5 million PE to produce the whole map. It is these designs that we designate as the 'optimum' design in each bin for comparison.

Airfoil height (Z up ) is the most influential feature, with thinner airfoils typically corresponding to lower drag values (and higher performance). The location of the airfoils highest point (X up ) has a more nuanced effect, dependent on the height of the airfoil. The best performing airfoils were not those at the extremes of the feature space, as a certain lift must be maintained, thinner airfoils are not always better. Though this map is only a prediction, in the following sections we present evidence that it presents a faithful representation of the performance potential of the design space.


Vassilis Vassiliades
email: vassilis.vassiliades@inria.fr

Konstantinos Chatzilygeroudis
email: konstantinos.chatzilygeroudis@inria.fr

Jean-Baptiste Mouret
email: jean-baptiste.mouret@inria.fr

Comparing multimodal optimization and illumination

Keywords: illumination algorithms, multimodal optimization, MAP-Elites, quality diversity, behavioral diversity, novelty search

published or not. The documents may come

INTRODUCTION

Illumination [START_REF] Mouret | Illuminating search spaces by mapping elites[END_REF] or quality diversity [START_REF] Pugh | ality Diversity: A New Frontier for Evolutionary Computation[END_REF] (QD) algorithms refer to a new type of evolutionary algorithms (EAs) capable of returning a large set of solutions that are as diverse and as high-performing as possible. ese algorithms originated in the eld of evolutionary robotics with the introduction of the novelty search algorithm (NS) [START_REF] Lehman | Abandoning objectives: Evolution through the search for novelty alone[END_REF] which suggests to look for individuals that are behaviorally di erent from previously encountered ones. at is, in NS there is an explicit distinction between the genotype space (in which the EA directly operates, e.g., a space of bit strings), the phenotype space (e.g., the space of neural networks derived from the genotype space) and the behavior space (e.g., the possible behaviors of individuals over their lifetimes, such as the end-locations of robots controlled by the neural networks of the phenotype space).

NS continually explores the behavior space, without considering the task performance. However, in many cases we are o en interested in having pressure for performance (e.g., nding the Permission to make digital or hard copies of part or all of this work for personal or classroom use is granted without fee provided that copies are not made or distributed for pro t or commercial advantage and that copies bear this notice and the full citation on the rst page. Copyrights for third-party components of this work must be honored. For all other uses, contact the owner/author(s). GECCO'17, Berlin, Germany © 2017 Copyright held by the owner/author(s). 978-x-xxxx-xxxx-x/YY/MM. . . $15.00 DOI: 10.1145/nnnnnnn.nnnnnnn fastest controller that reaches a certain location and not just any controller). NS with Local Competition (NSLC) [START_REF] Lehman | Evolving a Diversity of Virtual Creatures rough Novelty Search and Local Competition[END_REF] addresses this issue using a multiobjective approach: it ranks individuals according to their novelty (as in NS) and their local performance (i.e., how many from the k closest neighbors the individual outperforms). e Multi-dimensional Archive of Phenotypic Elites (MAP-Elites) [START_REF] Cully | Robots that can adapt like animals[END_REF][START_REF] Mouret | Illuminating search spaces by mapping elites[END_REF] algorithm o ers a di erent solution: it discretizes the behavior space into a number of niches, storing in every niche only the elite individual throughout the evolutionary simulation.

Niching techniques have traditionally been used in EAs with the purpose of multimodal optimization (e.g., see [START_REF] Das | Real-parameter evolutionary multimodal optimization -a survey of the state-of-the-art[END_REF][START_REF] Harik | Finding Multimodal Solutions Using Restricted Tournament Selection[END_REF][START_REF] Mahfoud | Niching methods for genetic algorithms[END_REF][START_REF] Pétrowski | A clearing procedure as a niching method for genetic algorithms[END_REF]), i.e., for discovering the multiple optima of the underlying genotype or phenotype space. In contrast to such approaches, the primary goal of illumination algorithms is not optimization but diversity [START_REF] Pugh | ality Diversity: A New Frontier for Evolutionary Computation[END_REF]. Other di erences between the way the two approaches are typically used are the following: (1) the number of solutions returned with multimodal optimization is in the order of tens or hundreds (e.g., see [START_REF] Harik | Finding Multimodal Solutions Using Restricted Tournament Selection[END_REF]), whereas with illumination algorithms it is in the order of hundreds or thousands [START_REF] Cully | Robots that can adapt like animals[END_REF][START_REF] Pugh | ality Diversity: A New Frontier for Evolutionary Computation[END_REF]; [START_REF] Das | Real-parameter evolutionary multimodal optimization -a survey of the state-of-the-art[END_REF] illumination is performed in behavior space [START_REF] Cully | Robots that can adapt like animals[END_REF][START_REF] Pugh | ality Diversity: A New Frontier for Evolutionary Computation[END_REF], whereas multimodal optimization is performed either in genotype or phenotype space [START_REF] Das | Real-parameter evolutionary multimodal optimization -a survey of the state-of-the-art[END_REF].

It is currently unknown how multimodal optimization algorithms behave when they are set to return as many solutions as illumination algorithms and to operate in behavior space (though there are works that use speciation in behavior space [START_REF] Trujillo | Speciation in behavioral space for evolutionary robotics[END_REF]). is is what we investigate in this short paper.

Multimodal optimization has a long history and the purpose of this study is not to provide a comprehensive evaluation of the di erent algorithms that exist. Instead, we select two representative algorithms that are simple to implement and can be applied in the behavior space: the "clearing" method [START_REF] Pétrowski | A clearing procedure as a niching method for genetic algorithms[END_REF] and "restricted tournament selection" (RTS) [START_REF] Harik | Finding Multimodal Solutions Using Restricted Tournament Selection[END_REF]. Clearing requires from the user to provide a clearing radius which de nes the niche of an elite individual, whereas RTS restricts competition among a user-de ned number of randomly selected individuals from the population.

EXPERIMENTAL SETUP

We compare 2 illumination (NSLC, MAP-Elites) and 2 multimodal optimization algorithms (RTS, Clearing) in a maze navigation task (Fig. 1 upper le ) in which a simulated mobile robot (diameter: 20 units) is controlled by an arti cial neural network, whose architecture and parameters are evolved [START_REF] Mouret | Novelty-based multiobjectivization[END_REF]. e robot starts from the bo om of a maze (size: 1000 × 1000 sq. units) and needs to reach the goal point at the center.

anks to the openings, this maze permits 16 families of trajectories towards the center (thus, at least 16 behaviorally distinct optima). e tness function is the smallest Euclidean distance between the center of the robot and the goal GECCO'17, July 2017, Berlin, Germany Vassiliades et. al.

Median

Figure 1: Best distance to the goal (center) for each algorithm in the environment used during evolution (le most) and all 16 evaluation environments each of which permits a single path to the goal. e box plots show the median and the interquartile range over 30 solutions, apart from the last column which is calculated from the medians over all 17 environments. Overall, Clearing is comparable to MAP-Elites, which in turn is slightly better than NSLC. RTS has the worst overall performance.

over the robot's lifetime [START_REF] Mouret | Novelty-based multiobjectivization[END_REF], which is set to 3000 simulation steps.

e behavior descriptor of each individual (i.e., a point in behavior space) is the end location (2D) of the robot [START_REF] Lehman | Abandoning objectives: Evolution through the search for novelty alone[END_REF][START_REF] Mouret | Novelty-based multiobjectivization[END_REF].

We evaluate the quality of the solutions produced by the illumination and multimodal optimization algorithms by assessing the performance of each individual found in the nal archive/population in environments that are modi ed versions of the one used during evolution: if the archive is made of diverse and high-performing individuals, then it should contain individuals with all the kinds of trajectory, including some that work in the modi ed environments; in the extreme opposite, if all the individuals of an archive have the same behaviors, none of them will work in the modi ed environments. In our scenario, we assess whether the resulting sets contain controllers that nd all 16 paths to the goal by evaluating all the solutions produced by each algorithm in 16 di erent environments (Fig. 1, top), each of which permi ing only one path to the goal.

In addition, we measure the QD-score for each algorithm [START_REF] Pugh | ality Diversity: A New Frontier for Evolutionary Computation[END_REF], which is calculated by mapping an archive's behavior descriptors to a MAP-Elites grid, keeping the best performing one in a cell, and summing the tness scores from all cells.

RESULTS AND CONCLUSION

We use 30 independent evolutionary runs of 200k evaluations. MAP-Elites (71 2 = 5041 cells), NSLC (max archive size = 5041) and Clearing (pop. size = 5041) return solutions with a median tness of less than 10 units in all evaluation environments. RTS (pop. size = 5041) has the worst overall performance with a median distance of more than 10 units in 5 environments and more than 100 units in 2 environments. is indicates that RTS might have not found all optima, or it might have found them and eventually lost them. Interestingly, Clearing has comparable performance to MAP-Elites and NSLC and be er performance in 2 environments.

e QD-scores for a typical archive/population of all algorithms, calculated in the initial environment using a 32 × 32 grid, are the following (lower is be er): MAP-Elites: 129493.8; NSLC: 79396.1; RTS: 14504.6; Clearing: 1152. [START_REF] Lehman | Abandoning objectives: Evolution through the search for novelty alone[END_REF].

is shows that Clearing has a be er QD-score, followed by RTS, then NSLC and nally MAP-Elites. However, this does not agree with the results of Fig. 1, according to which RTS archives are substantially less diverse than MAP-Elites archives, and illustrates that the QD-score does not capture everything about behavioral diversity: Fig. 1 shows that in all the environments, MAP-Elites found at least a single, highperforming individual, whereas the QD-score shows that, in the initial environment, on average, the tness of the individual found by MAP-Elites for each bin is lower than the one found by RTS.

e similarity in performance between MAP-Elites and Clearing could be explained by the fact that they use a xed niche size, despite their di erent niche shapes (MAP-Elites: rectangular; Clearing: spherical) or whether they operate in a bounded (MAP-Elites) or unbounded (Clearing) space. Clearing has be er performance than RTS which is the opposite of the ndings in multimodal optimization [START_REF] Singh | Comparison of Multi-modal Optimization Algorithms Based on Evolutionary Algorithms[END_REF].

Overall, this study shows that when some multimodal optimization algorithms (e.g., Clearing, but not RTS) are provided with (1) large population sizes and (2) distance metrics that operate in behavior space, they can be as e ective as illumination algorithms. More investigation is warranted to understand the links between the two approaches.

ACKNOWLEDGMENTS

is work was supported by the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (Project: ResiBots, grant agreement No 637972).


Vassilis Vassiliades
email: vassilis.vassiliades@inria.fr

Konstantinos Chatzilygeroudis
email: konstantinos.chatzilygeroudis@inria.fr

Jean-Baptiste Mouret
email: jean-baptiste.mouret@inria.fr

A comparison of illumination algorithms in unbounded spaces

Keywords: •Computing methodologies → Evolutionary robotics, • eory of computation → Evolutionary algorithms, illumination algorithms, MAP-Elites, quality diversity, behavioral diversity, novelty search

Illumination algorithms are a new class of evolutionary algorithms capable of producing large archives of diverse and high-performing solutions. Examples of such algorithms include Novelty Search with Local Competition (NSLC), the Multi-dimensional Archive of Phenotypic Elites (MAP-Elites) and the newly introduced Centroidal Voronoi Tessellation (CVT) MAP-Elites. While NSLC can be used in unbounded behavioral spaces, MAP-Elites and CVT-MAP-Elites require the user to manually specify the bounds. In this study, we introduce variants of these algorithms that expand their bounds based on the discovered solutions. In addition, we introduce a novel algorithm called "Cluster-Elites" that can adapt its bounds to non-convex spaces. We compare all algorithms in a maze navigation problem and illustrate that Cluster-Elites and the expansive variants of MAP-Elites and CVT-MAP-Elites have comparable or be er performance than NSLC, MAP-Elites and CVT-MAP-Elites.

INTRODUCTION

Evolutionary illumination algorithms are a new class of algorithms capable of producing large archives of diverse and high-performing solutions [START_REF] Mouret | Illuminating search spaces by mapping elites[END_REF][START_REF] Pugh | ality Diversity: A New Frontier for Evolutionary Computation[END_REF]. Inspired by the phenomenon of species diversi cation in nature (e.g., see [START_REF] Lewitus | Natural constraints to species diversi cation[END_REF]), these algorithms have been introduced in the eld of evolutionary robotics with the purpose of encouraging diversity in what is known as the behavior space [START_REF] Lehman | Exploiting Open-Endedness to Solve Problems rough the Search for Novelty[END_REF][START_REF] Mouret | Novelty-based multiobjectivization[END_REF][START_REF] Trujillo | Discovering several robot behaviors through speciation[END_REF]. is space describes the possible behaviors of individuals over their Permission to make digital or hard copies of part or all of this work for personal or classroom use is granted without fee provided that copies are not made or distributed for pro t or commercial advantage and that copies bear this notice and the full citation on the rst page. Copyrights for third-party components of this work must be honored. For all other uses, contact the owner/author(s). GECCO'17, Berlin, Germany © 2017 Copyright held by the owner/author(s). 978-x-xxxx-xxxx-x/YY/MM. . . $15.00 DOI: 10.1145/nnnnnnn.nnnnnnn lifetimes: for example, a point in this space, i.e., a behavior characterization/signature, could be the nal positions of simulated robots whose controllers are evolved [START_REF] Lehman | Exploiting Open-Endedness to Solve Problems rough the Search for Novelty[END_REF]. In contrast, the genotype space is the space in which the evolutionary algorithm operates (e.g., a space of bit strings) and the phenotype space encodes the possible controllers (e.g., neural networks) that are derived from the genotype space.

Novelty Search (NS) [START_REF] Lehman | Exploiting Open-Endedness to Solve Problems rough the Search for Novelty[END_REF] is the rst algorithm that suggested to abandon any tness objective and continually explore for novel behaviors by de ning novelty as sparseness, i.e., the average distance to the n nearest neighbors, in behavior space. NS with Local Competition (NSLC) [START_REF] Lehman | Evolving a Diversity of Virtual Creatures rough Novelty Search and Local Competition[END_REF] improved upon NS based on the observation that it is more bene cial to explore globally and optimize locally: this local optimization is achieved using a secondary objective.

e Multi-dimensional Archive of Phenotypic Elites (MAP-Elites) [START_REF] Cully | Robots that can adapt like animals[END_REF][START_REF] Mouret | Illuminating search spaces by mapping elites[END_REF] algorithm proposed a conceptually simpler approach: it discretizes the behavior space into a grid of k cells, storing in each cell the elite solution over the evolutionary generations. is algorithm has recently been extended using a Centroidal Voronoi Tessellation (CVT) by the CVT-MAP-Elites algorithm to deal with highdimensional behavior spaces [START_REF] Vassiliades | Scaling Up MAP-Elites Using Centroidal Voronoi Tessellations[END_REF]. Intuitively, CVT-MAP-Elites partitions the behavior space by uniformly distributing k centroids: these centroids correspond to the centers of the cells in MAP-Elites if both algorithms use the same number of cells.

Both MAP-Elites and CVT-MAP-Elites assume knowledge of the bounds of the behavior space. More speci cally, they enclose the behavior space inside a bounding hyperrectangle and make the assumption that the user knows the ranges of this rectangle. In contrast, NSLC does not make such an assumption. As the original spirit of illumination algorithms is constant exploration and diversi cation, such user-de ned knowledge is a limitation of the MAP-Elites family of algorithms. us, in this study we ask two questions:

(1) Can expansive versions of MAP-Elites and CVT-Map-Elites be as e ective as their non-expansive counterparts, in spite of the fact that they do not know the bounds? (2) Would an algorithm that allocates centroids using the actual shape of the behavior space, instead of the bounding rectangle like expansive MAP-Elites, be more e ective?

NEW ALGORITHMS 2.1 Expansive MAP-Elites

In the "expansive" variant of MAP-Elites (Appendix A, Alg.1), the behavior characterizations of the o spring at every generation de ne the bounds of the space. As the newly calculated bounds change the width of the cells, solutions that already exist in the archive are taken out and treated as new solutions along with the o spring. A side-e ect of this procedure could be an initial increase in the archive size (due to lling a restricted space), with some subsequent decrease, due to the expansion of the bounds and the mapping of more than one solutions to a single cell. In contrast, in the case of standard MAP-Elites, the archive can only grow.

Expansive CVT-MAP-Elites

For the "expansive" version of CVT-MAP-Elites (Appendix A, Alg.2), we perform a similar procedure to the above. We (re)compute the CVT based on the newly-found bounds, taking out the existing solutions and treating them as new ones (because the centroids will ll a larger bounding hyperrectangle). It is worth noting that we perform this procedure periodically and not at every generation (as in expansive MAP-Elites), in order to reduce the computational load induced by the repeated CVT construction.

Cluster-Elites

In this paper, we introduce a variant of MAP-Elites and CVT-MAP-Elites called "Cluster-Elites" (Appendix A, Alg.3) that a empts to maximally spread a number of centroids on a potentially nonconvex manifold on which the behavioral descriptors reside, rather than inside a hyperrectangle de ned by the ranges of the sampled behavioral descriptors. Since the shape of this manifold is unknown and potentially high-dimensional, we cannot use approaches like alpha-shapes [START_REF] Edelsbrunner | On the shape of a set of points in the plane[END_REF] (which are generalizations of convex-hulls), but we instead resort to methods that use nearest-neighbor calculations. Clustering algorithms such as the "k-medoids" [START_REF] Kau Man | Clustering by means of Medoids[END_REF] aim to partition a dataset into k clusters by choosing the "centroid" of each cluster to be the point that minimizes the within-cluster sum of squares. Such algorithms, however, cannot be used in our case because of the problem of sampling bias: denser regions place more emphasis on allocating "resources" (centroids) there, whereas we would like to have a set of uniformly-spread points that is not highly a ected by the density of sampled solutions. In addition, sparser regions could provide the stepping stones needed for discovering be er solutions.

Cluster-Elites a empts to address these issues by continually sampling the behavior space and maximally spreading its available resources in the space spanned by the sampled solutions, while keeping in each centroid the locally est solution. More speci cally, in each generation, Cluster-Elites rst creates a set that contains a copy of the o spring and the current centroids. It then iteratively computes and removes the densest solution from the resulting set, until the size of this set reaches the desired number of well-spread centroids k, which is progressively increased over the generations. Finally, it stores at each centroid position the est solution among its local neighborhood by considering all initial solutions of the current generation (o spring and previous centroids).

EXPERIMENTAL SETUP 3.1 Task

We compare the original MAP-Elites and CVT-MAP-Elites, with their expansive variants, NSLC and Cluster-Elites. We use a maze navigation task (see Fig. 1 upper le ) where a simulated mobile robot (radius: 10 units) is controlled by an arti cial neural network, whose structure and weights are evolved [START_REF] Mouret | Novelty-based multiobjectivization[END_REF]. e robot starts from the bo om of the arena (size: 1000 × 1000 sq. units) and needs to reach the goal at the center. is arena permits 16 families of trajectories towards the goal due to the openings (thus, at least 16 behaviorally distinct optima). e tness function is the smallest Euclidean distance between the center of the robot and the goal over the robot's lifetime [START_REF] Mouret | Novelty-based multiobjectivization[END_REF], which is set to 3000 simulation steps.

e behavior characterization of each individual is the nal (x,y) location of the robot [START_REF] Lehman | Abandoning objectives: Evolution through the search for novelty alone[END_REF][START_REF] Mouret | Novelty-based multiobjectivization[END_REF].

Evaluating the quality of the archives

We evaluate the quality of the archives produced by the algorithms by measuring the performance of their solutions in 16 modi ed versions of the environment used during evolution, each of which corresponds to a di erent family of trajectories (see Fig. 1). If an archive is made of diverse and high-performing individuals, then it should contain individuals with every type of trajectory, including some that work in the modi ed environments; in the extreme opposite, if all the individuals of an archive have the same behavior, none of them will work in the modi ed environments.

RESULTS

We use 30 independent evolutionary runs of 200k evaluations (990 generations). For MAP-Elites and our expansive variant we use 71 discretization intervals per dimension (thus, 71 2 = 5041 cells), for CVT-MAP-Elites and the expansive version we use 5041 centroids, and we set the maximum archive size of NSLC to be 5041. For Cluster-Elites we use an initial number of centroids k init = 50 and increase it by adding k incr = 5 more centroids at every generation, resulting at 5k centroids at the nal generation. For the calculation of the densest points in Cluster-Elites, we empirically set the number of nearest neighbors to d + 2, where d is the dimensionality of the behavior space (i.e., d = 2).

All algorithms return solutions with a median tness of less than 10 units (radius of the robot) in all evaluation environments. In the 8 th and 14 th evaluation environments MAP-Elites, CVT-MAP-Elites and NSLC display a large variance, whereas the newly introduced algorithms have a lower median distance to the goal and less variance. e expansive variants of MAP-Elites and CVT-MAP-Elites nd "good" bounds in this environment from the 0 th generation, which become more re ned and stop changing a er approximately 100 generations. e bounds do not extend to 0 and 1000, as in the non-expansive variants, due to the outer border; this means that more cells are allocated inside the arena.

We have also calculated the quality-diversity (QD) score [START_REF] Pugh | ality Diversity: A New Frontier for Evolutionary Computation[END_REF] by mapping an archive's behavior descriptors to a 32 × 32 MAP-Elites grid, keeping the best performing one in a cell, and summing the tness scores from all cells. e QD-scores for a typical archive of all algorithms, calculated in the initial environment are the following (lower is be er): NSLC: 79396.1; Cluster-Elites: 79572.5; expansive MAP-Elites: 83770.2; expansive CVT-MAP-Elites: 94492.9; CVT-MAP-Elites: 119168.8; MAP-Elites: 129493.8.

CONCLUSION

Overall, our results illustrate that Cluster-Elites and the expansive variants of MAP-Elites and CVT-MAP-Elites have comparable performance with their " xed-bounds" counterparts and NSLC, without requiring knowledge of the bounds. Moreover, Cluster-Elites is a promising algorithm that demands further investigation. In particular, experiments with complex tasks in which the points in behavior space lie on highly non-convex manifolds could highlight the bene ts of Cluster-Elites over the bounding rectangle approach followed by MAP-Elites and CVT-MAP-Elites. In addition, combining NSLC with Cluster-Elites, i.e., by reducing the novelty archive in a manner similar to Cluster-Elites, might have advantages over both algorithms.

A PSEUDOCODES

(b min , b max ) = ([in f , ..., in f ], [-in f , ..., -in f ]) 4:
X ←random solutions() Initialization: I random x X ←random solutions() Initialization: I random x

MedianFigure 1 :

 1 Figure 1: Best performance (distance to the goal) for each algorithm in the maze environment used during evolution (le most column) and all 16 evaluation environments each of which permits a single path to the center (goal). e box plots show the median and the interquartile range over 30 solutions, apart

7 : 7 : 3 :

 773 , b max ) ←update ranges(B, b min , b max ) B ←normalize(B, b min , b max ) , b max ) ←update ranges(B , b min , b max ) 14: B ←normalize(B, b min , b max ) 15: B ←normalize(B , b min , b max ) 16: (X , P, B) ←-(X ∪ X , P ∪ P , B ∪ B ) 17: (X, P, B) ←-({}, {}, {}) Clear the archive 18: (X, P, B, X , P, B) 19:return archive (X, P, B) 20: procedure (X, P, B, X , P, B)21:for i = 0 → |X | do 22: c ←get cell index(B[i]) 23: if P[c] = null or P[c] < P[i] then 24: (P[c], X[c], B[c]) ←-(P[i], X [i], B[i]) , b max ) = ([in f , ..., in f ], [-in f , ..., -in f ]) , b max ) ←-update ranges(B, b min , b max ) C ←-CVT(k, b min , b max ) Get centroids from CVT 8: (X, P, B, X , P, B, C) , b max ) ←update ranges(B , b min , b max ) 15: C ←-CVT(k, b min , b max ) 16: (X , P, B) ←-(X, P, B) 17: (X, P, B) ←-({}, {}, {}) Clear the archive Algorithm 3 Cluster-Elites algorithm 1: procedure C E (k init , k incr , k max ) 2: k = k init , k f = 0 (X, P, B) ←create empty archive() 4:

ACKNOWLEDGMENTS

is work was supported by the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (Project: ResiBots, grant agreement No 637972).

5:

(P, B) ←evaluate(X ) I evaluations 6:

(X, P, B, X , P, B)

(X, P, B, X , P, B) Cluster-Elites procedure 8:

if k >= k max then

23: procedure (X, P, B, X , P, B)

24:

27:

else 28: return elites (X ϵ , P ϵ , B ϵ )