# Understanding fitness landscapes in morpho-evolution via local optima
networks

Sarah L. Thomson  Edinburgh Napier UniversityEdinburghScotland, UK
<EMAIL_ADDRESS>,  Léni Le Goff  Edinburgh Napier
UniversityEdinburghScotland, UK<EMAIL_ADDRESS>,  Emma Hart  Edinburgh
Napier UniversityEdinburghScotland, UK<EMAIL_ADDRESS>and  Edgar Buchanan
University of YorkYorkEngland, UK<EMAIL_ADDRESS>
(2024)

###### Abstract.

Morpho-Evolution (ME) refers to the simultaneous optimisation of a robot’s
design and controller to maximise performance given a task and environment.
Many genetic encodings have been proposed which are capable of representing
design and control. Previous research has provided empirical comparisons
between encodings in terms of their performance with respect to an objective
function and the diversity of designs that are evaluated, however there has
been no attempt to explain the observed findings. We address this by applying
Local Optima Network (LON) analysis to investigate the structure of the
fitness landscapes induced by three different encodings when evolving a robot
for a locomotion task, shedding new light on the ease by which different
fitness landscapes can be traversed by a search process. This is the first
time LON analysis has been applied in the field of ME despite its popularity
in combinatorial optimisation domains; the findings will facilitate design of
new algorithms or operators that are customised to ME landscapes in the
future.

fitness landscape analysis; evolutionary robotics; local optima networks;
indirect representation

††copyright: acmcopyright††doi: 10.1145/nnnnnnn.nnnnnnn††isbn: 978-x-xxxx-
xxxx-x/YY/MM††conference: The Genetic and Evolutionary Computation Conference
2024; July 15–19, 2024; Melbourne, Australia††journalyear: 2024

##  1\. Introduction

Ever since the pioneering work of Karl Sims almost 30 years ago (Sims, 1994),
evolutionary approaches to simultaneously optimise both the designs (i.e. the
body) and controllers of robots have gathered pace. We term this morpho-
evolution (ME) in this article. A variety of approaches now exist ranging from
those that apply evolution purely in simulation (Lipson et al., 2016), through
hybrid approaches that combine hardware and software evolution (Hale et al.,
2019) to those that evolve purely in hardware (Brodbeck et al., 2015).
Evolutionary frameworks that facilitate ME must consider the challenging
question of how to represent both designs and controllers in a form that can
be evolved. Direct encodings such as directed trees (Veenstra and Glette,
2020) provide a one-to-one mapping between a genotype encoding the design and
control that is realised as the phenotype, and have been postulated to
facilitate the fine-tuning of parameters (Veenstra et al., 2022). On the other
hand, indirect encodings such as Compositional Pattern Producing Networks
(CPPNs) (Stanley, 2007) and L-Systems (Lindenmayer and Jürgensen, 1992) result
in many-to-one genotype-phenotype mappings which have several potential
benefits, e.g. in being able to produce repeating patterns. This
characteristic is generally useful in robotics as it can facilitate the
creation of symmetrical designs, e.g. with equal numbers of actuators on each
side of a body. However, when using an indirect encoding, the genotype and the
phenotype spaces are not necessarily isometric: small changes in the genotype
can lead to major changes in the phenotype, and vice-versa (Hart and Le Goff,
2022; Angus et al., 2023)). Several authors have also noted that the
evolutionary progress using an indirect encoding is slower than that of a
direct one (e.g. (Le Goff et al., 2020)). While there have been several
studies that have compared direct and indirect encodings for jointly
optimising design and control in terms of the level of performance reached
(i.e. objective fitness) and the diversity of designs of evolved robots
(Veenstra et al., 2017; Miras et al., 2018b; De Carlo et al., 2020; Veenstra
and Glette, 2020; Veenstra et al., 2022), to the best of our knowledge there
has been no attempt to explain these empirical observations. That is, there
has been no attempt to study the fitness landscapes induced by a choice of
encoding to shed light on the ease by which these landscapes can be traversed
during a search process. Understanding these landscapes could help explain
empirical observations regarding performance or diversity, but could also
facilitate the design of new algorithms or operators that are customised to
the type of landscapes found.

The novel contribution of this paper is to apply Local Optima Network (LON)
analysis (Ochoa et al., 2014) to analyse the structure of the fitness
landscapes created by different encodings when jointly optimising the body and
control of a simulated robot. LONs are a technique first introduced in 2008
(Ochoa et al., 2008) which enables the information contained in a search space
to be compressed into a single graph. Analysis of the graph enables the
calculation of a set of metrics which characterise the structure of a
landscape. The technique has proved very popular in the combinatorial
optimisation community (Ochoa and Veerapen, 2018; Ochoa and Chicano, 2019;
Muñoz-Herrera and Suchan, 2022), and more recently with continuous functions
(Mitchell et al., 2023) but has not been used in evolutionary robotics. We
pose the question _”To what extent can applying LON analysis to fitness
landscapes induced by different encodings in morpho-evolution provide an
explanation for (1) observed performance and (2) the extent to which the
search-space is explored”._ To answer this, we apply LON analysis to
landscapes induced by three different encodings (a direct encoding producing a
tree, L-System, CPPN) for an ME locomotion task, shedding new light on the
effect of solution encoding on search-space _navigability_. We define
navigability as the ability for the search to a) consistently escape local
optima, and b) explore a large number of genotypes. This is the first time
that LON analysis has been applied within Evolutionary Robotics (ER).
Furthermore, it is also novel in applying LON analysis to representations in
which there is no direct equivalence between genotype and phenotype: there is
only one other attempt at this that we are aware of in which LON analysis is
applied to Linear Genetic Programming (Hu et al., 2023). The results indicate
that L-System indirect encoding leads to better fitness landscape
navigability; this helps to explain its superior performance in terms of
fitness in some empirical studies (Veenstra and Glette, 2020; Veenstra et al.,
2022). In contrast, the CPPN encoding is shown to induce lower navigability —
in fact, it seems to have landscapes which exhibit large numbers of low-
quality local optima which are difficult to escape.

##  2\. Related Work

