Premature convergence in morphology and control co-evolution: a study

2.1. Fitness sharing

This diversity preservation technique involves scaling an individual’s fitness inversely proportional to its similarity with respect to the other individuals in the population so that the fitness function is modified to penalize those individuals that are similar. In Mouret & Doncieux, (2009), similarity is defined as

ϕ ( δ ( i , j ) ) = { 1 − δ ( i , j ) σ s h if δ ( i , j ) < σ s h 0 otherwise .

(1)

Thus, the modified fitness of an individual i is defined as

f i ′ = f i c i

(2)

c i = ∑ j = 1 n ϕ ( δ ( i , j ) ) ,

(3)

2.2. Speciation

Another approach to preserving diversity in evolutionary algorithms is to divide the population into species based on the similarity between individuals, such that:

• Each species has a representative genome that comes from the previous generation.

Speciation can be carried out in combination with fitness sharing and crossing constraints, as happens in the NEAT algorithm (Stanley & Miikkulainen, 2002). However, it should be noted that it is not necessary to do so, and each one can be used as a separate technique.

2.3. Crossover operator restriction

By limiting the crossover operator, this technique allows only individuals within the same species to be crossed or, alternatively, permits the crossing of only those individuals belonging to different species.

2.4. Diversity as an objective

Using a multi-objective framework, it is possible to include a diversity metric as an additional objective to maximize. This diversity can be applied in both the variable and objective space and is typically problem-dependent.

Previous research has demonstrated that the best results are achieved by measuring the average distance d between an individual and the rest of the population. Consequently, for the purposes of this article, the diversity D of an individual x with respect to the population in generation n (P _n ) is defined as

D ( x ) = 1 | P n | ∑ y ∈ P n d ( x , y ) .

(4)

In the context of the co-evolution of morphology and control of virtual creatures, the distance between individuals can be defined in various ways depending on the representation and/or subsequent phenotypic expression they have. The choice of distance metric can have a significant impact on the co-evolutionary process and the resulting solutions. As mentioned in Mouret & Doncieux (2009), if the genotype is represented as a neural network, calculating the distance between two genotypes is an extremely complex problem since computing the distance between graphs is NP-Complete in itself. However, this problem can be reduced by choosing a vector v⁽ⁱ⁾ that expresses desirable characteristics of the network, subject to a notion of distance that is suitable for the problem at hand.

For example, in Mouret (2011), several algorithms are compared, including the integration of diversity and fitness under a multi-objective paradigm, in the task of evolving controllers (neural networks) for two-wheeled mobile robots to reach a goal in an obstacle environment.

The distance between two individuals can be defined in terms of what is called the behavior vector. For example, the behavior vector b of an individual x can defined as

b ( x ) = p r ( T ) ,

(5)

In Nordmoen et al. (2020), both morphology and control of modular robots are evolved for the task of walking the farthest distance possible. The morphology is given by a direct encoding based on trees, where each node in the tree is a module and each edge is a connection between two modules. There are two types of modules, non-mobile and joints, each supporting up to 5 and 3 child nodes, respectively. The control system consists of all individual controllers of each joint, and each controller is defined by 5 parameters. In this case, it is desired to protect morphological diversity, so the behavior vector of an individual x in the space Q of all robots that can be represented with the mentioned encoding is defined as

b ( x ) = ( m x , j x ) ,

(6)

d b ( x , y ) = 1 − e | b ( x ) − b ( y ) | .

(7)

The diversity is defined in the same way as in equation (4). It should be noted that the distance defined in this case is a function d b : R 2 → [ 0,1 ] × [ 0,1 ] ⊂ R 2 , so the multi-objective formulation is defined as

max x ∈ Q ( F ( x ) ) , where F ( x ) = [ f ( x ) D m ( x ) D j ( x ) ] .

2.5. Divergent search

The ideal of evolutionary computation to replicate the complexity and sophistication achieved in nature and use it as an optimizer is undoubtedly an extremely ambitious goal. However, this ambition has led to the search for more complex encodings and algorithms. According to this ideal, it should be enough to use selection pressure correctly toward an arbitrarily ambitious objective, but the results, especially in evolutionary robotics, tell us otherwise (Pugh et al., 2016).

Despite strong criticism from purists in the fields of machine learning and optimization, labeling it as an ad hoc method (Pugh et al., 2016), a question stands out as motivation to continue exploring the alternative of artificial evolution: how has nature been able to generate intelligence, and if we look at our best efforts to create artificial intelligent systems with traditional optimization, despite being better than humans in some tasks, they often fall short in areas where human cognition excels?

Divergent search proposes that perhaps the idea of evolution has been applied in a somewhat incorrect way. Natural evolution is divergent, rather than convergent. This idea gives rise to a new vision: evolution as a machine for diversification instead of optimization. This new perspective suggests that the problem lies in the abstraction of evolution that has been used all this time. If we look at natural evolution, it is clear that there is no optimization involved, and there is no clear objective except for a notion of reward for individuals that are not only better but also different (evolutionary niches: some organisms that are sufficiently different from their predecessors may enjoy reduced competition, making them more likely to survive).

2.6. Novelty search

It is a radicalization of the previous idea, completely disregarding optimization toward a particular goal (e.g., the fastest robot, the most efficient, the best at solving mazes), and directing the search solely through a notion of novelty as a consequence of diversity (Pugh et al., 2016). In Novelty Search, no fitness function based on a specific goal is optimized to explore the search space. Instead, a novelty metric of the individual is rewarded, allowing us to know how far it is from the rest of the population (either in the search or solution space) (Lehman & Stanley, 2011a). A novelty metric that allows us to measure the dispersion at any point in the space in question is

N ( x ) = 1 k ∑ i = 0 k d ( x , μ i ) .

(8)

This refers to the average distance of the individual x to its k nearest neighbors. Here, d(x, μ _i ) represents a certain distance metric dependent on the domain. The greater the average distance to its k nearest neighbors, the more dispersed the point x is located in a particular region (Lehman & Stanley, 2011a).

Regarding the implementation, to direct the search process toward more dispersed regions, a novelty archive Λ ⊂ Q of previously produced solutions that meet a minimum novelty threshold is usually incorporated. The novelty of each individual x ∈ P _j is evaluated by choosing its k nearest neighbors U from the solutions in the population of the j-th generation P j ⊂ Q and Λ (Omelianenko, 2019), that is, U = { μ 0 = x , μ 1 , … , μ i , … , μ k } = K N N ( x , P j ∪ Λ , d , k ) .

Novelty search has been successful in some domains, especially those in which the fitness evaluation space is considered “deceptive” (problems where objective-based evaluations fail), as objective-based evaluation is susceptible to getting trapped in local minima in which the evolutionary process converges to a single solution, losing diversity. An example of a problem where NS has been successfully applied is in the production of agents capable of navigating mazes. In particular, navigation in an old city whose map is unknown can be thought of as a deceptive problem, given that they often have an irregular path structure (Omelianenko, 2019). Novelty search has demonstrated that the evolutionary ability for diversification can be harnessed to find solutions close to global optima. In contrast to the traditional notion in evolutionary computation of survival of the fittest, novelty search emphasizes survival of the most novel (Lehman & Stanley, 2011a).

Similarly to adding diversity as an additional objective, novelty can also be incorporated as an extra objective. This was demonstrated by Mouret (Mouret, 2011). In this study, the multi-objective regimes of diversity+fitness and novelty+fitness were compared against novelty and fitness as single objectives respectively in evolving neural network controllers for mobile robots in the task of solving a maze. The best results were obtained using novelty+fitness, although only marginally better than diversity+fitness.

2.7. Quality-diversity

The problem with the Novelty Search approach is that it uses diversity as a means to achieve the usual end of finding a global minimum (Pugh et al., 2016), that is, as an optimizer. Moreover, the idea of completely eliminating objectives means that there is no notion of good and bad solutions, which seems to be crucial in most applications. In contrast to the search for global optima, Quality-Diversity (QD) methods seek to emulate the way evolution works in nature, so that in a single run, diverse ways of solving the same problem are found, all equally viable. The objective of QD algorithms is to find a maximally diverse collection of individuals with respect to the space of possible solutions, where each member is maximally fit in its own evolutionary niche (Pugh et al., 2016). QD algorithms balance the search for diversity, which finds many evolutionary niches, with local searches for quality within each niche. Within a niche, competition is focused on individuals, thus producing the best possible individuals locally within the niche. This type of search is capable of finding diverse solutions while improving existing ones (Pugh et al., 2016).

In QD, the diversity measurement of an individual x and its position in the space of features or behaviors B is given through its characterization in that space with a feature vector b(x), which represents an abstraction of the actions and/or characteristics of the phenotype of x in the context of the problem being solved. It is worth mentioning that this abstraction is equivalent to that of the behavior vector used in Mouret & Doncieux (2009). Other names for this in the QD literature are: Behavior Space and Behavior Characterization (BC) (Pugh et al., 2016), and Descriptor Space and Behavioral Descriptor (BD) (Chatzilygeroudis et al., 2020; Cully & Demiris, 2018). The choice of the features of interest is up to the user, but it should seek to be representative of the problem and how potential solutions solve it (Chatzilygeroudis et al., 2020). The idea of being able to represent solutions through a feature vector is similar to a dimensionality reduction in machine learning, in the sense of re-expressing data from a high-dimensional and difficult-to-handle space to a lower-dimensional subspace that is easier to handle. The big difference is that in this case, it is a model-free projection where the dimensions and characteristics represented in them are chosen manually, while in machine learning and deep learning, algorithms such as the well-known principal component analysis are used, which perform a linear projection from a high-dimensional space to a lower-dimensional one, or deep neural networks, which can act as feature selectors in addition to learning models. In certain models of deep learning, such as autoencoders, a pair of neural networks are used and trained to compress data from a high-dimensional feature space to a low-dimensional latent space (encoder network), and to decompress from the latent space to the original feature space (decoder network). The encoder component can be seen as a non-linear generalization of PCA, with the latent space as an approximation of the principal components (Gaier et al., 2020).

The general problem to be optimized in QD algorithms is to find, for each point b ∈ B , the parameters x with the highest fitness value (Chatzilygeroudis et al., 2020), that is

x * = arg max x F ( x ) s . t . b = b ( x ) ∀ b ∈ B .

(9)

In general, QD algorithms work by sampling with a certain granularity B , and for each sampled region, obtaining the best solutions. In these algorithms, diversity has a higher priority than quality, and therefore, a greater number of regions in the space are visited, even if they may have less potential in terms of quality.

In QD, the feature space B is divided into t niches, N = N₁, N₂, …, N _t . Where N must cover B . Each point b _j belongs to a niche N _i , according to a proximity rule. The objective of QD algorithms is to maximize a quality metric Q with respect to each niche (Pugh et al., 2016).

The promise of QD algorithms is that, in addition to being competent in performing tasks traditionally oriented toward optimization, they can surpass what is possible with traditional optimization algorithms (Pugh et al., 2016). These types of algorithms have also been called “illumination algorithms” because they reveal (or “illuminate”) the best achievable performance in each region of the phenotype space.

It should be noted that these types of algorithms are different from simply using a multi-objective framework and adding novelty or diversity as an additional objective. The distinction lies in the locality of the quality search in QD. Unlike the purely multi-objective case, where the notion of quality of an individual x is global, in QD, the notion of quality is local to the niche N _i , where the phenotype of x is according to b(x), and therefore relative to the quality of the phenotypes found so far in N _i .

2.8. Novelty Search with Local Competition

The NSLC algorithm aims to use a multi-objective framework, specifically with two objectives: quality and novelty. However, the approach to multi-objective optimization is somewhat special in this case. As explained, in QD algorithms, global optimization is not sought, as would happen in a traditional multi-objective approach. Instead, the goal is to find niches that are as diverse as possible through the novelty objective exposed by Novelty Search, and within each niche obtain high-quality solutions through the concept of “local competition,” which consists of evaluating the individual according to the number of individuals that it outperforms in terms of fitness within its own niche, that is, it locally looks for the best within the same niche. There is no promotion of a concept of “best” globally, and this algorithm does not focus on inherently better niches.

In the work of Lehman and Stanley (Lehman & Stanley, 2011b), where this algorithm was first proposed and applied to morphology and control co-evolution, Pareto-based non-dominated sorting is used with two objectives: movement quality evaluated as a certain distance traveled (in this case the Euclidean distance) by an individual in a given simulation time, and novelty evaluated as a subroutine through Novelty Search in a three-dimensional morphological space (in this case, these morphological characteristics make up a behavioral characterization) composed of height, mass, and number of active joints (actuators). The novelty metric is given by the average Euclidean distance to the k nearest neighbors (in the same way as in pure NS), and the quality objective is given by the proportion of those same k nearest neighbors that are surpassed with respect to the movement quality metric. In NSLC, niches emerge when calculating the k nearest neighbors to the individual, so it can be said that the quality of a certain individual is calculated locally to the niche as the proportion of individuals in the niche that are surpassed by it.

The niche N of an individual x is defined as the k nearest neighbors to the individual in the B space according to a distance measured with respect to b: d _b (., .), that is

N x = { x = μ 0 , μ 1 , … , μ i , … , μ k } , = K N N ( x , P j ∪ Λ , d b , k ) .

(10)

NSLC can be formulated as a two-objective multi-objective problem, as follows

max x ∈ Q { F ( x ) = [ Q ( x ) N ( x ) ] } , Q ( x ) = k − ‖ { μ i ∀ i ∈ { 0 , … , k } | μ i ∈ N x , x ∈ Q & F ( x ) < F ( μ i ) } ‖ , N ( x ) = 1 k ∑ i = 0 k d b ( x , μ i ) .

(11)

In Lehman & Stanley (2011b), a substantial trade-off was found: despite obtaining more diverse results in terms of morphology, the fittest individuals obtained by NSLC were slightly worse than the fittest individuals from a traditional evolutionary algorithm. One possible hypothesis proposed is that unlike what happens in the natural world, in artificial evolution, populations with fixed sizes are used, and this limits the ability of each niche to explore for better individuals since exploring new niches necessarily involves redirecting resources away from niches that have already been discovered. This is another example of the complexity ceiling that the field has encountered.

2.9. MAP-Elites

Unlike NSLC, where niches emerge as part of the NS search process and the novelty archive it maintains, in the MAP-Elites (Multi-dimensional Archive of Phenotypic Elites) algorithm, the behavioral space is explicitly subdivided into niches (bins/containers) created by discretizing B (Pugh et al., 2016). In this case, each bin stores the fittest individual (elite) produced so far, and unlike NSLC, the concept of population is discarded in each generation in favor of the set of elites discovered in each bin (Mouret & Clune, 2015). In this sense, it is not evolving a population, but the archive of elites (Cully & Demiris, 2018). The search process in MAP-Elites involves selecting an elite from some niche with uniform probability in each generation, and applying some variation operator (mutation or crossover) to the genotype. Afterward, the phenotype is obtained along with its corresponding feature vector, and evaluated. Then, based on each value of the phenotype’s feature vector, it is determined which niche it belongs to. Finally, the fitness of the elite in the niche to which the new individual belongs is compared. If the fitness of the new individual is greater, it is mapped to the niche as the new elite individual (Mouret & Clune, 2015). As is often the case, this process is repeated until some problem-dependent stopping criterion is met. This process can be massively parallelized (Mouret & Clune, 2015) if a pessimistic locking scheme is considered, where each process blocks the choice of the bin, the mutation, and the potential substitution in the new bin. Batch processing can also be considered, where parallelism is implicit, similar to that found in population-based evolutionary algorithms.

It is noteworthy that in Chatzilygeroudis et al. (2020 and Cully & Demiris (2018), the authors propose a framework for the generalization of QD algorithms. This framework abstracts the main components common to all QD algorithms and provides an example of how it can be used to generate new QD algorithms. Additionally, the authors distinguish between two main types of archives present in QD algorithms: the NSLC archive, referred to as an unstructured archive or simply an archive, and the MAP-Elites archive, referred to as a “grid.”

2.10. Multi-BC QD

One important aspect of a suitable feature space is its alignment with the concept of quality. Alignment in the context of a BC is defined as the degree to which finding diversity also leads to finding higher fitness (Pugh et al., 2016). It has been found that seeking diversity using a BC with higher alignment results in better outcomes simply because the search for diverse solutions also leads to better quality solutions, thereby mitigating the problem of premature convergence. Two examples where the effects of having aligned and unaligned BCs can be distinguished are:

• In the domain of maze-solving agents, the goal is to find the topology and weights of a neural network that controls the agent through inputs provided by proximity sensors (for obstacle detection) and navigation sensors (that inform the robot the direction of the goal) and outputs that control linear and angular acceleration. An associated BC simply based on the final position of the agent in the maze will be highly aligned, since an agent that continuously finds new positions will eventually reach the final position in the maze (Pugh et al., 2016).

The degree of alignment is extremely difficult to measure a priori; however, it can be anticipated by considering the following two properties:

1. Each behavior is associated with a reduced range of fitness values.

In Pugh et al. (2016), it is suggested that when diversity is sought using unaligned BCs, QD algorithms fail to obtain optimal solutions, and in more complex problems, they may completely fail to make any progress toward even viable solutions. This is why they propose the use of new variants of MAP-Elites and NSLC algorithms capable of handling multiple BCs. The following are the most representative variants.

3.1. Representation

We chose to adopt Lipson et al.’s approach (Lipson et al., 2016) as the basis for the morphology and control representation in our study. This approach employs a diploid genotype that encodes both morphology and control using two CPPN (Compositional Pattern Producing Network) neural networks (Stanley, 2007). These networks are structured as weighted, connected, acyclic graphs, where each node or neuron corresponds to a perceptron (Mcculloch & Pitts, 1943).

The use of two different CPPNs allows for variations in one aspect without affecting the other, direct measurement of genetic diversity and complexity of the neural networks of each aspect, and proposing relevant indicators to detect and measure the phenomenon of premature convergence of morphology relative to control. Additionally, as already observed in the literature (Cheney et al., 2014), the use of CPPNs as an indirect encoding tends to lead to better evolvability, as they are capable of producing morphologies and controllers with regular patterns.

For the CPPNs used, we have a region Ω ⊂ R 3 which refers to an arbitrarily bounded three-dimensional design space, which in our case is

Ω = { x ∈ [ 0,1 ] y ∈ [ 0,1 ] z ∈ [ 0,1 ] .

(12)

3.1.1. Morphology

For the experiments conducted in this article, we chose to evolve robots with soft voxel-based morphologies (softbots). We chose these morphologies for their promise in terms of versatility, ease of evolution, and potential to be transferred from virtual to physical worlds. We used morphologies with two types of voxels, active and passive. The active voxels act as muscles or actuators, each performing expansions and contractions at the same frequency but with different phase offsets, allowing for movements similar to the propagation of a wave. A valid morphology is created in a discrete design space with a fixed resolution of 5 × 5 × 5, and is composed of a combination of any number of contiguous voxels.

The neural network CPPN _m has a fixed number of inputs and outputs. It has four inputs (x, y, z, r), corresponding to the three Euclidean coordinates of the three-dimensional design space, plus a radial coordinate (with the hope of encouraging the appearance of radial patterns). Similarly, it has two outputs, ( p , m ) ∈ R 2 , where p indicates the presence or absence of the voxel in the morphology, and m indicates the material (active or passive). The way in which a softbot s morphology is “grown” from its CPPN is by making individual queries to the neural network for each element of the design space. If p ∈ (−∞, 0), the voxel remains empty, and if p ∈ [0, ∞), the voxel belongs to the morphology and is either passive if m ∈ (−∞, 0) or active if m ∈ [0, ∞).

We can say that the CPPNs in these experiments are functions of (x, y, z), and not r, by considering as always having an input neuron with and activation function r ( x , y , z ) = x 2 + y 2 + z 2 .

Henceforth, the evaluation of CPPN _m can be expressed as a function m → s : Ω → R 2

m → s ( x , y , z ) = C P P N m ( x , y , z ) = [ p ( x , y , z ) m ( x , y , z ) ] .

(13)

Which furthermore, gives place to the expression of the morphology as

M s ( x , y , z ) = { V o x e l n i l if p ∈ ( − ∞ , 0 ) V o x e l p if p ∈ [ 0 , ∞ ) & m ∈ ( − ∞ , 0 ) V o x e l a if p ∈ [ 0 , ∞ ) & m ∈ [ 0 , ∞ ) ( x , y , z ) ∈ Ω .

(14)

3.1.2. Control

The control system is composed of all the active voxels within the morphology, which is defined as

A ( s ) = { M s ( x , y , z ) | ( x , y , z ) ∈ Ω & M s ( x , y , z ) = V o x e l a }

(15)

Once the morphology is obtained, all the control system’s voxels are taken, and given the coordinates of each one Ω A = { ( x , y , z ) | ( x , y , z ) ∈ Ω y M s ( x , y , z ) = V o x e l a } , a second network CPPN _c corresponding to the control is queried. The evaluation CPPN _c (x, y, z) yields the control parameters for each active voxel.

In CPPN _c , the inputs are the same as in CPPN _m , and the outputs are also two, (ϕ, τ) ∈ [0, 2π] × [−1, 1], where ϕ corresponds to the phase shift that the actuator voxel has with respect to a global frequency f _g , and τ corresponds to the period. For the case of ϕ, it can take scaled values of [0, 2π] radians, which is achieved by always connecting the output u to a neuron with a logistic activation function scaled by 2π

ϕ ( u ( x , y , z ) ) = 2 π 1 + e − u ( x , y , z ) .

(16)

The period τ is used to obtain the global frequency f _g , averaging the frequency obtained for each active voxel f (x, y, z) = 1/τ(x, y, z), ( x , y , z ) ∈ Ω A

f g = 10 1 ‖ Ω A ‖ ∑ ( x , y , z ) ∈ Ω A f ( x , y , z ) + 1.5 .

(17)

The value of τ is obtained by connecting the output v of the neural network to a logistic activation function biased in the interval [−1, 1], where v ( x , y , z ) : Ω A → R . The evaluation of τ is defined as follows

τ ( v ( x , y , z ) ) = 2 1 + e − v ( x , y , z ) − 1 .

(18)

Moreover, the calculation of f _g has discontinuities at 1 / ‖ Ω A ‖ ∑ ( x , y , z ) ∈ Ω A f ( x , y , z ) = − 2 / 3 and τ(v (x, y, z)) = 0. Therefore, all softbots whose f _g calculation leads to numerical errors of divisions by zero will be evaluated with a fitness of zero.

Henceforth, the evaluation of CPPN _c can be expressed as a function c → s : Ω A → { [ 0,2 π ] × [ − 1,1 ] }

c → s ( x , y , z ) = C P P N c ( x , y , z ) = [ ϕ ( x , y , z ) τ ( x , y , z ) ] .

(19)

Which gives place to the expression of the controller’s parameters as

C s ( x , y , z ) = [ f = f g ϕ = ϕ ( x , y , z ) ] , ( x , y , z ) ∈ Ω A .

(20)

The control system is like a mass-spring-damper system, where masses correspond to the centroids of voxels (nodes), and springs to the edges formed from the Delaunay triangulation, thereby the resting length L₀ is determined. Actuation is achieved by altering the resting length according to the following signal

L ( t ) = [ 1 + A sin ( 2 π t τ 1 + ϕ 1 ) + A sin ( 2 π t τ 2 + ϕ 2 ) ] ⋅ L 0 .

where τ₁, τ₂ and ϕ₁, ϕ₂ are the periods and phases associated with the centroids at the end of the spring. In other words, modifying the resting length serves as an input signal to the system. The controllers are open-loop, as they do not take into account any error signals or sensor feedback.

3.1.3. Softbot definition

Per the definitions for morphology and control presented in Sections 3.1.1 and 3.1.2, it is possible to define a softbot s, such as the ones evolved in these experiments, based on its morphology (Equation (13)) and control (Equation (19)) genotypes, as well as its morphology (Equation (14)) and control (Equation (20)) phenotypes.

The genotype G _s of a softbot s is defined as

G s = ( C P P N m , C P P N c , m → s ( x , y , z ) , c → s ( x , y , z ) ) .

(21)

And the phenotype P _s of s, is defined as

P s = ( M s ( x , y , z ) , C s ( x , y , z ) ) .

(22)

With which, a softbot s is defined as

s = ( G s , P s ) .

(23)

3.2. Task

The task chosen for all experiments in this work was locomotion in a straight line, undoubtedly the most common task in the literature on morphology and control co-evolution. However, it is not the only one, and in practice, locomotion toward an objective (Lan et al., 2018; Lan et al., 2021; Sims, 1994;), locomotion with obstacles and in different environments (Stensby et al., 2021), moving in confined spaces (Cheney et al., 2015), or competing to seize an object (Miconi & Channon, 2006; Sims, 1994) may be more of interest. Nonetheless, the chosen task is comparatively simple and reduced, and the interest of this research is not to explore different environments or tasks (although they are equally valid research avenues), but mainly to investigate the incidence of premature convergence of morphology with respect to control, and to benchmark different algorithms with mechanisms for protecting morphological diversity according to relevant indicators.

3.3. Simulated environment

The softbots were tested on a flat surface within the Voxelyze physics engine (Hiller & Lipson, 2014), which was accelerated by GPU using CUDA, and referred to as Voxcraft (Liu et al., 2020).

3.4. Algorithms and parameters

Firstly, we shall define the unrestricted softbot search space S ′ . The design space, per its resolution, is composed of 5 × 5 × 5 voxels, and hence has a volume of 5³ = 125 voxels. The space S ′ is composed of all softbots which can be expressed in that design space. In order to be able to define the actual search space which was explored in the experiments S , a restriction to S ′ was put in place in order to avoid trivial morphologies.

With such restrictions, S is defined as

S = { s | s ∈ S ′ y 7 ≤ M a ( s ) + M p ( s ) ≤ 112 ⌈ 0.1 ⋅ ( M a ( s ) + M p ( s ) ) ⌉ ≤ M a ( s ) ≤ ⌊ 0.9 ⋅ ( M a ( s ) + M p ( s ) ) ⌋ } ,

(24)

M a ( s ) = ‖ A ( s ) ‖ , M p ( s ) = ‖ P ( s ) ‖ , M p ( s ) , M a ( s ) : S ′ → [ 0,125 ] ⊂ N , S ⊂ S ′ .

(25)

Functions P ( s ) and A ( s ) , are defined as

A ( s ) = { s . P s . M s ( x , y , z ) | ( x , y , z ) ∈ Ω & s . P s . M s ( x , y , z ) = V o x e l a } P ( s ) = { s . P s . M s ( x , y , z ) | ( x , y , z ) ∈ Ω & s . P s . M s ( x , y , z ) = V o x e l p } .

(26)

The experiments were conducted using five different algorithms: (i) Single-Objective Evolutionary Algorithm (SO), (ii) Fitness-Novelty MOEA (QN), (iii) Novelty Search with Local Competition (NSLC), (iv) Multi-dimensional Archive of Phenotypic Elites (MAP-Elites), and (v) Multi-BC NSLC (MNSLC), with 20 runs of each algorithm to obtain statistically significant results.

To identify the parameters with which the algorithms had to be run, especially those with diversity protection mechanisms, runs were performed on benchmark problems with multimodal functions instead of directly on the morphology and control co-evolution problem due to its high computational cost. It was of particular interest to find a set of parameters that balanced the fact that a higher number of individuals per generation would take advantage of GPU simulations. The QN, NSLC, and MAP-Elites algorithms were run on the Rastrigin, Ackley, Eggholder, and the fourth Schaffer Function benchmark problems, paying special attention to the parameters used for the latter, as the domain for which it is defined (x₁, x₂ ∈ [−100, 100]) is similar in size to the misaligned BC used for the number of passive and active voxels M _p , M _a ∈ [0, 125]. Once the parameters for these functions were identified, some test runs were done with fewer generations (200) to further refine the obtained parameters. By doing all this, the following parameters were established for all experiments:

• population_size: 100

Two archives and two grids are used in the experiment. The two archives are employed to measure aligned and unaligned novelty. The two fitness and aligned novelty grids are used to assess coverage and QD-scores. In the case of the fitness grid, MAP-Elites also employs it for the search process instead of a population.

3.5. Implementation

PyMOO (Blank & Deb, 2020), a multi-objective evolutionary algorithms framework implemented in Python, was used due to its strong object-oriented approach, modularity, and extensibility, allowing for clean, cohesive, and simple integration of the specific implementations of Novelty Search, Novelty Search with Local Competition, and MAP-Elites into a generic evolutionary algorithm, as well as the variations of these algorithms, and the evaluation function based on Voxcraft simulations.

3.5.1. Single-Objective Evolutionary Algorithm (SO)

We took as a baseline the simple single-objective evolutionary algorithm (μ + λ). The problem to be solved was

max s ∈ S ( F ( s ) = ‖ r → s ( t = T ) ‖ ) , where F : S → R

(27)

This algorithm does not provide any kind of protection for morphological diversity and serves as a starting point to identify the phenomenon of premature convergence of morphology with respect to control and compare it with other algorithms that do provide protection for morphological diversity.

3.5.2. Fitness-Novelty MOEA (QN)

In this experiment, the problem to be solved was

max s ∈ S ( F ( s ) ) , where F ( s ) = [ F ( s ) N u ( s ) ] and F : S → R 2 .

(28)

The NSGA-II implementation from the PyMOO framework was used.

We aimed to answer the following questions using this algorithm:

• How does adding morphological novelty as a second objective affect the phenomenon of premature convergence of morphology with respect to control?

• How does adding morphological novelty as a second objective affect the fitness of the solutions?

• How does this multi-objective approach that uses global fitness compare to its QD counterparts that consider local fitness (NSLC and MNSLC) in terms of the phenomenon of premature convergence and the fitness of the solutions?

3.5.3. Novelty search with local competition

This is the same algorithm as the one described in Section 2.8, the problem to be solved was

max s ∈ S { F ( s ) = [ Q ( s ) N u ( s ) ] } , where F : S → N × R .

(29)

The NSGA-II implementation from the PyMOO framework was used, replacing the crowding distance with the genetic diversity calculation described in equation (54).

The purpose of including this algorithm in the experiments is to provide a response to the following questions:

• How does using local competition instead of global fitness affect the phenomenon of premature convergence of morphology compared to the control?

• How does using local competition instead of global fitness affect the fitness of solutions?

• How does using genetic diversity instead of crowding distance affect the results?

• How does it compare to a QD algorithm with fitness-based selection pressure (ME) and another Multi-BC algorithm in terms of the phenomenon of premature convergence and the fitness of solutions?

3.5.4. Multi-dimensional Archive of Phenotypic Elites (MAP-Elites)

MAP-Elites was implemented using the generic genetic algorithm proposed in Eiben and Smith (2003), taking into account the following considerations:

1. Algorithm Initialization: The initialization of MAP-Elites relies on the discretization of the feature space into bins or niches and the initialization of a number of individuals, denoted as P, which can be considered as the initial population. Each individual is placed in its corresponding bin based on b(x).

• Discretization of B : Given a feature space B of dimension n, defined as

B = { b ∈ R n | l i ≤ b i ≤ u i ∀ i ∈ { 1 , … , n } } ,

(30)

where l _i and u _i are the lower and upper limits of the i-th feature, defining the possible values for each feature. The discretization b _i of a certain feature b _i , with a certain granularity or number of points q _i , is given by

b i = { l i + k i Δ b i | k i ∈ { 0 , … , q i − 1 } } ,

(31)

where : Δ b i = u i − l i q i − 1 .

The discretization B of the space B is defined as follows

B = b 1 × b 2 ⋯ × b i ⋯ × b n

• Population initialization and placement of individuals in the discretized space B is done as follows: A random sample of individuals of size P is taken. P genotypes are initialized randomly, expressed as phenotypes, and evaluated. In order to place an individual x in the discretization B, the feature vector b(x) is taken into consideration. The bin in which x should be placed is determined based on the range of values to which it belongs in each characteristic. The number of bins or subdivisions per characteristic is determined by s _i = q _i − 1. Therefore, we have

Δ b i = u i − l i s i .

(32)

If we want the index of the bin where an individual x falls in the i-th dimension, given the value of b(x) _i , we can simply solve for k _i in the resulting expression obtained by substituting equation (32) into equation (31). Considering that the bin indices start at zero, the expression that determines the index of the bin in the i-th dimension where the individual x is located, given its feature vector b(x), is as follows

J i = min ( ⌊ ( b ( x ) i − l i ) s i u i − l i ⌋ , s i − 1 ) .

(33)

With this, it is possible to convert a real-valued feature vector into an integer-valued feature vector, where each entry determines the bin or stencil it belongs to, based on the value of the corresponding feature. In other words

( b 1 , b 2 , … , b i , … , b n ) → ( J 0 , J 1 , … , J i , … , J n ) .

With this mapping, we can obtain the index of the bin to which the individual x belongs by

I b ( x ) = J 0 + ∑ j = 1 n ( J j ∏ i = 0 j − 1 s i ) .

(34)

2. Survival: Each individual in a population of offspring of size P is evaluated and placed in a bin based on its feature vector b. If the offspring is fitter than the elite of the bin or the bin is empty, the offspring becomes the new elite.

3. Selection: In the original algorithm, each individual in a parent population of size P is chosen with uniform probability from the elites present in the grid. In these experiments, we made the modification described in (Cully & Demiris, 2018) to select parents proportionally to their fitness.

The modification in the selection process was made in order to investigate the following research questions:

• How does the introduction of global selection pressure affect a QD algorithm like MAP-Elites?

• Does it still maintain the desirable characteristics of QD algorithms?

• How does it compare to a Multi-BC QD algorithm?

3.5.5. Multi-BC NSLC (MNSLC)

It’s the same algorithm described in section 2.10.1. The problem being addressed is the following

max s ∈ S { F ( s ) = [ Q ( s ) N u ( s ) N a ( s ) ] } , where F : S → N × R 2 .

(35)

The behavior or aligned novelty function N _a (s) is described in detail later in Section 4.1.2. With the introduction of this algorithm to the experiments, we ask the following two research questions:

• Can the addition of a third novelty objective to NSLC in a BC that is assumed to be aligned with quality mitigate the problem of slow progress in QD algorithms when the diversity of interest (morphological) is misaligned with quality?

• Are any desirable characteristics of QD algorithms lost, especially those related to premature convergence of morphology with respect to control?

4.1. Task performance indicators

4.1.2. Behavioral or aligned novelty

A notion of diversity that appears to be quite aligned with fitness or quality is defined based on the three-dimensional simulation space in which the softbots perform, since it would be precisely measured in the feature space (BC) where the robots’ fitness is measured.

Drawing direct inspiration from (Mouret, 2011), we define the behavioral distance between two softbots s₁ and s₂ as the L₂ norm of the vector formed by their final positions after a simulation of duration T, that is

d r → ( s 1 , s 2 ) = ‖ r → s 1 ( T ) − r → s 2 ( T ) ‖ , d r → : S 2 → R ( + ) .

(36)

Using the aforementioned distance, it is possible to define behavior novelty, or alternatively, aligned novelty (given that it is defined based on the distance over an aligned BC). To calculate the aligned novelty N _a of a softbot s an aligned novelty archive is required. By using equation (36) and substituting it into 8, we obtain

N a ( s ) = 1 k ∑ i = 0 k d r → ( s , μ i ) , N a : S → R ( + ) .

(37)

4.1.3. QD indicators

In general, the primary metric for QD algorithms is the QD indicator. This indicator uses a grid, and regardless of the quality metric (Q) used in the selection criteria for elites in each bin, it is summed for all bins, resulting in an indicator that can inform the total or accumulated local quality of the solutions found so far or the quality of the diversity found so far.

Let Q (x _i ) be the quality of individual x _i (or any attribute of x _i ) which has the highest quality found so far in bin i of the grid G, where G(i) = x _i . The QD indicator is defined as

Q D − s c o r e = ∑ i = 0 | G | − 1 Q ( x i = G ( i ) ) x i ∈ Q .

(38)

It is worth mentioning that this kind of indicator, along with coverage (which will be discussed later), is always computed with respect to an archive (grid) since they are derived from the MAP-Elites algorithm. This is different from all the indicators mentioned earlier, which were based on the population. It is important to keep in mind that in MAP-Elites, the primary structure being evolved is the grid, not the population as in other algorithms.

We propose two QD indicators: QD _ff indicator (Fitness Grid/Fitness) and QD _anan indicator (Aligned Novelty Grid/Aligned Novelty).

For the indicator QD _ff , a grid is introduced where the main criterion for retaining the elite in each bin is fitness. In each evaluated generation, a MAP-Elites selection cycle is performed, and quality is obtained by substituting Q _i (x) = F(x) in (38)

Q D f f = ∑ i = 0 | G f | − 1 F ( s i = G f ( i ) ) s i ∈ S .

(39)

On the other hand, for the indicator QD _anan , a grid based on the proposed aligned behavioral characterization (BC) is introduced, where the selection criterion for each MAP-Elites cycle was the aligned novelty N _a instead of fitness F. Similar to the previous indicator, the notion of quality for each bin, in this case, aligned novelty, are summed up to give a metric that tells us the quality found in terms of morphological diversity in the solutions of an algorithm, where quality is understood as a notion of diversity supposedly aligned with quality

Q D a n = ∑ i = 0 | G a n | − 1 N a ( s i = G a n ( i ) ) s i ∈ S .

(40)

It should be noted that the aligned novelty in this grid is calculated not only with respect to the members of the grid and the current population but also with respect to the members of the aligned novelty archive. This is because the archive maintains a history of the solutions generated throughout the evolutionary run, not just the elites. If only the elites are considered against the current population, members of the current generation could be mistakenly considered more novel than they actually are, and members of the grid could be considered less novel, or vice versa.

4.2. Phenotypic diversity indicators

4.2.1. Morphological diversity

If we want to introduce diversity in the co-evolution of morphology and control of virtual creatures, we need a notion of distance, which can vary greatly. Furthermore, if we decide to measure distances between the neural networks that encode morphology (graphs), and we use a graph-defined distance, such as an edit distance, we would find that such a distance for directed and cyclic graphs is an NP-Complete problem in itself. However, the notion of vector and feature space (BC) can be introduced by making a model-free dimensionality reduction. We took inspiration from Nordmoen et al. (2020) to define the morphological feature vector of a softbot s as

b m ( s ) = ( M p ( s ) , M a ( s ) ) , b m : S → N 2 .

(41)

where M _p and M _a are the number of passive and active voxels, respectively, in the morphology of the softbot s. With this, the morphological distance between two softbots s₁ and s₂ is defined as

d m ( s 1 , s 2 ) = ‖ b m ( s 1 ) − b m ( s 2 ) ‖ , b m : S 2 → R ( + ) .

(42)

Thus, the morphological diversity of a softbot s with respect to the population in generation n (P _n ) can be defined by substituting equation (42) into equation (4)

D m ( s ) = 1 | P n | ∑ s ′ ∈ P n d m ( s , s ′ ) , D m : S 2 → R ( + ) .

(43)

4.2.2. Morphological or unaligned novelty

In morphology and control co-evolution, the search for novel morphologies often does not lead to better walking patterns. Therefore, the BC defined earlier, on which the morphological distance was defined, is a misaligned BC. Similarly, if we define novelty using this distance, we will have a morphological or unaligned novelty. To calculate this unaligned novelty N _u of a softbot s, we use equation (42) and substitute it into equation (8)

N u ( s ) = 1 k ∑ i = 0 k d m ( s , μ i ) , N u : S → R ( + ) .

(44)

4.2.3. QD _fun indicator (Fitness Grid/Unaligned Novelty)

The QD _fun indicator quantifies the accumulated unaligned novelty present in the fitness grid. It is computed by setting Q _i (x) = N _u (x) in equation (38), and is evaluated with respect to the fitness grid, population, and the archive of unaligned novelties

Q D f u n = ∑ i = 0 | G f | − 1 N u ( s i = G f ( i ) ) , s i ∈ S .

(45)

4.2.4. Coverage

This metric is defined as the ratio between the number of discovered cells and the total number of cells in the grid. It can be expressed as

C = | { G ( i ) ∀ i ∈ { 0 , … , | G | − 1 } | G ( i ) ≠ n i l } | | G | .

(46)

It is worth noting that for a given algorithm, any number of grids defined in the same BC will always have the same coverage, as it quantifies the proportion of the BC space that has been visited at least once, regardless of the quality metric being used to update the grid.

4.3. Genetic diversity indicators

These indicators arise in the implementation of NSLC (Novelty Search with Local Competition), as there is a need to define a more rigorous diversity measure among genotypes than what is originally defined in the original paper of the algorithm (Lehman & Stanley, 2011b).

This indicator is derived from the implementation of NSLC, as there is a need to define a more rigorous diversity among genotypes compared to the original paper of the algorithm (Lehman & Stanley, 2011b). In that work, it is stated that there is a need to replace the crowding distance of NSGA-II as a mechanism to promote diversity along Pareto fronts with a genotype-based diversity indicator. The authors argue that a diversity defined in the space of local competition and novelty objectives may not necessarily reflect the diversity of interest (in this case, the diversity of morphologies). For instance, it is possible to have two individuals with exactly the same local competition and novelty scores. However, since these scores are entirely relative to the k nearest individuals, these two individuals could be located in completely different regions of the morphology space (i.e., they could have distinct morphologies). However, the crowding distance would penalize these two individuals heavily, potentially discarding them. Lehman and Stanley (Lehman & Stanley, 2011b) then proceed to define a genetic diversity based on a distance metric that counts the nodes and edges of the graphs that represent the morphology and control of the evolving robots. For the genotype distance defined in this article, the fact that we are working with genotype representations based on neural networks was leveraged to define a genotype distance based on the outputs of the CPPNs that encode the morphology and control. Considering that a neural network is a universal function approximator (Cybenko, 1989), a CPPN with m input neurons and n output neurons can be seen as a function

F : R m → R n .

(47)

〈 f | g 〉 | a b = ∫ a b f ( x ) g ( x ) d x .

(48)

the distance between f and g can be seen as d (f, g) = ‖f − g‖.

If the norm of the function f is defined in terms of the inner product as

‖ f ‖ = 〈 f | f 〉 1 2 .

(49)

‖ f − g ‖ = 〈 f − g | f − g 〉 1 2 = ∫ a b ( f − g ) ( f − g ) d x = ∫ a b ( f ( x ) − g ( x ) ) 2 d x .

(50)

Extending to vector fields F , G : R m → R n , we can define a distance between vector fields over a region of space Ω ∈ R m as

‖ F − G ‖ = 〈 F − G | F − G 〉 1 2 = ∫ Ω ( F − G ) T ( F − G ) d Ω .

(51)

For the CPPNs used in the experiments, the region Ω corresponds to the design space. For a design space with a resolution of R = X × Y × Z, we have that

Δ Ω i = b 1 − a 1 X b 2 − a 2 Y b 3 − a 3 Z .

(52)

As a result, the previously defined distance can be expressed in discrete form as

d ( F , G ) ≈ ∑ i = 0 R − 1 ( F ( x i , y i , z i ) − G ( x i , y i , z i ) ) T ( F ( x i , y i , z i ) − G ( x i , y i , z i ) ) Δ Ω i , ( x i , y i , z i ) ∈ Ω .

(53)

4.3.1. Genetic diversity

If we use the morphology and control neural networks as the gene expression function g → s for a softbot s, where g → s : R 3 → R 4 , based on equation (53), we can define the genetic diversity of s with respect to the population in generation n (P _n ) by substituting equation (53) into equation (4) as

D g ( s ) = 1 | P n | ∑ s ′ ∈ P n d ( g → s , g → s ′ ) , D g : S 2 → R ( + ) .

(54)

4.3.2. Genetic diversity of control

If we now consider only the neural network controller as the gene expression function for the controller of a softbot s, denoted as c → s : R 3 → R 2 , we can apply a similar approach as in equation (54)

D g c ( s ) = 1 | P n | ∑ s ′ ∈ P n d ( c → s , c → s ′ ) , D g c : S 2 → R ( + ) .

(55)

4.3.3. Genetic diversity of morphology

If we now do the same but instead we take the morphology neural network of a softbot s, we will have that

D g m ( s ) = 1 | P n | ∑ s ′ ∈ P n d ( m → s , m → s ′ ) , D g m : S 2 → R ( + ) .

(56)

5.1. Task performance

5.2. Phenotypic diversity

5.2.1. Morphological diversity

The impact of the loss of morphological diversity becomes evident when comparing algorithms where explicit diversity preservation is absent (SO) with those where such strategies are present (QN, ME, MNSLC, NSLC), as shown in Figure 2(a). The differences between these situations are clear and contrasting based on this indicator. The information presented in Figure 2(a) reveals the following:

• The SO algorithm exhibits a decreasing trend that reaches a stable state at much lower values.

• In ME, a decreasing trend can be observed in the first few generations, which could be attributed to the fact that much of the morphological search space has not yet been explored in these early stages of optimization. However, this behavior ceases after a few generations, and a slow increase in morphological diversity can be observed, stabilizing at a value of approximately 35.

• A very rapid growth pattern can be observed in all algorithms that use NSGA-II and explicitly target N _u

• In QN, a higher variance can be distinguished, as well as a decreasing pattern immediately after the initial period of rapid growth. This behavior is more pronounced than in NSLC and MNSLC, almost reaching the same value at which ME stabilizes.

• For the QD algorithms with a local selection pressure (NSLC and MNSLC), a slightly decreasing trend can be observed, which stops almost immediately and quickly stabilizes within a range of values.

OPEN IN VIEWER

Figure 2.

Convergence, KDE, and Box plots for phenotypic diversity indicators: (a)-(c) Morphological diversity D _m , (d)-(f) Unaligned novelty N _u , (g)-(i) Fitness Grid/Unaligned Novelty Q D f u n , and (j)-(l) Coverage C .

In broad terms, the algorithms that exhibit the lowest loss of morphological diversity are the “purely” QD ones (NSLC and MNSLC), followed by QN, which integrates global quality rather than local quality and aligned novelty (or morphology) into a multi-objective framework. Next is ME, which does not include any objective for aligned novelty or explicit selection pressure toward it, but rather finds morphological diversity through a grid defined in the aligned behavior space and by maximizing Q D f f . In ME, the selection pressure toward global fitness becomes evident, with the loss of population-level diversity it represents compared to its “pure” QD counterparts (NSLC and MNSLC).

5.3. Genetic diversity indicators

5.3.1. Genetic diversity D _g

When observing the behavior of this indicator in Figure 3(a)-(c) and comparing it with other diversity indicators, there seems to be a certain relationship between genetic diversity as defined in this article, and phenotype diversity. Regarding this indicator, each experiment presents a characteristic and easily distinguishable behavior at first glance. As shown in Figure 3(a), no experiment exhibits sudden growth at the beginning of the evolutionary runs, that is, in all experiments, there is some degree of loss of genetic diversity, a phenomenon that is prevalent in the literature. With regard to Figure 3(a), it can also be observed that:

• Starting with SO, a rapid loss of genetic diversity is observed, closely aligned with all the previously mentioned population-level phenotypic diversity indicators (N _u and D _m ).

• QN exhibits the second-largest loss of genetic diversity, while ME exhibits the third-largest. This is likely due to the use of a grid in the morphological feature space, which is closely related to the behavior found in population-level phenotypic diversity indicators. Once again, it appears that the introduction of selection pressure proportional to global quality separates ME from the other QD algorithms.

• In the case of NSLC and MNSLC, a slightly different behavior is observed, with one significantly outperforming the other. In this case, MNSLC outperforms NSLC, which can perhaps be explained by the fact that aligned novelty search leads to greater diversity in the controllers. Additionally, a slightly increasing behavior can be observed after the initial loss of diversity.

OPEN IN VIEWER

Figure 3.

Convergence, KDE, and Box plots for genetic diversity indicators: (a)-(c) Genetic diversity D _g , (d)-(f) Genetic diversity of control D _gc , and (g)-(i) Genetic diversity of morphology D _gm .

It is noteworthy that the two experiments that exhibit local competition (or quality) are also the ones that exhibit the highest genetic diversity. It also seems that QD algorithms (ME, NSLC, and MNSLC) are the ones that exhibit the highest genetic diversity. Furthermore, it appears that NSLC and MNSLC, which represent more purely the characteristics of QD algorithms, such as division into niches (unstructured in this case) in the feature space and the notion of local quality, are the ones that exhibit the highest genetic diversity. On the other hand, ME, by favoring a selection pressure proportional to global quality instead of a uniform selection, loses genetic diversity, as was the case with N _u and D _m .

5.4. Results of other indicators and visualizations

In this section, we present other indicators and visualizations in addition to those presented previously. Although they are not relevant for the selection of the winning algorithm, they do represent interesting results.

6.1. Hypothesis test

Since the data does not meet the assumptions for parametric statistics, namely, normality, we used the non-parametric Wilcoxon Rank-Sum test for conducting pairwise comparisons. Additionally, considering the high computational cost of the experiments, we avoided conducting the full set of experiments for each hypothesis (5 ⋅ 4 ⋅ 11 = 220), and instead, we only executed twenty trials per experiment. Therefore, the Wilcoxon Rank-Sum test was performed with Bonferroni correction. We choose α = .05/220 = .000227 as the minimum p − value threshold to reject the null hypothesis.

6.2. Expressing the preferences

With a method in place to compare algorithms, we can arrange them in order based on this comparison. In the simplest case, we have 2 algorithms, 1 problem, and 1 indicator.

• Perform multiple independent runs

– A 1 ≺ A 2 . I f μ 1 > μ 2 ; p < α

– A 2 ≺ A 1 . I f μ 2 > μ 1 ; p < α

– A 1 ∼ A 2 . I f p ≥ α

A matrix can be constructed to express the preferences of (P₁, I₁) for n candidates (algorithms) A = A 1 , … , A n , where P₁ is a problem and I₁ is an indicator:

OPEN IN VIEWER

	A 1	⋯	A n
A 1	–	↑	↓
⋮	↓	–	↓
A n	↑	↔	–

The matrix can be read row-wise, from left to right as:

• A _i ≺ A _j = A _i ↑ A _j .

The table can be summarized into a count:

• +1 for the preferred algorithm.

Using this counting method, we can order them from highest to lowest, and such ordering would express the preference for (P₁, I₁). If we wanted to generalize for

P = { P 1 , … , P m }

I = { I 1 , … , I l }

The mean and standard deviation of all the indicators for each algorithm were calculated and are presented in Table 1. By combining these results with the hypothesis tests, we determined the outcomes of the comparisons for each indicator. Next, using the results from all comparisons, we obtained the scores for each algorithm across all indicators (see Table 2). Finally, we determined the preferences for each indicator by tabulating the frequency of occurrences for each permutation (see Table 3).

OPEN IN VIEWER

Table 1.

Indicators by algorithm.

Indicator	ME	MNSLC	QN	NSLC	SO
N a	4.61 (2.21)	2.41 (1.28)	3.00 (1.55)	1.19 (.52)	1.62 (1.04)
D gc	2.76 (.86)	4.92 (1.19)	2.25 (.80)	4.39 (1.05)	1.20 (.49)
C	.37 (.02)	.36 (.02)	.35 (.02)	.36 (.02)	.28 (.02)
F	6.10 (2.81)	3.94 (2.53)	5.96 (3.20)	2.68 (1.21)	5.91 (3.07)
D g	3.63 (.83)	5.84 (1.36)	3.04 (.86)	4.90 (.99)	1.41 (.56)
D m	34.72 (2.71)	51.01 (8.68)	42.56 (15.20)	52.24 (7.70)	4.86 (3.20)
D gm	2.03 (.47)	2.47 (1.61)	1.81 (.56)	1.80 (.41)	.50 (.33)
Q D a n a n	1678.72 (881.86)	493.25 (359.01)	803.88 (547.56)	335.00 (182.46)	269.88 (203.65)
Q D f f	1295.69 (705.72)	240.92 (101.44)	371.71 (217.51)	193.60 (82.90)	125.29 (76.88)
Q D f u n	5929.06 (453.52)	6278.15 (333.42)	5254.00 (334.64)	6233.63 (418.72)	3967.39 (475.09)
N u	15.54 (1.48)	24.58 (3.68)	15.65 (5.25)	25.25 (3.06)	2.68 (1.86)

OPEN IN VIEWER

Table 2.

Scores by algorithm for each indicator.

Indicator	ME	MNSLC	NSLC	QN	SO
N a	4	0	−4	2	−2
D gc	0	4	2	−2	−4
C	4	2	0	−2	−4
F	4	−2	−4	1	1
D g	0	4	2	−2	−4
D m	−2	2	4	0	−4
D gm	2	4	−2	0	−4
Q D a n a n	4	0	−2	2	−4
Q D f f	4	0	−2	2	−4
Q D f u n	0	4	2	−2	−4
N u	−2	2	4	0	−4

OPEN IN VIEWER

Table 3.

Permutation incidence.

1	2	3	4	5	I
ME	MNSLC	NSLC	QN	SO	1
ME	QN	MNSLC	NSLC	SO	2
ME	QN	MNSLC	SO	NSLC	1
ME	QN	SO	MNSLC	NSLC	1
MNSLC	ME	QN	NSLC	SO	1
MNSLC	NSLC	ME	QN	SO	3
NSLC	MNSLC	QN	ME	SO	2

6.3. Choosing the winner

To choose the winner, we applied the Borda count method, which is obtained as follows for each algorithm

C B ( M E ) = 5 ( 5 ) + 4 ( 1 ) + 3 ( 3 ) + 2 ( 2 ) + 1 ( 0 ) = 42

C B ( M N S L C ) = 5 ( 4 ) + 4 ( 3 ) + 3 ( 3 ) + 2 ( 1 ) + 1 ( 0 ) = 43

C B ( N S L C ) = 5 ( 2 ) + 4 ( 3 ) + 3 ( 1 ) + 2 ( 3 ) + 1 ( 2 ) = 33

C B ( Q N ) = 5 ( 0 ) + 4 ( 4 ) + 3 ( 3 ) + 2 ( 4 ) + 1 ( 0 ) = 33

C B ( S O ) = 5 ( 0 ) + 4 ( 0 ) + 3 ( 1 ) + 2 ( 1 ) + 1 ( 9 ) = 14

The Borda count results (see Table 4) indicate that the MNSLC algorithm is the winner, closely followed by ME. The order of the experiments based on this count is as follows: MNSLC, ME, NSLC, QN, SO.

OPEN IN VIEWER

Table 4.

Borda count.

Position	ME	MNSLC	NSLC	QN	SO
1	5	4	2	0	0
2	1	3	3	4	0
3	3	3	1	3	1
4	2	1	3	4	1
5	0	0	2	0	9
Borda count	42	43	33	33	14

Another way to choose the winner is through the Condorcet method, where:

• A matrix of n × n is created.

Using the Condorcet method, we again find that the winner is MNSLC (see Table 5). However, as it is evident, the victory of MNSLC was marginal against ME, under both criteria. In terms of the Borda count, MNSLC scored 43, against the 42 of ME, with ME being in the first place five times, one more than MNSLC, which would have made ME the winner by the majority criterion. The difference was marked by the two more times that MNSLC was in the fourth place. With the Condorcet method, as mentioned before, the victory of MNSLC was also marginal, with six wins out of the 11 contests, compared to five of ME. If we examine Table 1, it becomes clear that the victories of MNSLC were almost exclusively focused on indicators that consider purely diversity (D _m , N _u , D _g , D _gc , D _gm , Q D f u n ), with the exception of coverage ( C ) . This means that MNSLC can find greater diversity, both morphological and control-related, but ME is superior in all criteria that consider global fitness and local quality, finding not only greater fitness (F), but also accumulated quality ( Q D f f ) , aligned accumulated novelty ( Q D a n a n ) , and aligned novelty (N _a ).

OPEN IN VIEWER

Table 5.

Condorcet method.

	ME	MNSLC	NSLC	QN	SO
ME	0	5	6	9	11
MNSLC	6	0	9	7	10
NSLC	5	2	0	6	9
QN	2	4	5	0	10
SO	0	1	2	0	0

Now, let’s go back to the questions raised in Section 3.5.4. With the above analysis, we have answered how ME compares with a Multi-BC QD algorithm. To answer how the introduction of global pressure affects a QD algorithm like MAP-Elites and if it still maintains desirable characteristics of QD algorithms, we need to compare it against NSLC.

Returning to examine the Borda count in Table 4, we can observe that ME significantly outperforms NSLC. However, examining the Condorcet method in Table 5, we realize that ME only marginally outperforms NSLC, with six wins compared to NSLC’s five, winning in almost the same indicators as MNSLC, with the addition of D _gm . This brings us to a similar conclusion as with ME and MNSLC. Considering the desirable characteristics of a QD algorithm, we find that ME was superior to both MNSLC and NSLC in terms of the quality found in diversity, but generally inferior in terms of the notions of diversity considered.

When considering the results of MNSLC versus NSLC, it seems that the victory was quite considerable in both Borda count and Condorcet method. Examining the matches one by one in Table 5, we can see that MNSLC won 9 out of the 11 matches against NSLC, while NSLC won the remaining ones. Furthermore, if we examine the two indicators in which NSLC won, we can observe that the only advantage NSLC has is that it achieves slightly more morphological diversity and novelty (D _m , N _u ), specifically 1.23 and .67 on average, respectively, small but statistically significant differences. In both indicators, NSLC is in first place, followed immediately by MNSLC. Therefore, we can say that, although there is a loss of diversity and novelty relative to NSLC, and both indicators are key to premature convergence of morphology, this loss is minimal, albeit statistically significant. Similarly, taking into account the 9 indicators in which MNSLC was superior, we can observe that not only was MNSLC superior in all those related to fitness performance in the task, but it was also superior in two out of the four related to phenotype diversity and in all those related to genotype diversity. Moreover, the phenotypic diversity indicators in which it performed better were those related to the fitness grid. Considering the observations regarding the comparison of MNSLC and NSLC, it seems that the addition of aligned novelty, firstly, and as already mentioned in the results, is indeed an appropriate notion of aligned novelty, and secondly, effectively mitigates the problem of slow fitness progress in QD algorithms, specifically in NSLC, without compromising performance in terms of morphological diversity.

When comparing the performance of QN with its “purely” QD counterparts (MNSLC and NSLC), we find that MNSLC outperforms QN both in Borda count and in the Condorcet method, while QN and NSLC have equal performance in Borda, and NSLC outperforms QN in Condorcet by one encounter. We can notice that the 4 victories of QN over MNSLC were exclusively in the indicators that consider global or local quality, and similarly, the 7 victories of MNSLC over QN were exclusively in the indicators that consider genotypic and phenotypic diversity. Regarding the 5 victories of QN over NSLC, we have a situation very similar to the previous one, except that a very marginal victory in genetic diversity of morphology (D _gm ) is added to the QN victory count. Therefore, we can say that taking QN as a reference, the use of local competition leads to more diverse results in both genotype and phenotype, while the use of global fitness leads to results with higher global quality and higher quality on the diversity found.

We have already talked about the relative advantages and disadvantages among the different algorithms that present diversity protection, now we will talk about the absolute advantages and disadvantages that algorithms with diversity protection present compared to not presenting any type of diversity protection, using the simple genetic algorithm (SO) as a reference.

As far as the performance of QN compared to SO is concerned, the former outperforms the latter by more than double the number of Borda points, and QN wins all its encounters against SO, except for one, where there was a tie (the null hypothesis could not be rejected), which corresponds to the fitness, that is, the effect of adding morphological novelty as a second objective to the fitness of solutions to the virtual creature morphology and control optimization problem in the walking pattern task, when the first objective is global quality, is null. And while there is no effect on fitness, there are improvements in all other fitness-related indicators ( Q D f f ; Q D a n a n , N _a ), as well as all indicators related to genetic and phenotypic diversity. That is, the addition of morphological novelty not only does not affect global fitness, but it also improves related criteria such as aligned novelty, aligned novelty found in morphological diversity, and quality in morphological diversity. It also improves in all aspects related to diversity, including morphological diversity, and therefore, the robustness against the premature convergence phenomenon of morphology regarding control.

Comparing NSLC against SO, we found again an advantage of more than double the Borda count points, and NSLC wins 9 of its 11 encounters against SO. SO outperforms NSLC in the two population indicators related to fitness (F and N _a ), and loses in any other indicator that incorporates some notion of diversity, from quality indicators in diversity ( Q D f f ; Q D a n a n ) to genotype and phenotype diversity indicators. Therefore, we can affirm that adding morphological novelty and changing the global fitness to local fitness leads to more diverse results in both genotype and phenotype, with higher quality over the discovered diversity, but with worse global fitness.

Moving on to MNSLC versus SO, we now have the best algorithm versus the worst one with respect to the Borda count ordering. This time, MNSLC wins 10 out of 11 encounters, while SO outperforms MNSLC only in one indicator, the global fitness, so it is noteworthy the effect of adding aligned novelty as a third objective, since MNSLC now outperforms SO in that indicator. With these observations, we can affirm that adding novelty in morphology and alignment as second and third objectives respectively, and changing global fitness for local fitness, leads to more diverse genotype and phenotype results, with higher quality over the found diversity, and with higher aligned novelty, but at the expense of worse global fitness.

Finally, comparing ME and SO, we observe that ME ranks second by a narrow margin, with SO coming in last. In this case, ME outperforms SO in each and every indicator, indicating that adding a selection pressure proportional to global fitness to MAP-Elites results in an algorithm that not only finds solutions that are more diverse and of higher quality in terms of diversity but also solutions that are globally fitter.

NSLC exhibits the worst progress in terms of fitness, which is consistent with the findings in the literature (Pugh et al., 2016). However, it is difficult to determine with this indicator alone whether it is only slow progress or a complete stall. On the other hand, MNSLC manages to outperform NSLC, but it appears that, when examining the last generations in the convergence graph, it would not grow much in subsequent generations after 1000. In the other algorithms, we can see that there is still a growing trend in those last generations, with ME presenting the clearest trend. However, even SO shows signs that it will continue to grow beyond generation 1000. Therefore, it is possible that if the experiments had been allowed to run for more generations, the convergence of SO to a final fitness value would have been observed.

As mentioned several times, behavioral novelty is indeed aligned with fitness; however, it is possible to obtain aligned novelty with high fitness through a not-so-novel trajectory, if the creature covers a greater distance than any previous solutions within its own niche. Conversely, aligned novelty can also be low if the creature does not achieve such improvement. Another way to achieve aligned novelty is through novel trajectories, which are not necessarily associated with high fitness values and may be associated with changes in morphology that lead to erratic behaviors. Therefore, the reason why ME and QN outperform MNSLC in aligned novelty may be because both algorithms are capable of finding solutions that are both morphologically novel and have higher global fitness than MNSLC, covering both scenarios in which high aligned novelty can occur, that is, through diverse trajectories and consistent progress in global fitness. One way to examine this would be to explore the angles of the trajectories of each individual and the frequency of times per generation that the individual with maximum fitness was surpassed. This last metric could also be related to progress relative to fitness, and therefore to convergence to the final fitness value. Thus, it could also be tested whether higher aligned novelty is related to greater progress in terms of fitness.

Regarding the main findings that the examination of morphological diversity provides, it can be mentioned that there is undoubtedly a clear effect in the introduction of diversity protection mechanisms, which allows for much more diverse populations (i.e., with higher values of average morphological distance), and that not only is there a slower loss of morphological diversity, but also that morphological diversity is maintained, thus avoiding premature convergence of morphology with respect to the control. Regarding unaligned novelty, something very similar can be mentioned, except that in this case, it is not maintained for any algorithm, but rather, at best decreases slowly. It should be noted that this indicator is closely related to the parameters chosen for the archive, whereas at least the test problems, plots with logarithmic shape (with small values of k) or plots that reach a maximum value only to immediately start decreasing (when the archive had reached its maximum size) could also be found.

The coverage provides information about the ability of each algorithm to explore different morphologies, and in consistence with all indicators that incorporate some type of diversity, SO is found to perform much worse than the other experiments. ME outperforms all other experiments in this regard, despite having a probability of bin selection proportional to fitness, however, a uniform selection probability would likely achieve greater coverage.

Similarly, a loss of genetic diversity was found to varying degrees in all algorithms, with SO showing the greatest loss of genetic diversity, consistent with diversity indicators of the phenotype at the population level. The experiments that showed the least loss of genetic diversity were the QD algorithms. Just as in novelty and morphological diversity, ME had the worst performance of the three QD algorithms. These indicators are population-based, so it makes sense to argue that when introducing a parent selection proportional to global quality, rather than uniform, certain individuals have a greater probability of being selected as parents of the next generation, sharing characteristics that make them fitter, and therefore similar to one another, which in turn results in less diversity and novelty.

It is worth noting that both MNSLC and NSLC, after performing the non-dominated sorting of NSGA-II, choose the survivors of the last front with a ranking based on genetic diversity, instead of crowding distance, which can also play a role not only in genetic diversity, but also in novelty and morphological diversity. If we compare with QN, we can argue something similar about how global selection pressure affects diversity, and the possible role that genetic diversity may play as a cutoff criterion in the last front of NSGA-II. It does not seem too risky to assume that the higher genetic diversity, morphological diversity, and morphological novelty presented by MNSLC and NSLC are due not only to local competition but also to the use of genetic diversity in NSGA-II; however, the necessary data collection was not carried out to determine this conclusively.

The genetic diversity of control presents the same ordering as genetic diversity and seems to cause a greater diversity of behavior, which in turn causes a higher aligned novelty (with MNSLC being the best in this indicator). Furthermore, the two experiments with local competition present greater diversity in the controllers and it could be argued that maximizing both local competition score and aligned novelty increases the genetic diversity of the controllers.

It is interesting to see that NSLC and QN are so similar in terms of morphological genetic diversity. It seems that morphological genetic diversity ignores local competition, and is only interested in morphological phenotypic diversity, and the increase in aligned novelty that it entails.