Recent years have seen a significant increase in efforts focused on the joint
optimisation of both design and control of robots using evolutionary methods
(Auerbach et al., 2014; Brodbeck et al., 2015; Lipson et al., 2016; Miras et
al., 2020; Hale et al., 2019; Li et al., 2023). Much of this work is inspired
by the seminal work of Pfiefer & Bongard (Pfeifer and Bongard, 2006) ‘How the
Body Shapes the Way We Think’ which argues that thought is not independent of
the body but is tightly constrained, and at the same time enabled, by it —
suggesting the evolution of design and cannot control be considered
independently. Several evolutionary frameworks have been proposed to
simultaneously evolve design and control (Miras et al., 2018b; Hale et al.,
2019; Gupta et al., 2021; Miras et al., 2020; Li et al., 2023). Regardless of
the choice of framework, a decision must be made regarding how to encode both
design and control on a single genome. A number of encodings have been
suggested, which can be categorised as direct or indirect. In the former,
typical examples include the use of directed trees (Auerbach et al., 2014;
Veenstra et al., 2017; Veenstra and Glette, 2020) which for example specify
how modules are connected to create a design, as well as the control
parameters of each module. Indirect encodings require a mapping from genotype
to phenotype. The use of CPPN (Stanley, 2007; Cheney et al., 2014; Kriegman et
al., 2020) is common here: a particular advantage of this type of encoding is
that it can produce repeating patterns, a characteristic which is often
beneficial in producing symmetric robots which facilitate locomotion. An
alternative indirect encoding is an L-System (Lindenmayer and Jürgensen, 1992)
in which an evolved grammar produces rules that specify the robot (Miras et
al., 2018a; Miras and Eiben, 2019). A number of studies have examined the
effect of the choice of encoding on performance (i.e. objective fitness) in
joint optimisation of design and control (Veenstra et al., 2017, 2022). In
(Veenstra and Glette, 2020), the effect of the encoding on solution diversity
is also examined. Miras et. al. (Miras, 2021) compare two generative encodings
in terms of how the choice of encoding influences the phenotype of evolved
robots. In (Pigozzi et al., 2023), the extent to which the encoding influences
learning ability in joint optimisation. However, there have been very few
attempts to relate observed empirical performance metrics to the underlying
fitness landscape induced by the choice of encoding, particularly in morpho-
evolution. Some early literature in fitness landscape analysis in ER focused
solely on the controller of the robot (Dittrich et al., 1999). In (Smith et
al., 2002) a simple form of morpho-evolution was investigated in which the
sensor placement of a fixed robot body is evolved along with the controller.
The authors investigate ruggedness and neutrality within this landscape. Naya-
Varela et. al. (Naya-Varela et al., 2023) propose an evolutionary algorithm
for growth-based morphological development in which the robot structure
‘grows’ overtime following a fixed developmental schedule. In this work, the
authors employed a graph-based tool called Search Trajectory Networks (STN)
(Ochoa et al., 2021) which enables the trajectory of the best solutions found
by a population-based algorithm to be visualised over time. The STN reveals
information about how the algorithm traverses the underlying fitness
landscape. However, this algorithm does not fall in the class of ‘morpho-
evolution’ algorithms considered in this article as the design changes
according to a pre-determined schedule, rather than being evolved. Outside of
ER, LONs have been used to bring new insights into the optimisation dynamics
of a wide variety of different problems e.g. in combinatorial optimisation
(Ochoa et al., 2008; Pavelski et al., 2021), neural architecture search (NAS)
(Ochoa and Veerapen, 2022), in algorithm-selection (Bożejko et al., 2018) and
in continuous landscapes (Ochoa et al., 2014). All of the aforementioned
examples apply LON analysis to encodings in which there is a one-to-one
mapping between genotype and phenotype. The only example of which we are aware
in which LON analysis is applied to a scenario in which there is a many-to-one
mapping is a recent analysis of Linear Genetic Programming (Hu et al., 2023).
In summary, this article adds to the field of ER by using LONs for the first
time to understand fitness landscapes, in the context of jointly optimising
design and control using different encodings.

##  3\. Preliminaries

This section provides a brief primer on LONs and defines the terms that are
used throughout the remainder of the article in this respect, given that this
is the first time they have been used in the context of Evolutionary Robotics.
LONs (Ochoa et al., 2008) are a means to study the global structure of a
fitness landscape. In a LON, nodes represent local optima in a fitness
landscape, and edges between nodes are search transitions between them. By
_search transition_ , we mean that two locations in the configuration space
were linked during a metaheuristic search: in our study, a search transition
between nodes is an application of perturbation followed by hill climbing. It
is therefore a sequence of search operations applied to the source node which
results in the destination node. LONs can be created by sampling local optima
in a fitness landscape using Iterated Local Search (ILS) (Lourenço et al.,
2003). We now define some key terminology:

##### Neighbourhood.

The _neighbourhood_ of a solution, $s_{i}$, are the solutions which are
adjacent to $s_{i}$ according to a neighbourhood function: $N(s)$. In this
work, the notion of adjacency depends on the design encoding type: the
particulars of these are detailed in Section 4.4.1.

##### LON nodes.

A local optimum has superior or equal fitness to its neighbours according to a
fitness function _f_. In this work, we do not exhaustively search the
neighbourhood: this would be computationally infeasible. Instead, we consider
that a solution $lo_{i}$ is a local optimum if it has superior or equal
fitness to its _sampled_ neighbourhood $SN$. Formally: $\forall{n\in
SN(lo_{i})}:$ $f(lo_{i})>=f(n)$ (assuming maximisation, as is the case for
this study) where $SN(lo_{i})$ is the sampled neighbourhood, $n$ is a
particular neighbour. The nodes in a LON, $LO$, are the local optima as just
defined.

##### LON edges.

There is an edge from local optimum $lo_{i}$ to local optimum $lo_{j}$, if
$lo_{j}$ can be obtained after applying a random perturbation to $lo_{i}$
followed by local search, and $f(lo_{j})\geq f(lo_{i})$. In LON terminology,
these are called _escape_ edges (Verel et al., 2012). The edges are called
monotonic because they record only non-deteriorating, directed connections
between local optima. Edges are weighted with the frequency of transition: the
number of times during searches that $lo_{j}$ was reached by applying
perturbation then local search to $lo_{i}$. The set of edges is denoted by
$E$.

##### Local optima network (LON).

A local optima network, LON = $(LO,E)$, consists of nodes $lo_{i}\in LO$ which
are the local optima, and edges $e_{ij}\in E$ between pairs of nodes $lo_{i}$
and $lo_{j}$ with weight $w_{ij}$ iff $w_{ij}>0$. A LON which only includes
neutral or improving transitions between local optima is a _monotonic_ LON, or
MLON; we construct these for the present work. For the interested reader,
detailed descriptions of MLONs can be found in previous literature (Ochoa et
al., 2017).

##  4\. Methodology

###  4.1. Task and Environment

We conduct LON analysis on a joint optimisation task taken from Veenstra _et
al._ (Veenstra et al., 2022) which focused on empirical performance
comparisons of robots evolved using different encodings. The goal is to evolve
a robot that needs to traverse a virtual landscape consisting of a set of
gentle hills. The objective is to maximise the distance travelled. We use an
environment called gym2D which facilitates robot design evolution for OpenAI
Gym (Brockman et al., 2016); its foundation was the bipedal walker environment
111https://github.com/openai/gym/blob/master/gym/envs/box2d/bipedal_walker.py.
The beginning section of the environment we use is depicted in Figure 1. A
virtual robot begins at the left of a rectangular box containing a virtual 2D
landscape and the objective is to move as far to the right of the box as
possible; the horizontal axis for the box is associated with a scale between 0
and 100 (0 is no movement; 100 is reaching the end) and the fitness of a robot
is how far along this axis it travels. There is also a ”kill-switch”
implemented to reduce the computational cost: if a robot travels slower than a
minimum speed (set at 0.04) then the attempt terminates and the robot is
ascribed a default fitness of 5.0. The same environment configuration is used
across all the experiments. A robot can be seen attempting to travel across
the landscape in Figure 1.

Figure 1. An example robot attempting to move horizontally through the virtual
environment. Each shape is a module, and each module has its own controller.

###  4.2. Robots

A robot is defined by a set of modules combined to produce a 2D shape. Modules
are either circular or rectangular and each has parameters defining its size.
The allowable size ranges for modules are as follows: a width and height
between 0.5 and 1, and (in the case of circle shapes) a radius between 0.25
and 0.5. Modules also each have a controller, as described shortly. A list of
eight modules is kept in memory and these are selected to be used as nodes to
build the robot. For each individual run, a new initial module list is
generated of four circles and four rectangles. The radius of the circles is
randomly generated between the bounds just mentioned; similarly, the width and
height of the rectangles are randomised. The parameters of the modules are
mutated during the evolutionary runs. Each virtual robot has one controller
per module. A controller is encoded as a sine wave with parameters: amplitude
($\alpha$, frequency $\theta$, phase $\delta$, and offset $\varepsilon$,

(1) $y(t)=\alpha sin(\theta t+\delta)+\varepsilon$

The ranges for the parameters are: $\alpha$ and $\delta$ can vary between -1
and +1; $\theta$ between -0.1 and +0.1; and $\varepsilon$ between $-\pi$ to
$+\pi$. These parameters are also subject to mutation (details of mutation are
found later in Section 4.4.1. Finally, the maximum depth of the tree which
represents a robot is seven, while the maximum size (number of allowable
modules) is 40. These parameters are taken directly from the literature
(Veenstra et al., 2022).

###  4.3. Encodings

The phenotype of every virtual creature in Gym2D is represented as a tree (see
the left sub-plot of Figure 2). We evaluate three different encodings to
produce this tree. The first encoding is direct, therefore exactly represents
the tree. Two encodings are indirect therefore there is a mapping from
genotype to phenotype. The same three encodings were used in (Veenstra and
Glette, 2020). For the indirect encoding, the tree starts from an axiom and
then expand iteratively until reaching a maximum depth. Each node can have a
maximum of three nodes attached. The following section describes how each
encoding generates a tree.

Figure 2. On the left: directed tree created by an L-System encoding after two
iterations. The nodes labelled $C_{1}$ and $R_{1}$ show a circle and rectangle
module, respectively. On the right: an example of a virtual 2D robot; its
graph is visualised in the top right. Figure is from (Veenstra et al., 2022)

##### Direct encoding.

The direct robot design encoding is a directed tree structure which represents
the topology of the robot’s design. For each node that makes up a robot, the
encoding describes: the node’s index in the tree; the index of the module from
the module list it represents; the index of its parent node; and the parent
connection coordinates.

Figure 3. Diagram representing how an L-system works

##### L-System.

For the L-System encoding, a parametric approach is used where symbols and
grammar rules are used (see figure 3). A set of symbols is specified via an
alphabet. The alphabet is composed of eight symbols, each of which corresponds
to a module in the modules list. A set of rules, known as _grammar rewrite
rules_ , are applied to the symbols in order to construct the tree. A rule is
composed of an input symbol and substitution symbols, traditionally, the input
symbols are replaced by the substitution symbols. In this system, starting
from an axiom symbol, the rules are used to add new leaves (each leaf
corresponds to one of the substitution symbols) to the tree from a previously
created node (the input symbol).

Figure 4. Diagram representing a CPPN network

##### CPPN

CPPNs are neural networks specifically designed to produce patterns where the
inputs are spatial coordinates. Each network node is a function representing a
basic pattern such as Gaussian for axial symmetry, sinusoidal for repetition,
or sigmoid for central symmetry. As CPPN is similar to a feed-forward neural
network it is classically trained using neuro-evolution. In this study, we
make use of the Neuro-evolution through Augmenting Topologies (NEAT) (Stanley
and Miikkulainen, 2002). To create the tree, the CPPN is queried incrementally
for each connection site, i.e. where a node can be potentially attached. The
CPPN takes three inputs: the depth of the connection site in the tree, the
parent index in the module list and the angle of connection. Thus, the CPPN
produces 6 outputs: whether a module should be attached or not, the type of
module, i.e. the index in the modules list, and four outputs corresponding to
the controller parameters (see figure 4).

###  4.4. LON Construction

Following the convention in the LON literature (see, for example: (Treimun-
Costa et al., 2020; Veerapen et al., 2017; Thomson et al., 2023; Pavelski et
al., 2021)), we sample local optima by executing runs of Iterated Local Search
(ILS). ILS (Lourenço et al., 2003) is a single-point meta-heuristic which
conducts repeated cycles of applying large perturbations followed by local
search. Although algorithms with recombination are typical in evolutionary
robotics and neuro-evolution, here we use mutation only — following the
majority of work in LONs — but return to the question of whether crossover
should be used later. At each iteration during the LON construction process, a
local optimum is obtained; an acceptance condition then dictates whether it is
accepted as the new incumbent solution, or whether the previous local optimum
should be obtained. In this way, a process reminiscent of a Markov chain is
followed. We gather data from 30 runs of ILS, each commenced from random
starting solutions during the LON construction process. ILS combines local
search with random perturbations. The local search employs a first-improvement
pivot rule. This means that whenever an improving (or equal) neighbour is
found, it immediately becomes the incumbent genotype (Tari and Ochoa, 2021);
this is in contrast to best improvement, where the best possible neighbour is
found. This would require too much computation here. In the experiments, the
local search process stops when there has been 100 operations applied with no
improvement. Each of the runs terminates after 30 iterations with no
improvement in local optimum quality or after 100 iterations in total. The
local search conducts random mutation on both the robot controller and the
robot design according to pre-defined probabilities described in the next
section. During the runs, every accepted local optimum is hashed and logged as
a network node, and any transition between local optima is logged as network
edges. When a local optimum is added, its fitness, an index, and its hashed
phenotype are noted too, in addition to the hashed genotype. Edges are
recorded as the combination of the source node index and the destination node
index, and weight is logged as well. Edge weights are the number of times a
transition was followed. The LON is an amalgamation of the nodes and edges
from the separate runs.

####  4.4.1. Mutation operators

We use the mutation settings which were tuned in previous literature (Veenstra
and Glette, 2020); these are provided in Table 1. Perturbation is the mutation
operation applied three times. Only improving or equal local optima are
accepted during the search. As mentioned, both the controllers and the designs
are mutated and each has a separate application rate (see Table 1). For the
design, the approach to mutation depends on the encoding type:

##### Direct encoding.

Mutation on the direct robot design encoding has three possible operations:
removing a node, adding a node to the tree structure, and modifying a module’s
node. For each possible location for adding a node, a node is added according
to the mutation rate; for removals, the mutation rate is halved (owing to the
fact that sub-nodes are removed alongside the original node; that is, it is a
more dramatic mutation). Finally, the module’s nodes are modified by mutating
the shape (width, height, radius) and connection angle using Gaussian
mutation.

##### L-System.

Mutation on the L-System robot design encoding takes place in the following
way: the rules are changed by adding or removing substitution symbols and the
shape (width, height, radius) and angle of the modules in the list are mutated
using Gaussian mutation.

##### CPPN

For the CPPN encoding, the mutation is the default NEAT mutation from the
neat-python222https://neat-python.readthedocs.io/en/latest/ library. This has
several operations modifying the network structure: adding a connection;
deleting a connection; adding a node; deleting a node; and replacing bias. All
of these are carried out according to the robot design mutation rate. There
are a large number of parameters associated with the CPPN encoding. We use the
values provided in a configuration
file333https://github.com/FrankVeenstra/gym_rem2D/blob/master/ModularER_2D/NeuralNetwork/config
and associated with previous literature (Veenstra et al., 2022). As well as
the CPPN mutation, the shape (width, height, radius) and angle of the modules
in the list are mutated according to the mutation probability.

##### Controller mutation.

For all three encodings, the same controller representation is used (recall
Equation 1). For the direct and L-System encodings, controllers’ parameters
are mutated according to a mutation rate (controller rate in Table 1) using
Gaussian mutation. The standard deviation of the Gaussian distribution is set
at 0.2 (the same as in (Veenstra and Glette, 2020)). For the CPPN, the process
is different: in that case, each of the four parameters which define a
controller is an output from the network (as described in Section 4.3 and with
Figure 4); they are therefore affected when the NEAT mutation is carried out
on the CPPN.

Table 1. Mutation rates in the LON construction algorithm

parameter | direct | L-System | CPPN
---|---|---|---
rate (controller) | 0.32 | 0.16 | 0.02
rate (design) | 0.16 | 0.04 | 0.02

###  4.5. Experimental Setup

The majority of our experimental setup matches that of Veenstra _et. al_ , who
compared morphological encodings (Veenstra and Glette, 2020; Veenstra et al.,
2022). We begin from the associated online
repository444https://github.com/FrankVeenstra/gym_rem2D as a base for our
implementation; from now on, we will refer to this as _veenstra-repository_.
All parameters associated with the encodings and environment are from there.
For each of the three encodings, we conducted 30 independent ILS runs to log
local optima and their edges; these were amalgamated to form a single LON
object for each encoding.

Statistics are calculated using the Mann-Whitney u-test (Mann and Whitney,
1947) to reject the hypothesis that two distributions are statistically the
same where * means $p<0.05$, ** means $p<0.005$, *** means $p<0.0005$ and ****
means $p<0.00005$ and $p$ represent the probability.

Analysis of the LONs includes considering their weakly connected components.
In directed networks, a _weakly connected component_ (WCC) is a group of nodes
which are all reachable from each other, even if edges are uni-directional or
in different directions; WCCs are identified in the analysis using Tarjan’s
algorithm (Tarjan, 1974) and are referred to as _components_ from now on.

##  5\. Results

###  5.1. Visualisation

Figures 5(a)-5(d) display visualisations of the extracted LONs; there is a
separate Figure for each of the three robot design encoding types. In the
plots, larger size and darker colour of nodes indicate better fitness. Node
size is proportional to fitness. For colour, we took the approach of
separating nodes into three fitness levels according to the _overall
distribution of fitness_ across all three LONs: nodes with fitness in the
first quartile (below 9.59) are low-fitness and are coloured very pale purple.
Nodes within the interquartile range are middle-fitness and are light purple
in colour. Finally, nodes with fitness in the upper quartile — which is 24.48
or above — are high-fitness and are dark purple. Self-loops are represented by
edges which curve out of and back into the left of a node. The force-directed
network layout algorithm _kamada kawai_ (Kamada et al., 1989) in R has been
used to layout the LONs. Let us first consider the direct encoding LON shown
in Figures 5(a) (with self-loops) and 5(b) (without self-loops). It is
immediately evident from the latter that this is a network containing several
isolated components; in fact, there are 30, equal to the number of ILS runs
used to construct the LON. This is not a typical phenomenon seen in other
domains where LON analysis has been used, e.g. (Treimun-Costa et al., 2020;
Mitchell et al., 2023; Cleghorn and Ochoa, 2021) which usually result in more
densely connected graphs. This is likely due to the vastness of the search
space. In evolutionary robotics, the search space can be essentially infinite
(most controllers are encoded as real-valued numbers) and the goal is not
necessarily to find a global optimum. Instead, algorithms are typically
applied to find a _good enough_ robot that functions well.

(a) Direct encoding with self-loops

(b) Direct encoding without self-loops

(c) L-System encoding

(d) CPPN encoding

Figure 5. LONs for the three robot design encodings. larger size and darker
colour of nodes indicate better fitness. Node size is proportional to fitness.
For colour: nodes with fitness in the first quartile of the overall fitness
distribution across the three LONs (this is a fitness of less than 9.59) are
low-fitness and are coloured very pale purple. Nodes within the interquartile
range are middle-fitness and are light purple in colour. Finally, nodes with
fitness in the upper quartile — which is 24.48 or above — are high-fitness and
are dark purple. Self-loops are represented with edges which curve out of and
back into the right of a node.

In Figure 5(b) we can see that the 30 components are chains — each node has
either zero or one outgoing directed edges, and either zero or one incoming
directed edges. We note from Figure 5(a) that the direct encoding LON exhibits
several self-loops: these are recorded when an escape attempt from a local
optimum is made, but returns to the same local optimum. Self-loops are not
present for the LONs of the other two encodings, indicating that it may
sometimes be difficult to escape local optima when using the direct encoding
(i.e., there are strong attractors present). Additionally, the chains of the
direct encoding (Figure 5(b)) seem to be shorter when compared to the other
two encodings. This indicates a lack of _navigability_ associated with this
encoding: it may struggle to explore genotypes (indeed, we can see from Table
2 that this encoding resulted in far fewer unique designs being considered
than the L-System). Despite this, we notice from Figure 5(b) that the direct
encoding chains often terminate with a good-fitness (dark purple) local
optimum. Looking across Figures 5(a)-5(d), the divergence in fitness
distributions can be seen. The CPPN LON in Figure 5(d) contains a lot of nodes
which are very pale purple, indicating their fitness is in the first quartile
of the three-LON fitness distribution. Comparing that observation with Figure
5(c), where the LON for the L-System encoding is shown, we note that there are
more chains which contain the dark purple; that is, high-fitness (in the upper
quartile) nodes are reached in more searches. The chains for this LON appear
to be longer than those evident in the direct encoding LON (Figure 5(a)), and
they very often end up with fitness in the third quartile (above 24.48). These
observations, when taken together, indicate the _navigability_ associated with
an L-System encoding.

Another observation from the CPPN LON in Figure 5(d) is that some of the
chains are lengthy but do not display any notable improvement in fitness. This
tells us that the fitness landscape associated with using this encoding is not
easily navigable: there appear to be large poor-quality (where poor-quality is
having a fitness in the first quartile: less than 9.59) local optima plateaus
which are difficult to escape (under the studied search operations and
parameters). Looking at the direct and L-System encodings in Figures 5(a) and
5(c), and remark that they appear to be associated with better fitness
landscape navigability than the CPPN: the chains in their LONs more often
display incremental improvements in fitness, which can be seen in the increase
in size of node throughout a chain. This is particularly evident with the
L-System LON in Figure 5(c). It seems that the search becomes trapped at low-
quality local optima with these encoding less often. This helps to explain the
higher performance which has been observed using L-System design encoding in
the literature (Veenstra et al., 2022).

###  5.2. Features

In this Section we will first consider statistics about the ILS runs conducted
with the different encodings. Table 2 shows the rate of mutation acceptance;
the rate of design acceptance (this is computed as the number of accepted
mutations which were associated with a changed design); the number of unique
designs present; and the total number of attempted mutations. From Table 2, we
notice that the mutation acceptance is lowest for direct encoding and the CPPN
encoding has the lowest design acceptance. The highest rates are associated
with the L-System encoding. The other two encodings tried out substantially
lower numbers of unique designs than the L-System, despite having a larger
number of attempted mutations. These facts, taken together, help to ratify the
observations in the previous section that the direct encoding struggles to
widely explore and find fitness improvement, while the L-System encoding
facilitates wide exploration and moving towards better fitness.

Table 2. ILS run statistics for the three encodings

metric | direct | L-System | CPPN
---|---|---|---
_mutation acceptance_ | 10.52% | 47.65% | 33.77%
_design acceptance_ | 32.15% | 38.10% | 19.38%
_unique designs_ | $165\,800$ | $227\,411$ | $147\,443$
_attempted mutations_ | $352\,438$ | $297\,964$ | $328\,970$

Table 3 reports measurements for the LONs associated with each encoding type.
We introduce now the metrics whose names or meanings may not be intuitive:
_components_ is the number of weakly connected components in the network. The
metric _path length_ is the average edge distance between any two nodes in the
network (excluding nodes which are unreachable to each other); _degree_ is the
average degree of a node, including self-loops; and _infeasible_ is the
percentage of nodes which had fitness indicating that the robot activated the
killswitch and can therefore be considered infeasible.

Table 3. LON statistics

metric | direct | L-System | CPPN
---|---|---|---
_nodes_ | 269 | 360 | 546
_edges_ | 248 | 330 | 516
_components_ | 30 | 30 | 30
_path length_ | 4.37 | 6.32 | 12.66
_degree_ | 1.84 | 1.83 | 1.89
_infeasible_ | 0% | 0.003% | 34.43%

The plots in Figure 6 show distributions for run and LON metrics across the
three three encoding types. For each, there is a boxplot with a horizontal
line indicating the median. We also overlay as a swarm the actual data points.
For each pair of distributions, the u-test result is shown. Figure 6(a)
displays the distribution of all fitnesses in the sampled LON nodes (local
optima) of each of the three robot design encoding choices: direct, L-System,
and CPPN.

(a) All fitnesses in the LON

(b) Maximum fitness per run

(c) Number of evaluations

(d) Fitness delta

(e) Chain length

Figure 6. Distribution for various run and LON metrics across the three three
encoding types. If there is a significant difference (according to a u-test)
between two distributions, this is indicated with a horizontal bar and
$p$-value level: *$p<0.05$; **$p<0.005$; ***$p<0.0005$; ****$p<0.00005$

We can notice from visual inspection that the three distributions have very
different medians, and the swarms show that the overall shape of the
distributions are divergent to one another as well. In terms of IQR and
median, the L-System has the best fitness range, with the direct encoding
being next-best. This was also shown in a previous study (Veenstra et al.,
2022). The CPPN encoding is especially distinct from the other two; when
including outliers, it has a much wider range, and includes more lower-quality
fitnesses than those of the other two. Another interesting phenomenon is that
in the case of L-System and CPPN, there appear to be groups of local optima
with the same fitness: at approximately 25 for the L-System and direct
encodings, and approximately 10 for the CPPN. Despite having the lowest IQR (a
finding which ratifies those in the literature (Veenstra and Glette, 2020)),
the CPPN encoding resulted in the highest fitness sampled. Figure 6(b) shows
distributions for the _maximum_ fitness found in each run. We can see from
visual inspection that the CPPN is associated with the widest range. Indeed,
it has both the highest value and the lowest. The encoding which shows the
highest (best) IQR and median is the L-System type. The plot in Figure 6(c)
shows, for each of the three encodings, the distribution for the number of
fitness evaluations used by the 30 runs. Notice from comparing the boxes
visually that the distributions are similar for all three encodings.

In Figure 6(d) we present the data concerning fitness deltas. A delta is here
defined as the fitness increase associated with a directed LON edge (recall
that the edges are always directed towards equal or better fitness, so the
delta is always zero or above). For the plot, we compute the _median_ delta
within each component in the LON; these correspond to the 30 runs. Those
medians are the data shown in the Figure. Observing the three encodings, we
note that the direct encoding is associated with the largest fitness deltas.
This can be seen with its higher IQR and median line. Additionally, the
largest value in across all encodings is achieved by the direct encoding. The
L-System distribution has the lowest and narrowest IQR, implying that the
fitness deltas are the smallest (but perhaps most consistent in magnitude) of
the three. In terms of significant difference, the L-System and CPPN pair are
statistically the same; however, the direct/L-System and direct-CPPN pairs
differ, with $p<0.005$ and $p<0.05$ respectively. Finally, the distributions
of chain lengths in the LON (there is one for each of the 30 runs) are
captured in Figure 6(e). A chain is the path of local optima followed by one
run of the ILS. We observe from the plot that the direct encodings had the
shortest chains; this was followed by the L-System, and then CPPN with the
longest chains. The direct-L-system pair has an associated $p<0.05$, the
L-system-CPPN pair are $p<0.05$ different and the direct-CPPN are
statistically different with $p<0.0005$.

###  5.3. The Elephants in the Room

One consideration of our approach is that we did not use a crossover operator.
The reader may wonder whether the ILS used here actually explores promising
regions of the search space, or whether instead all the sampled local optima
are poor when compared to those that might be found by crossover-based
approaches. To this end, we note that the fitness levels reached as shown in
Figure 6(b) are similar to those shown in previous literature — where an
evolutionary algorithm with crossover was used (Veenstra et al., 2022) with
the same virtual environment and both direct and L-System encoding. It is
therefore reasonable to assume that the ILS approach reaches fitness levels
approximately equivalent to those achieved by a crossover-based algorithm, and
that at least some of the local optima sampled are of acceptable quality.
Another factor to consider when interpreting the results is that we used
mutation parameters from the literature ((Veenstra and Glette, 2020)). Those
were tuned for an algorithm which also included crossover, so it is probable
that they are not the optimal settings when applied in an iterated local
search. This being said, being consistent with the parameters from previous
work is a reasonable place to start when beginning a new line of research. In
this study, we looked at different possible genotype encoding to explore one
phenotype space: a tree-based structure. Therefore, the present findings are
limited to this phenotype space. In future work, it will be interesting to
look at how different encodings explore different design representations such
as voxel-based.

##  6\. Conclusions

We have conducted a study with the aim of providing insight into performance
differences associated with different robot design encodings which have been
observed in the literature. To this end, we carried out fitness landscape
analysis with local optima networks (LONs). Three design encodings for a
virtual 2D robot were considered: a direct tree structure representation, an
L-System, and a type of neural network (a CPPN). We sampled the LONs
associated with each of these, visualising them and comparing their features.
The results showed that an L-System encoding is associated with good fitness
landscape navigability: search has the ability to escape local optima and
discover a large number of unique designs. In contrast, although being able to
produce a large number of unique designs, the CPPN encoding appears to be
linked with low-quality local optima from which it is hard to escape towards
better fitness. This observation indicates that CPPN encoding is not suited
for creating the structure used to define a robot morphology in this study: a
tree-based structure. We also found that a direct encoding can lead to search
struggling to leave local optima; however, good fitnesses can sometimes be
found when it does. These findings provide insight into previously-observed
performance differences between encodings in the literature. We therefore
propose that LONs can be a valuable tool for explaining phenomena seen in
observed in evolutionary robotics search can inform the result of better
algorithms and operators in future. The code and data associated with this
work will be published in an online Zenodo repository upon acceptance.

## References

  * (1)
  * Angus et al. (2023) Mike Angus, Edgar Buchanan, Léni K Le Goff, Emma Hart, Agoston E Eiben, Matteo De Carlo, Alan F Winfield, Matthew F Hale, Robert Woolley, Jon Timmis, et al. 2023\.  Practical hardware for evolvable robots.  _Frontiers in Robotics and AI_ 10 (2023).
  * Auerbach et al. (2014) Joshua Auerbach, Deniz Aydin, Andrea Maesani, Przemyslaw Kornatowski, Titus Cieslewski, Grégoire Heitz, Pradeep Fernando, Ilya Loshchilov, Ludovic Daler, and Dario Floreano. 2014.  Robogen: Robot generation through artificial evolution. In _Artificial Life Conference Proceedings_. MIT Press One Rogers Street, Cambridge, MA 02142-1209, USA journals-info …, 136–137.
  * Bożejko et al. (2018) Wojciech Bożejko, Andrzej Gnatowski, Teodor Niżyński, Michael Affenzeller, and Andreas Beham. 2018.  Local optima networks in solving algorithm selection problem for tsp. In _International Conference on Dependability and Complex Systems_. Springer, 83–93.
  * Brockman et al. (2016) Greg Brockman, Vicki Cheung, Ludwig Pettersson, Jonas Schneider, John Schulman, Jie Tang, and Wojciech Zaremba. 2016.  Openai gym.  _arXiv preprint arXiv:1606.01540_ (2016).
  * Brodbeck et al. (2015) Luzius Brodbeck, Simon Hauser, and Fumiya Iida. 2015.  Morphological evolution of physical robots through model-free phenotype development.  _PloS one_ 10, 6 (2015), e0128444.
  * Cheney et al. (2014) Nick Cheney, Robert MacCurdy, Jeff Clune, and Hod Lipson. 2014.  Unshackling evolution: evolving soft robots with multiple materials and a powerful generative encoding.  _ACM SIGEVOlution_ 7, 1 (2014), 11–23.
  * Cleghorn and Ochoa (2021) Christopher W Cleghorn and Gabriela Ochoa. 2021.  Understanding parameter spaces using local optima networks: a case study on particle swarm optimization. In _Proceedings of the Genetic and Evolutionary Computation Conference Companion_. 1657–1664.
  * De Carlo et al. (2020) Matteo De Carlo, Eliseo Ferrante, and AE Eiben. 2020.  Comparing indirect encodings by evolutionary attractor analysis in the trait space of modular robots. In _Proceedings of the 2020 Genetic and Evolutionary Computation Conference Companion_. 73–74.
  * Dittrich et al. (1999) Peter Dittrich, Andre Skusa, Wolfgang Banzhaf, and Wolfgang Kantschik. 1999.  Dynamical properties of the fitness landscape of a GP Controlled random morphology robot. In _Proceedings of the 1st Annual Conference on Genetic and Evolutionary Computation-Volume 2_. Citeseer, 1002–1008.
  * Gupta et al. (2021) Agrim Gupta, Silvio Savarese, Surya Ganguli, and Li Fei-Fei. 2021.  Embodied intelligence via learning and evolution.  _Nature communications_ 12, 1 (2021), 5721.
  * Hale et al. (2019) Matthew F Hale, Edgar Buchanan, Alan F Winfield, Jon Timmis, Emma Hart, Agoston E Eiben, Mike Angus, Frank Veenstra, Wei Li, Robert Woolley, et al. 2019\.  The are robot fabricator: How to (re) produce robots that can evolve in the real world. In _Artificial Life Conference Proceedings_. MIT Press One Rogers Street, Cambridge, MA 02142-1209, USA journals-info …, 95–102.
  * Hart and Le Goff (2022) Emma Hart and Léni K Le Goff. 2022.  Artificial evolution of robot bodies and control: on the interaction between evolution, learning and culture.  _Philosophical Transactions of the Royal Society B_ 377, 1843 (2022), 20210117.
  * Hu et al. (2023) Ting Hu, Gabriela Ochoa, and Wolfgang Banzhaf. 2023.  Phenotype Search Trajectory Networks for Linear Genetic Programming. In _European Conference on Genetic Programming (Part of EvoStar)_. Springer, 52–67.
  * Kamada et al. (1989) Tomihisa Kamada, Satoru Kawai, et al. 1989\.  An algorithm for drawing general undirected graphs.  _Information processing letters_ 31, 1 (1989), 7–15.
  * Kriegman et al. (2020) Sam Kriegman, Douglas Blackiston, Michael Levin, and Josh Bongard. 2020.  A scalable pipeline for designing reconfigurable organisms.  _Proceedings of the National Academy of Sciences_ 117, 4 (2020), 1853–1859.
  * Le Goff et al. (2020) Léni K Le Goff, Emma Hart, Alexandre Coninx, and Stéphane Doncieux. 2020.  On pros and cons of evolving topologies with novelty search. In _Artificial Life Conference Proceedings 32_. MIT Press One Rogers Street, Cambridge, MA 02142-1209, USA journals-info …, 423–431.
  * Li et al. (2023) Wei Li, Edgar Buchanan, Léni K Le Goff, Emma Hart, Matthew F Hale, Matteo De Carlo, Robert Wooley, Alan F Winfield, Jonathan Timmis, and Andy M Tyrrell. 2023.  Evaluation of frameworks that combine evolution and learning to design robots in complex morphological spaces.  _IEEE Transactions on Evolutionary Computation_ (2023).
  * Lindenmayer and Jürgensen (1992) Aristid Lindenmayer and Hans Jürgensen. 1992.  Grammars of development: discrete-state models for growth, differentiation, and gene expression in modular organisms.  In _Lindenmayer systems: Impacts on theoretical computer science, computer graphics, and developmental biology_. Springer, 3–21.
  * Lipson et al. (2016) Hod Lipson, Vytas Sunspiral, Josh Bongard, and Nicholas Cheney. 2016.  On the difficulty of co-optimizing morphology and control in evolved virtual creatures. In _Artificial life conference proceedings_. MIT Press One Rogers Street, Cambridge, MA 02142-1209, USA journals-info …, 226–233.
  * Lourenço et al. (2003) Helena R Lourenço, Olivier C Martin, and Thomas Stützle. 2003.  Iterated local search.  In _Handbook of metaheuristics_. Springer, 320–353.
  * Mann and Whitney (1947) Henry B Mann and Donald R Whitney. 1947.  On a test of whether one of two random variables is stochastically larger than the other.  _The annals of mathematical statistics_ (1947), 50–60.
  * Miras (2021) Karine Miras. 2021.  Constrained by design: Influence of genetic encodings on evolved traits of robots.  _Frontiers in Robotics and AI_ 8 (2021), 672379.
  * Miras and Eiben (2019) Karine Miras and Agoston E Eiben. 2019.  Effects of environmental conditions on evolved robot morphologies and behavior. In _Proceedings of the Genetic and Evolutionary Computation Conference_. 125–132.
  * Miras et al. (2020) Karine Miras, Eliseo Ferrante, and Agoston E Eiben. 2020.  Environmental influences on evolvable robots.  _PloS one_ 15, 5 (2020), e0233848.
  * Miras et al. (2018a) Karine Miras, Evert Haasdijk, Kyrre Glette, and AE Eiben. 2018a.  Effects of selection preferences on evolved robot morphologies and behaviors. In _Artificial Life Conference Proceedings_. MIT Press One Rogers Street, Cambridge, MA 02142-1209, USA journals-info …, 224–231.
  * Miras et al. (2018b) Karine Miras, Evert Haasdijk, Kyrre Glette, and Agoston E Eiben. 2018b.  Search space analysis of evolvable robot morphologies. In _Applications of Evolutionary Computation: 21st International Conference, EvoApplications 2018, Parma, Italy, April 4-6, 2018, Proceedings 21_. Springer, 703–718.
  * Mitchell et al. (2023) Paul Mitchell, Gabriela Ochoa, and Romain Chassagne. 2023.  Local optima networks of the black box optimisation benchmark functions. In _Proceedings of the Companion Conference on Genetic and Evolutionary Computation_. 2072–2080.
  * Muñoz-Herrera and Suchan (2022) Sebastián Muñoz-Herrera and Karol Suchan. 2022.  Local Optima Network Analysis of Multi-Attribute Vehicle Routing Problems.  _Mathematics_ 10, 24 (2022), 4644.
  * Naya-Varela et al. (2023) Martín Naya-Varela, Andrés Faíña, and Richard J Duro. 2023.  Guiding the Exploration of the Solution Space in Walking Robots Through Growth-Based Morphological Development. In _Proceedings of the Genetic and Evolutionary Computation Conference_. 1230–1238.
  * Ochoa and Chicano (2019) Gabriela Ochoa and Francisco Chicano. 2019.  Local optima network analysis for MAX-SAT. In _Proceedings of the Genetic and Evolutionary Computation Conference Companion_. 1430–1437.
  * Ochoa et al. (2021) Gabriela Ochoa, Katherine M Malan, and Christian Blum. 2021.  Search trajectory networks: A tool for analysing and visualising the behaviour of metaheuristics.  _Applied Soft Computing_ 109 (2021), 107492.
  * Ochoa et al. (2008) Gabriela Ochoa, Marco Tomassini, Sebástien Vérel, and Christian Darabos. 2008.  A study of NK landscapes’ basins and local optima networks. In _Proceedings of the 10th annual conference on Genetic and evolutionary computation_. 555–562.
  * Ochoa and Veerapen (2018) Gabriela Ochoa and Nadarajen Veerapen. 2018.  Mapping the global structure of TSP fitness landscapes.  _Journal of Heuristics_ 24, 3 (2018), 265–294.
  * Ochoa and Veerapen (2022) Gabriela Ochoa and Nadarajen Veerapen. 2022.  Neural Architecture Search: A Visual Analysis. In _International Conference on Parallel Problem Solving from Nature_. Springer, 603–615.
  * Ochoa et al. (2017) Gabriela Ochoa, Nadarajen Veerapen, Fabio Daolio, and Marco Tomassini. 2017.  Understanding phase transitions with local optima networks: number partitioning as a case study. In _Evolutionary Computation in Combinatorial Optimization: 17th European Conference, EvoCOP 2017, Amsterdam, The Netherlands, April 19-21, 2017, Proceedings 17_. Springer, 233–248.
  * Ochoa et al. (2014) Gabriela Ochoa, Sébastien Verel, Fabio Daolio, and Marco Tomassini. 2014.  Local optima networks: A new model of combinatorial fitness landscapes.  _Recent advances in the theory and application of fitness landscapes_ (2014), 233–262.
  * Pavelski et al. (2021) Lucas Marcondes Pavelski, Marie-Éléonore Kessaci, and Myriam Delgado. 2021.  Local Optima Network Sampling for Permutation Flowshop. In _2021 IEEE Congress on Evolutionary Computation (CEC)_. IEEE, 1131–1138.
  * Pfeifer and Bongard (2006) Rolf Pfeifer and Josh Bongard. 2006.  _How the body shapes the way we think: a new view of intelligence_.  MIT press.
  * Pigozzi et al. (2023) Federico Pigozzi, Federico Julian Camerota Verdù, and Eric Medvet. 2023.  How the Morphology Encoding Influences the Learning Ability in Body-Brain Co-Optimization. In _Proceedings of the Genetic and Evolutionary Computation Conference_. 1045–1054.
  * Sims (1994) Karl Sims. 1994.  Evolving 3D morphology and behavior by competition.  _Artificial life_ 1, 4 (1994), 353–372.
  * Smith et al. (2002) Tom Smith, Andrew Philippides, Phil Husbands, and Michael O’Shea. 2002.  Neutrality and ruggedness in robot landscapes. In _Proceedings of the 2002 Congress on Evolutionary Computation. CEC’02 (Cat. No. 02TH8600)_ , Vol. 2. IEEE, 1348–1353.
  * Stanley (2007) Kenneth O Stanley. 2007.  Compositional pattern producing networks: A novel abstraction of development.  _Genetic programming and evolvable machines_ 8 (2007), 131–162.
  * Stanley and Miikkulainen (2002) Kenneth O Stanley and Risto Miikkulainen. 2002.  Evolving neural networks through augmenting topologies.  _Evolutionary computation_ 10, 2 (2002), 99–127.
  * Tari and Ochoa (2021) Sara Tari and Gabriela Ochoa. 2021.  Local search pivoting rules and the landscape global structure. In _Proceedings of the Genetic and Evolutionary Computation Conference_. 278–286.
  * Tarjan (1974) Robert Endre Tarjan. 1974.  A new algorithm for finding weak components.  _Inform. Process. Lett._ 3, 1 (1974), 13–15.
  * Thomson et al. (2023) Sarah L Thomson, Nadarajen Veerapen, Gabriela Ochoa, and Daan van den Berg. 2023.  Randomness in local optima network sampling. In _Proceedings of the Companion Conference on Genetic and Evolutionary Computation_. 2099–2107.
  * Treimun-Costa et al. (2020) German Treimun-Costa, Elizabeth Montero, Gabriela Ochoa, and Nicolás Rojas-Morales. 2020.  Modelling parameter configuration spaces with local optima networks. In _Proceedings of the 2020 Genetic and Evolutionary Computation Conference_. 751–759.
  * Veenstra et al. (2017) Frank Veenstra, Andres Faina, Sebastian Risi, and Kasper Stoy. 2017.  Evolution and morphogenesis of simulated modular robots: a comparison between a direct and generative encoding. In _Applications of Evolutionary Computation: 20th European Conference, EvoApplications 2017, Amsterdam, The Netherlands, April 19-21, 2017, Proceedings, Part I 20_. Springer, 870–885.
  * Veenstra and Glette (2020) Frank Veenstra and Kyrre Glette. 2020.  How different encodings affect performance and diversification when evolving the morphology and control of 2D virtual creatures. In _Artificial Life Conference Proceedings 32_. MIT Press One Rogers Street, Cambridge, MA 02142-1209, USA journals-info …, 592–601.
  * Veenstra et al. (2022) Frank Veenstra, Martin Herring Olsen, and Kyrre Glette. 2022.  Effects of encodings and quality-diversity on evolving 2D virtual creatures. In _Proceedings of the Genetic and Evolutionary Computation Conference Companion_. 164–167.
  * Veerapen et al. (2017) Nadarajen Veerapen, Fabio Daolio, and Gabriela Ochoa. 2017.  Modelling genetic improvement landscapes with local optima networks. In _Proceedings of the Genetic and Evolutionary Computation Conference Companion_. 1543–1548.
  * Verel et al. (2012) Sébastien Verel, Fabio Daolio, Gabriela Ochoa, and Marco Tomassini. 2012.  Local optima networks with escape edges. In _Artificial Evolution: 10th International Conference, Evolution Artificielle, EA 2011, Angers, France, October 24-26, 2011, Revised Selected Papers 10_. Springer, 49–60.