## Abstract

Several previous studies have focused on modelling and analysing the collective dynamic behaviour of population-based algorithms. However, an empirical approach for identifying and characterising such a behaviour is surprisingly lacking. In this paper, we present a new model to capture this collective behaviour, and to extract and quantify features associated with it. The proposed model studies the topological distribution of an algorithm's activity from both a genotypic and a phenotypic perspective, and represents population dynamics using multiple levels of abstraction. The model can have different instantiations. Here it has been implemented using a modified version of self-organising maps. These are used to represent and track the population motion in the fitness landscape as the algorithm operates on solving a problem. Based on this model, we developed a set of features that characterise the population's collective dynamic behaviour. By analysing them and revealing their dependency on fitness distributions, we were then able to define an indicator of the exploitation behaviour of an algorithm. This is an entropy-based measure that assesses the dependency on fitness distributions of different features of population dynamics. To test the proposed measures, evolutionary algorithms with different crossover operators, selection pressure levels and population handling techniques have been examined, which lead populations to exhibit a wide range of exploitation-exploration behaviours.

## 1 Introduction

Population-based algorithms explore a search space by sampling it to gather information and then redistributing the population within it based on such information. This biases the movement of the population in a way that is hoped to guide the algorithm to examine known good areas of the search space with more intensity and/or to try to locate new promising areas. Naturally, it is well known that whether an algorithm's bias is beneficial or not depends on the match between that bias and the fitness measure used to inform the search (Wolpert and Macready, 1997). Nonetheless, as the search mechanisms move, create and/or eliminate individuals, the population shows a collective behaviour which often changes dynamically over time.

Understanding and characterising the collective dynamic behaviour of population-based algorithms are thus very important. For example, identifying emergent collective properties related to an algorithm's dynamics allows the comparison of algorithms in terms of behaviours, not just performance (e.g., the number of fitness evaluations to reach an optimum). Therefore, this is not only a prerequisite for designing new population-based algorithms for specific applications, but it may also enable practitioners to make technically sound choices of algorithms, operators, parameters, and so forth.

### 1.1 Previous Empirical and Theoretical Studies

Various previous approaches have been tried to develop a better understanding of the collective behaviour of population-based algorithms. An approach is to practically study the dynamic behaviour of an algorithm with respect to a specific problem by defining and analysing suitable empirical measures (Bethke, 1981; Goldberg, 1989; Mitchell et al., 1992; Jones, 1995; Tomassini et al., 2005; Vanneschi, 2007; Graff and Poli, 2008; Poli and Vanneschi, 2007; Vanneschi et al., 2004). Typically, efforts have concentrated on trying to characterise landscape features and landscape/search-algorithm interactions with the objective of understanding what makes a problem easy or hard for evolutionary algorithms rather than on the characterisation of the behaviour of populations.

An alternative approach is to theoretically model the dynamics of search algorithms. Within evolutionary computation theory, approaches have included Markov chain formulations (Vose and Liepins, 1991; Davis and Principe, 1993; Poli et al., 2004; Mitavskiy and Rowe, 2006), schema theory (Holland, 1992; Stephens and Waelbroeck, 1999; Poli et al., 2004; Poli and McPhee, 2003a, 2003b; Stephens and Poli, 2007), computational complexity techniques (Droste et al., 2002; Jansen et al., 2005; Witt, 2006; Jansen and Wegener, 2002; Storch and Wegener, 2004; Neumann and Wegener, 2007), the statistical-mechanics formulation (Prügel-Bennett and Shapiro, 1994; Shapiro et al., 1994), and the characterisation of the fundamental limitations of search (Wolpert and Macready, 1997; Whitley and Watson, 2005; Schumacher et al., 2001; Igel and Toussaint, 2003; Poli et al., 2009). All of these have seen some success at mathematically modelling evolutionary algorithms.

The main questions that the aforementioned theoretical studies have tried to address are whether an algorithm is able to reach an optimal solution (i.e., convergence), how fast this solution can be reached (i.e., convergence speed), and/or modelling the changes of a population's state over time (i.e., population dynamics). Here we will focus on areas of theory that investigate population dynamics, namely schema theory, the Markov chain model, and the statistical mechanical model, to position our work in relation to them.

Schema theory captures the dynamics of a population at a somewhat coarse-grained level by describing how the population redistributes itself among a set of schemata, where the schema is a fixed subset of the search space. The model describes how the number of individuals in a schema changes over time depending on the number of individuals in a related set of schemata and their fitness. Markov chain models of EAs use a different, more microscopic representation for the state of the population and describe how the probability distribution over population states changes over time. This model can make precise predictions about population dynamics, but these require considerable mathematical and computational resources. An alternative model, statistical mechanics, has been proposed to represent the dynamics of large complex systems in terms of their statistical properties. Statistical mechanics models population dynamics by predicting the changes of the population's fitness distribution over time.

Note that the first two models focus on a genotypic representation of the population while the third model uses phenotypes (fitness) to represent population dynamics, but no model uses both. Finally, we should note that in all cases representation schemes are static, in the sense that once one has decided what form of schema, state space representation, or fitness-distribution representation to adopt, this remains constant throughout a run. As we will see later, our work departs from these traditions and proposes a dynamic genotype-phenotype representation of population dynamics.

### 1.2 Exploration and Exploitation

While the models presented above are of significant theoretical importance, they have little practical use in characterising realistic algorithms running on realistic problems and extracting features of the collective dynamic behaviour of populations.

Looking at what has been done in this area, we find that researchers often talk about two emergent behaviours of search algorithms: *exploration* and *exploitation*. Intuitively, exploration refers to a behaviour associated with the discovery of new good regions in the search space, while exploitation refers to the exploration previously discovered good regions (Holland, 1992; Eiben and Schippers, 1998; Macready and Wolpert, 1998). However, we should note that there is no precise definition in the literature of the notions of exploration and exploitation and no precise characterisation of the distinction between them.

Clearly, whatever definition one embraces, when discussing exploration and exploitation one often talks about both the search space and the fitness values. After all the genetic operators create new solutions by combining/altering existing (useful) information contained in the genotypic representation of individuals, while selection favours fitter solutions, thereby using the phenotypic value (fitness) of individuals. Possibly because of this and of the lack of a more precise definition, often people simply attribute the exploration features of an evolutionary algorithm entirely to the properties of the genetic operators used by the algorithm (e.g., crossover or mutation), while it is assumed that exploitation is achieved by selection.

This is an unsatisfactory interpretation as the reality appears to be much more complex: each search operator, including selection, has a certain bias in guiding the search, and the aggregation of these biases results in an emergent collective behaviour that exhibits a certain degree of exploration/exploitation at different times of the search.

### 1.3 Contributions of this Paper

In this paper, we develop an empirical coarse-grained model that captures and tracks a population's dynamics represented as topological distribution of various features describing populations from both a genotypic and a phenotypic perspective.^{1} The model uses a variable representation that represents the population's dynamics at a relatively high level of abstraction, which makes it easier to talk about collective behaviours of populations.

Our model represents the population using a grid of nodes, where each node points to a region of population activity and records information about genotypic density, fitness values, and distances for the individuals in such a region. During a run of an algorithm, as the regions of population activity change over time, the nodes in the model adjust to track these regions and continue gathering information about them. Thus, over a run, the model produces a sequence of grids that capture the population dynamics represented as a topological distribution of genotypic density, fitness values, and distances. This model shares some features with previous theoretical models of EAs: it is similar to a schema-based model in that it uses sets of genotypes as elements of the representation, it is similar to Markov chain models in that it tracks the changes of the population, and it is similar to statistical mechanics models as it keeps track of the fitness distribution (but it does so at a finer level of granularity, i.e., at the level of activity nodes).

Higher levels of abstraction in the description of the population's collective dynamic behaviour can be obtained by grouping the nodes in our model into activity regions and by categorising these regions according to the type of activity carried out within them. As we will see later, this allows us to define a variety of measures that can be used to characterise different aspects of the collective behaviour of the population and to analyse its coarse-grained dynamics.

In addition, by exploiting our coarse-grained representation of population dynamics, we contribute to clarifying the notions of exploration and exploitation. In particular, we do so by quantifying by how much the fitness distribution guides the search activities of an algorithm. Although other proposed measures can be qualitatively used to characterise the exploitation/exploration behaviour, we suggest a specific new entropy-based exploration indicator of the dynamic behaviour of population. This assesses the fitness-dependency of the topological distributions of different features of the collective dynamic behaviour of an algorithm. This is particularly easy to do with our model, as it is designed to track the topological distribution of fitness values and other genotypic features of populations.

Many implementations of the proposed coarse-grained representation of populations are possible. In this paper, a modified version of self-organising maps (SOMs; Kohonen, 2001) are used to track population movement in the search space and to mine information about the emergent collective behaviour of an algorithm as it operates on solving a problem.^{2} The model has been applied on different EAs that use different operators and techniques. The reported results show that the model was able to characterise their collective behaviours and reveal the effect of different biases the algorithms used.

The rest of this paper is organised as follows. Section 2 will present the proposed model in detail and will develop measures of the collective dynamic behaviour of populations. It will also aim at assessing the fitness-dependency of different aspects of population activities. In Section 3, we present the details of our implementation of the proposed model using SOMs. Section 4 reports and discusses the results of our experiments. Finally, we end the paper with some conclusions in Section 5.

## 2 Proposed Model of Collective Dynamic Behaviour

In this section, the formulation of the proposed model is presented. Then, we define features of the collective behaviour of an algorithm based on the model. Finally, the proposed exploitation indicator is defined.

### 2.1 Principles of the Proposed Model

The proposed model uses a grid of nodes to capture the topological distribution of individuals in a population-based algorithm at certain points in time.^{3} Each node is a vector representing the centre of a group of individuals in a localised area in the search space. In addition, each node also stores information describing the activities of the group in that area. To track the dynamic behaviour of a population over a run, the grid is updated to represent the motion of the population and record information about its activities as the algorithm works on solving a problem. Emergent features that characterise collective dynamic behaviour are extracted by analysing the changes in node positions and other information stored in the nodes.

Formally, the grid, *C ^{t}*, represents a snapshot of the population distribution at time

*t*. Let

*C*

^{0}be a grid of nodes that represent the distribution of the initial population , where is the population size, represents the position of individual

*i*in the search space, where

*D*is the search space dimension, and is the associated fitness value. Furthermore, let the time-ordered set represent a sequence of all the individuals created by an algorithm over a period of a run from the initial population, where

*N*is the total number of calls to the fitness function. The grid is updated every individuals produced by the algorithm (i.e., represents the sampling period). For example,

*C*

^{1}is computed starting from

*C*

^{0}and using the sequence of individuals . That is, the grid's node vectors are adjusted to track individuals of the sequence

*P*

_{1}and record information about them. Generally speaking, the grid

*C*is created by tracking and recording information about the set of individuals , where . More specifically, if is the subset of individuals associated with a node,

^{t}*C*, based on their distance from the node vector (as explained in Section 3), we have that the sets and representing the positions and fitness values of the individuals associated with

^{t}_{r}*C*provide information about the population activity happening in the region represented by that node.

^{t}_{r}*C*consists of a matrix of nodes, where each node,

^{t}*C*, is represented by the following tuple where

^{t}_{r}*t*is time, is the position in the grid, and the elements of the tuple are as follows:

is a vector representing the centre of a region of activity of the population. As the algorithm works on redistributing the population around activity regions, the node vector tracks the newly created individuals in its vicinity (the mechanism of updating vector positions is explained in Section 3).

(hit distance) is the sum of the distances between and the individuals, , associated with the node. This can be used as an indicator of how far away from the centre of activity an algorithm searches for new solutions.

is the mean fitness of the individuals associated with

*C*. This element gives an insight about the quality of the region of the search space represented by the node.^{t}_{r}is the best fitness value of an individual associated with

*C*since time^{t}_{r}*t*=0. Formally, if , then the best fitness for time*t*>0 is defined as . This quantity tells us how well the population is searching this area of the search space, and whether this activity has led to achieving any growth in fitness.(hits counter) is the number of individuals associated with the node. This number represents the volume of activity an algorithm allocates for a particular area of the search space.

As the algorithm starts exploring the fitness landscape and redistributing the population, some nodes (that represent identified areas of activity at a certain point in time) may be abandoned as the population moves toward different areas. At any point in time, we are only interested in the nodes that have observed some population activity. These nodes are called active nodes and they are defined as . Active nodes in this model represent spots in the search space where population activities have taken place. Nodes are, in a way, similar to schemata, in that both represent subsets of the search space. However, differently from schemata, our nodes have the ability to track the population dynamics over time, and they store information about these dynamic activities. This feature enables us to capture the distributions of activity density, observed fitness values, and other features of the topological distribution of a population.

A higher level of abstraction in describing population activity is obtained by aggregating active nodes into activity regions. An activity region is a set of adjacent active nodes the vectors of which point toward similar directions. Formalising this definition requires the introduction of some further notions.

*Nr*(

*C*), that returns a set of all immediate neighbours of a node that is a member of the set of active nodes : where is a distance threshold, being the largest distance between any two points in the search space. The function above returns the empty set in case

^{t}_{r}*C*has no neighbours. With this in hand, we can define a function , that returns a set of consecutive neighbours between

^{t}_{r}*C*and :

^{t}_{r}*C*, that is: Figure 1 illustrates the proposed model and how it relates to the individuals in the population. The figure also shows how regions of active exploration in the search space are represented by activity regions within the model.

^{t}_{r}The sets of all the nodes belonging to growth regions and non-growth regions are defined as and , respectively. Later, we will also use the notion of best region, .

Note that the different levels of coarse-grained representation proposed above preserve the topological distribution of activity density, fitness values, and other features. In the next two sections, we will define a set of features to represent the dynamic collective and emergent behaviour of a population based on these notions.

### 2.2 Extracting Features of Collective Dynamic Behaviour

Many features that characterise aspects of the collective dynamic behaviour of an algorithm can be extracted by monitoring changes in the activity regions of the model over time and analysing the corresponding recorded information. These features can be classified into three types.

*Activity-related features* describe the way an algorithm allocates its activities to the regions of population. The volume of activity can simply be measured using the number of individuals (hits counter ) that have been associated with nodes of a certain region. Region size (number of nodes in a region) is also used to describe aspects of population activity. These features characterise the population collective behaviour from a higher level (phenotype level). In evolutionary algorithms, selection methods and population handling techniques (e.g., niching techniques) have significant influence on these emergent features and also affect other search operators.

*Distance-related features* describe two aspects of the population collective behaviour from the genotype point of view: how widely the population distributes its activities away from the centre (node vector), and the bias in population movement. Two distance measures can be calculated: the sum of the distances between the node's vector and input individuals () and the displacement, which signifies the amount of shift in population centres of activity and it is calculated as the change of vector positions (the distance between and ).

*Correlated features* analyse the association between two different emergent features. Revealing and quantifying dependencies between emergent features could help in characterising what directs the search and to what extent. In this work, we analyse the dependency of the two types of features mentioned above on fitness values. The resultant measures are used as indicators to the exploitation behaviour of the algorithm.

The first two types of features are introduced in the next section, while the third type is presented in Section 2.3.

#### 2.2.1 Activity-Related and Distance-Related Features

Below we will define four features that characterise the dynamic collective behaviour of the population over time as an algorithm operates on the problem. For each feature we also define a related coarse-grained quantity, the average of that feature over the course of a run, to represent the overall behaviour of an algorithm across time.

*S*is the number of hits these nodes receive collectively at time

^{t}*t*. It measures the intensity with which an algorithm explores these nodes and the regions they represent. It is defined as: where

*S*is the set of all nodes in a set of regions and is the sampling period. This function can be applied to any set of nodes, but it makes most sense when applied to sets of nodes representing activity regions, particularly

^{t}*G*,

^{t}*NG*and

^{t}*Best*. In these cases, this measure tells us the way an algorithm tends to allocate its activities. For example, allocating more activities to growth regions implies that the algorithm exhibits exploitation behaviour as it tends to invest more individuals in regions where good fitness values have been discovered. The average of activity rate over a run of an algorithm is defined as , where is a time-ordered sequence of sets of nodes representing a certain type of region over the period of a run, represents the total number of sampling points, and

^{t}*N*is the total number of fitness function evaluations in a run.

While measures the portion of search activity in a certain type of region, tells us how concentrated these activities were. When , that implies that the activities were evenly distributed among the nodes of *s ^{t}* and each node has been hit by only one individual. On the other hand, if , the activities were more concentrated and at least one node in

*S*has been hit by more than one individual. The smaller , the more concentrated the activities. The average size ratio over a run is defined as , where is a time-ordered sequence of sets of regions of a certain type.

^{t}*M*is defined above and

*AN*is the set of all active nodes at time

^{t}*t*. High displacement values imply that the algorithm tends to move its areas of activity around the search space, while smaller values imply that the algorithm is searching in the same areas.

### 2.3 Using Fitness-Dependency as an Indicator to an Exploitation Behaviour

Population activities can be guided (or biased) by various factors that an algorithm uses to make decisions. In addition to the observed fitness values, for example, distance or similarity between individuals can be used (as in fitness-sharing techniques). In this section, we develop a measure to quantify the level of fitness-dependency in directing the search, and use it as an indicator to population exploitation behaviour.

*X*and

*Y*, and is the entropy of a random variable

*X*.

In information theory, entropy is a measure of the uncertainty or randomness associated with a random variable (Shannon, 1948), while mutual information is a quantity that measures the mutual dependence of the two random variables. In Equation (6), the mutual information can be expressed as , where is the conditional entropy of *Y* given *X*. This quantifies the uncertainty about *Y* left after knowing *X*. So, in effect, the dependency between *X* and *Y* represents how much uncertainty about *Y* has been removed by knowing *X*. The value of lies between zero and one, where 0 means that *X* and *Y* have no association, while 1 means that knowing *X* can totally predict *Y*.

*m*subintervals (

*m*=10 in this work), and we then divide the active nodes into classes according to their feature values. If for are the classes of the active nodes obtained by binning them based on feature , then the probability of each class is computed as

With these definitions in hand, we can use the uncertainty coefficient (Equation (6) to analyse the dependency of hit distance, displacement, and activity (hit count) on fitness. For simplicity, we will denote the corresponding values of the uncertainty coefficient as , respectively. These three measures, in particular, are used as an exploitation behaviour indicator at a certain point in time. The bigger the value they have, the more dependent on fitness the algorithm behaviour is at that time.

The average of fitness-dependency of activity over one algorithm run is defined as . Similarly, we define and . The average fitness-dependency rate is the quantity we will use as an exploitation indicator characterising the collective behaviour of an algorithm over a run.

So far, in this section, we have presented a way to assess the exploitation characteristics of an algorithm in terms of the extent to which the algorithm exploits the sampled fitness values to direct the search. However, developing an accurate measure that quantifies the exploitation behaviour is complicated. This is because algorithms may make use of information other than fitness to redistribute their activities; also, a precise definition of exploitation/exploration is lacking. The same can be said about assessing exploration behaviour. However, if we assume that any action made by an algorithm can be viewed as either exploitation or exploration depending on whether or not information influenced the decision leading to that action, then the proposed fitness-dependency measure quantifies a significant part of population exploitation behaviour.

## 3 Implementing the Proposed Model Using a Self-Organising Map

A SOM (Kohonen, 2001) is an artificial neural network that can be trained using unsupervised learning. After training, SOMs can be used for mapping (classifying) or visualising high dimensional data. SOMs consist of a set of nodes (or neurons) arranged, usually, in a two-dimensional grid. A node, *i*, has an associated vector, *m _{i}*, of the same dimension of the input space and is connected to its neighbours according to a neighbourhood radius, as explained later in this section. The training is done by feeding a SOM with a large number of training samples drawn from the data space. Each time a sample,

*x*, is fed into a SOM, the best matching node (BMN) is identified as the node whose vector has the smallest distance from the input sample, that is, . Then

*m*

_{BMN}is updated by moving it slightly in the direction of

*x*. The change to the BMN's vector results in changing the vectors of its neighbours as well (more on this below).

In this work, the dynamics of the grid nodes *C ^{t}_{r}* of the proposed model is controlled by a SOM learning algorithm. In other words, a SOM is responsible for adjusting the position of each node vector, .

*P*. Node vectors, , are then updated as follows: where

_{init}*x*is the input individual vector, is the learning rate, is the neighbourhood radius, and

*r*is the index of the BMN. During this stage, the initial learning rate is set to and the neighbourhood radius is set to .

_{BMN}^{5}In this phase, these parameters are decreased over time as standard in SOMs. For training iteration

*c*, the learning rate is computed as , and the neighbourhood radius is computed as , where

*I*is the total number of training iterations. Note that, although an individual from the initial population may be used several times in SOM training, its information (as described by the proposed model, Equation (1)) is recorded by the grid only once. The resulting grid,

*C*

^{0}, is topologically consistent and provides a representation of the initial population distribution in the search space.

In the second stage, the initial model is progressively modified (again using a SOM training algorithm) so as to track the collective dynamic behaviour of a population. In this stage, newly created individuals are fed into the SOM to adjust the node vectors. Here, we use a fixed nonzero learning rate and neighbourhood radius so that the nodes can track the population. In addition to the change that new individuals bring to the node vectors of the SOM, information of those individuals (as required in the model) is also recorded in the nodes. This process produces a sequence of grids *C ^{t}*, for , that represent population distributions and activities over time. Collectively, they capture the collective dynamic behaviour of the population over the period of a run.

## 4 Experimental Results

In this section, we first describe a set of algorithms used in the experimentation as well as the operators and techniques used within them and then we experimentally characterise their collective dynamic behaviour using the model and features proposed above.

### 4.1 Algorithms and Experiments Description

To test the proposed model and features, a variety of operators, selection pressure levels, and population handling techniques were utilised to produce drastically different collective behaviours. More specifically, we used both standard evolutionary algorithms^{6} as well as EAs that utilise fitness sharing, deterministic crowding, and random immigrant (details are in Appendix A). With each algorithm, arithmetic crossover, blending crossover (BLX-), and heuristic crossover were used (their details and settings are explained in Appendix B). All EAs used tournament selection.

We chose arithmetic crossover, BLX-, and heuristic crossover as they are representative of operators that exchange genetic material in a stochastic, deterministic, and fitness-biased manner, respectively. These crossover methods show different behaviours (in particular different levels of exploration and exploitation) because they handle the diversity of the population and they make use of the information available to them in different ways (Herrera et al., 2003). However, their behaviour is also a function of their parameter settings.

Uniform arithmetic crossover (Michalewicz, 1996) handles the genetic material in a deterministic fashion. However, for or , this operator exhibits explorative behaviour and increases the population diversity (Herrera et al., 2003). BLX- (Eshelman and Schaffer, 1992) creates new individuals by generating random values for their genes within ranges calculated from the parents’ gene values and . BLX- has a more stochastic behaviour and its exploration/exploitation behaviour can be tuned by adjusting (Nomura and Shimohara, 2001). Heuristic crossover (Herrera and Lozano, 1996; Wright, 1990), on the other hand, creates offspring close to the best parent, which makes this operator more exploitative because it does not only make use of the provided genetic material, but also uses the fitness values to bias the search.

Unfortunately, we cannot assess the effect of recombination operators on the collective dynamic behaviour based only on the way they combine/alter genetic materials to produce new individuals. That is because the quality of the operators’ output depends on the quality of the genetic material which is chosen by the selection mechanism. For this reason, we used tournament selection, where the selection pressure can easily be controlled by increasing/decreasing the tournament size. Larger tournaments induce high selection pressure and lead to high exploitation behaviour; whereas small tournaments lead to low pressure and high exploration.

Our experiments focus on characterising the collective dynamic behaviour of EAs that use different parameter settings for their crossover methods over different levels of selection pressure. The interaction between the biases of different operators creates different collective dynamic behaviours. While tournament selection uses fitness values to bias the search, other population handling techniques involve other elements. Fitness sharing and crowding use distances in handling the population, while the technique of random immigrants tries to counterbalance the bias of selection by adding more stochastic behaviour (randomness) into the search.

We adopted a naming system to refer to the algorithms and their settings when presenting results. An algorithm name consists of four parts representing the population handling technique, the crossover method, the crossover control parameter, and the tournament size. For example, StBLX-0.5-3 is a standard EA using BLX crossover with control parameter value 0.5 and tournament size 3. Full details of the algorithm naming system along with parameter settings, algorithm techniques and benchmark problems are provided in Appendix A.

In showing the results of the experiments, we used two types of figures. The first describes the change of a certain feature over time (fitness function evaluations; e.g., Figure 3, discussed later in this paper), while the second characterises the effect of using different selection pressures (different tournament sizes) on a certain feature of the collective behaviour (e.g., Figure 5, discussed later in this paper). In most of the figures of the second type, we compared the effect of selection pressure on different settings of each crossover (the first row of plots), and we compared the crossover methods with each other using different settings (the lower row of plots). We called these parameter settings balanced, exploitative, and explorative, depending on the qualitative behaviour they induced. We picked the most explorative settings (or the least exploitative in case of Heuristic crossover) from the range of used settings for the explorative parameter settings (BLX-0.75, Heur-0.4, and Arth-1.5). In a similar way, we chose the exploitative parameter settings (BLX-0.0, Heur-0.1, and Arth-0.2). The set of balanced parameter settings represents more moderate and common settings for each of the crossover methods (BLX-0.5, Heur-0.2, and Arth-0.4).

### 4.2 Characterising Emergent Features of Collective Dynamic Behaviour

#### 4.2.1 Illustration of the Operations of the Model

Figure 2 depicts the region formation in the model for four of the algorithms we will characterise later in this section: three standard EAs with different types of crossover and one EA with crowding and arithmetic crossover. The figure shows seven snapshots of the grid of nodes in the model, during the first 1,000 fitness evaluations of a run of each of the four algorithms. All runs had the same initial random seed (and, hence, the same initial population and grid of nodes).

As one can see, only a small subset of the model nodes are used at any given time. Note how the model of the crowding EA holds many more nongrowth nodes than the standard EAs and how these are distributed more evenly within the grid, indicating that correspondingly different regions of the search space are concurrently explored. Note also how the most exploitative of the four algorithms (StHeur-0.2-4) rapidly focuses its search on fewer and smaller (and, hence, more focused) growth regions.

#### 4.2.2 Activity Rate and Size Ratio

To characterise the effect of using different techniques on how the population distributes its focus among different regions of activity, we measured the activity rate (Equation (2)) in growth, nongrowth, and best regions, in a standard EA (StArth-0.4-4) and in EAs with fitness sharing (FSArth-0.4-4), deterministic crowding (CrArth-0.4) and random immigrants (RIArth-0.4-4). All the algorithms use tournament of size 4 (except for deterministic crowding which uses its own selection mechanism) and arithmetic crossover (Arth-0.4). Figure 3 shows the results.

By comparing Figures 3(a) and 3(b) one immediately sees that algorithms using niching techniques allocate little activity to the growth regions, while most of the population focuses on searching nongrowth regions. This is a sign of the explorative nature of this type of algorithm. We also see that both the standard evolutionary algorithm and the random immigrants version have different collective behaviour, as they initially focus on the growth regions more than on the nongrowth ones. Allocating more activities to pursue fitness growth is an indicator of exploitation behaviour. Naturally, the activity rate of these EAs in growth regions declines as the run proceeds as a result of the population's convergence and its losing the ability to achieve fitness growth.

Figure 3(c) depicts activity rate in the best region found by each algorithm. As we can see from the figure, after an initial increase, the activity rate in the best region declines over the course of a run. This phenomenon may contradict the standard interpretation of convergence, which is normally viewed as the population gathering in one region around the best-found individual (especially in standard EAs). At the beginning of a run, a standard EA forms a large activity region around the best-discovered area. As the run goes on, this region shrinks. At some point, the bias of the genetic operators (e.g., BLX-0.75 or Arth-1.5) starts balancing the focusing bias of selection. Then the operators scatter individuals away from the best region, thereby creating small temporary regions in its vicinity. More results and discussion on this issue will be provided later in this section.

In Figure 4, we show the effect of using different crossover methods on the activity rates in different regions of the population in a standard EA. A comparison of Figures 4(a) and 4(b) shows that a standard EA with Heur-0.2 crossover allocates more activities to the growth regions, which is, of course, due to the exploitative nature of this operator. Standard EAs with Arth-0.4 have the lowest interest in growth regions, while an EA with BLX-0.5 changes its character as the population loses diversity and converges. Interestingly, as shown in Figure 4(c), the three algorithms have approximately the same activity rate in the best region.

To analyse the combined effect of using different crossover methods and different selection pressures (tournament sizes) on a population's allocation of activities to regions, we used standard EAs utilising three crossover methods, with four parameter settings for each one, over different selection pressure levels. The parameters were chosen so that a range of exploitation and exploration (or less exploitation, in the case of Heur-) behaviours would be produced. Figure 5 depicts the results.

Let us first focus on Figures 5(a)–(c), which report the average activity rate in growth regions using BLX-, Heur-, and Arth- crossovers. When using exploitative parameter settings (BLX-0.0, BLX-0.25, Arth-0.2, Arth-0.4, Heur-0.1, …, Heur-0.4), the average activity rates in growth regions increases with selection pressure, until a certain point beyond which increasing the selection pressure further has little effect on the activity rate. This happens because highly exploitative operators force the population to converge too quickly onto one region in the search space, thereby losing the ability to discover new good regions and produce further fitness growth. Thus, stepping up selection pressure accelerates losing the ability of producing growth in exploitative crossover methods. The faster the population converges, the fewer growth regions can be found and, thus, the less activity in these regions. Because of this, with exploitative crossover and high selection pressure, the population can only rely on mutation to bring about some diversity that may help discover new regions and move the search away from local optima.

Conversely, as shown in Figures 5(a)–(c), increasing the selection pressure in EAs that use crossover with explorative parameter settings (BLX-0.5, BLX-0.75, Arth-1.2, and Arth-1.5) drives the population to discover more growth regions.^{7} With low selection pressure, explorative crossover has a bias that prevents the population from converging and concentrating its activities only on the good regions that it has so far discovered. Instead, crossover will push individuals away from these regions, thereby reducing the chance of finding better solutions in these regions and, thus, producing growth. However, also in explorative crossover, increasing the selection pressure makes the population converge, although less rapidly, and counterbalances their diversification effects. In this case, the population converges to a small region around the best solutions and the explorative crossover creates some small growth regions by scattering individuals around the best region.

To better illustrate these observations, in Figures 5(d)–(e) we show nine configurations from the previous plots, but this time, we have grouped them on the basis of how exploitative/explorative the crossover's parameters settings are. More specifically, Figure 5(d) shows the activity-rate averages for crossovers using balanced parameter settings, while Figures 5(e) and 5(f) report those for exploitative and explorative settings, respectively.

Search operators and selection methods have a major influence on the size of activity regions explored by an algorithm. For example, using a high selection pressure will shrink the population activity regions more than using a low selection pressure. Figures 6(a)–(g) illustrate this using the size ratio and the activity rate of the best region in an EA using Heur-0.2 crossover. Figure 6 shows that as the selection pressure increases, the activity rate becomes more concentrated. However, the effect of the selection pressure on the size ratio becomes weaker as the selection pressure increases. This is also shown in Figure 6(h), where size ratio and activity rate averages are plotted as a function of the tournament size.

Figures 7(a)–(c) show the effects of the selection pressure on the activity concentration in the best region in standard EAs with different types of crossover. Figure 7(d) compares the size ratios of these algorithms and shows that they are more or less the same, which implies that selection, more than crossover, has a significant effect on changing the size ratio of the best regions of the population.

#### 4.2.3 Hit Distance and Displacement

Hit distance and displacement provide information about the population dynamics from a geometric perspective.

Hit distance measures how far a population searches away from previously sampled regions in the fitness landscape. For example, a large hit distance implies that search operators produce/maintain genetic diversity. The hit distance varies as the run progresses and also depends on the diversity of the population and the nature of the operators used. Figure 8 shows how using different crossover methods or population-handling techniques affects the hit distance.

In particular, Figures 8(a)–(c) illustrate the effect of using different crossover methods. (We have used the same selection pressure level for all results, except for CrArth-0.4, which uses its own selection mechanism.) Note that the ones with an explorative nature tend to have a large hit distance, while the hit distance in the more exploitative ones is smaller and declines rapidly to a low level as the run proceeds. Despite the fact that selection pressure has an impact on population diversity, explorative crossover manages to maintain a higher level of hit distance than exploitative crossover. StArth-1.2-2, StArth-1.5-2, StBLX-0.5-2 and StBLX-0.75-2 (see Figures 8(a) and 8(b), respectively) start the run with a high hit distance (note that as the selection pressure forces the population to converge, the hit distance decreases) but maintain a higher level than other exploitative crossovers. On the other hand, exploitative crossover starts the run with different hit distances (according to how exploitative it is based on the parameter settings), but as the run goes on, the hit distances drop to the same level regardless of the parameter settings. This is particularly noticeable with StHeur--2 (Figure 8(c)) but the same can be said about both other algorithms (Figures 8(a) and 8(b)). In these algorithms, all the operators, apart from mutation, have an exploitative nature and there is no operator that can maintain the diversity of the population. Moreover, as the diversity goes to the lowest level, the population activity becomes limited to a certain region where the new individuals (hits) are created not far away from their parents. Figure 8(d) shows the hit distance of three different algorithms that use the same crossover (Arth-0.4) but with different population-handling techniques. The algorithms act on the diversity of population in different ways, which lead the hit distance to be different even if the exact same crossover method is used.

To explore in greater depth the effects of selection pressure on the hit distance of different crossover methods, we measured the average of hit distance over different selection pressure levels, for the three crossover methods studied in this paper. Figures 9(a)–(c) show that the selection pressure has a large effect on explorative crossover methods, while the exploitative crossover methods undergo a lesser decline. The significant change in hit distance observed in explorative crossover methods is due to the conflicting biases of crossover and selection, where increasing the selection pressure overwhelms the bias of crossover. Figures 9(d)–(f) compare the hit distance for crossover with balanced, exploitative, and explorative parameter settings, respectively.

As indicated in Section 2.2, the total displacement characterises a population's dynamics by measuring the amount of change that population activities cause to the positions of node vectors. Populations that have a high total displacement exhibit a highly explorative nature, as they tend to transfer their activity regions around the search space and not to concentrate on certain regions. They usually have the ability to escape local optima, but they may risk abandoning the global optimum due to their tendency to widen and diversify the search. Figure 10 shows the effect of using different crossover methods and different level of selection pressure on the population's total displacement. Algorithms with explorative crossover have a greater total displacement. Also, increasing the selection pressure affects them more than exploitative crossovers.

Figure 11 uses the U-Matrix method (Ultsch and Siemon, 1990; which colours node vectors based on the average distance of their neighbours) to illustrate the movement of node vectors as they track population activity. The figure shows samples of two runs, StHeur-0.2-2 and StBLX-0.5-2, where the darker areas represent node vectors that are close to each other, while the light areas represent node vectors that point to different areas in the search space. We can see that StHeur-0.2-2 has a higher displacement, as the node vectors move faster toward the best-found region. More explanation on this particular example will be provided below.

Comparing population displacement and hit distance reveals information about different operator biases and the way they explore the search space. Figure 12 shows examples of how hit distance and displacement of different EAs change over a run. In Figure 12(a) we note that StHeur-0.2-2 has a low hit distance compared to StBLX-0.5-2. However, StHeur-0.2-2 displacement (Figure 12(d)) is greater than that of StBLX-0.5-2; and this changes in a fashion that implies a rapid move to a region in the search space, that is then followed by a sudden slowing down of its dynamics. This suggests that Heur-0.2 crossover has a stronger bias and tends to move the population in the direction of good observed fitness values, while BLX-0.5 tends to explore around activity regions longer and does not guide the search in such a directed way. The same explanation applies to StBLX-0.0-2 and StBLX-0.75-2 in Figures 12(b) and 12(e), while Figures 12(c) and 12(f) show that StHeur-0.1-2 and StArth-0.2-2 have approximately the same displacement behaviour despite the fact that they have a different hit distance.

### 4.3 Characterising Exploitation Behaviour

Experiments were carried out to characterise exploitation behaviour of the same set of algorithms by assessing the fitness-dependency of our collective-behaviour features. Figure 13 shows the fitness-dependency of the activity of different algorithms over time. When a population is more diverse and has a wide range of fitness values, which is the case of the initial population, algorithms show a higher level of exploitation as they tend to be more selective and focus on regions of good fitness. As the run progresses and the activity regions become smaller and so does the range of fitness values, our exploitation indicator declines toward a lower value which depends on the operators used and their parameter settings. Algorithms using heuristic crossover exhibit higher exploitation behaviour, while algorithms with diversity maintaining/generating techniques (fitness sharing, deterministic crowding, and random immigrants) show the lowest level of exploitation.

Figure 14 shows the effect of changing the selection pressure on the average of fitness-dependency of activity . In particular, Figures 14(a)–(c) show for the three crossover methods with four different settings for each, while Figures 14(d)–(f) compare algorithms with three different parameter settings. Note that for most of these algorithms, as we increase the selection pressure, the bias of selection tends to overcome the bias of crossover and the collective behaviour tends to become more exploitative. It is clear from Figure 14 that when crossover has an exploitative bias (i.e., Heur-0.1, …, Heur-0.4, Arth-0.2, Arth-0.4, BLX-0.0, and BLX-0.25), increasing selection pressure forces the algorithms to have similar exploitation behaviour regardless of the parameter settings of crossover. That is because these operators have a bias consistent with selection. On the contrary, crossover methods with explorative parameter settings tend to be more resistant to the bias of selection as shown in Figure 14(g).

Figures 15 and 16 depict the effects of the selection pressure on fitness-dependency of hit distance () and displacement (), respectively. These effects are more marked for fitness-dependency of activity than for the other features (for more detail, see Table 1), although our exploitation behaviour indicators are still increasing functions of the selection pressure. All crossover methods, approximately, have their bias overcome by the selection bias as the selection pressure increases. However, algorithms using crossover with explorative parameter settings are more resilient to the selection bias, as shown in Figures 15(g) and 16(g). Note also that in Figures 15(f) and 16(f), the fitness-dependency of hit distance and displacement of algorithms with exploitative parameter settings vary to a lesser extent.

Dependency measure . | Algorithm . | Minimum . | Maximum . | Range . |
---|---|---|---|---|

StArth | 0.3262 | 0.6583 | 0.3321 | |

StBLX | 0.2674 | 0.6597 | 0.3923 | |

StHeur | 0.3505 | 0.6648 | 0.3142 | |

StArth | 0.3775 | 0.5742 | 0.1967 | |

StBLX | 0.3711 | 0.5599 | 0.1888 | |

StHeur | 0.4812 | 0.5701 | 0.0889 | |

StArth | 0.3365 | 0.5502 | 0.2137 | |

StBLX | 0.2914 | 0.5461 | 0.2547 | |

StHeur | 0.3892 | 0.5547 | 0.1655 |

Dependency measure . | Algorithm . | Minimum . | Maximum . | Range . |
---|---|---|---|---|

StArth | 0.3262 | 0.6583 | 0.3321 | |

StBLX | 0.2674 | 0.6597 | 0.3923 | |

StHeur | 0.3505 | 0.6648 | 0.3142 | |

StArth | 0.3775 | 0.5742 | 0.1967 | |

StBLX | 0.3711 | 0.5599 | 0.1888 | |

StHeur | 0.4812 | 0.5701 | 0.0889 | |

StArth | 0.3365 | 0.5502 | 0.2137 | |

StBLX | 0.2914 | 0.5461 | 0.2547 | |

StHeur | 0.3892 | 0.5547 | 0.1655 |

Table 1 summarises our exploitation indicators for standard EAs. The table shows the minimum, maximum, and range of each of our three fitness-dependency averages for the three crossover methods over all the possible settings and all selection pressure levels. From the table, we see that has the lowest minimum value and the biggest maximum value, thus, the largest range, among all dependency measures. On the other hand, we see that has the highest minimum value and the smallest range among all dependency measures. These two observations tell us that the activity is more driven by selection because increasing the selection pressure increases the fitness-dependency, while the hit distance is mainly dependent on crossover, as changing the selection pressure does not have a significant effect on fitness dependency. Displacement fitness-dependency ranges are in between these extremes, which suggests that both selection and crossover have an effect on fitness-dependency. By looking at each of the table parts separately, we note that heuristic crossover (Heur-) has the highest fitness-dependency, which is consistent with its fitness-exploiting nature in handling the genetic material. We also see that hit-distance fitness-dependency of the same crossover is the least affected by selection pressure. On the other hand, blending crossover (BLX-) has the lowest minimum fitness-dependency for all the three exploitation indicators, which can be attributed to its stochastic nature (although, by looking at the last column of the table, we see that selection pressure increases the fitness-dependency of BLX-).

## 5 Conclusions

In this paper, we have presented a model to represent the collective dynamic behaviour of population-based algorithms. As we saw in Section 1, previous models used either the genotypic or the phenotypic representation to capture such behaviour. Many such models are tailored to evolutionary algorithms (usually a simplified version of them) and have limited practical use. Here we have proposed an algorithm-independent empirical approach that uses a sequence of topological genotypic and phenotypic representations (at the same time) to capture the population collective behaviour. Our models allow one to represent the dynamics at different levels of abstraction, via the aggregation of the most elementary components of the model into more coarse-grained entities.

In addition, a number of measures have been proposed to represent different features of the collective dynamic behaviour of a population. These fall within two categories: activity-related features and distance-related features. Activity-related measures describe the way an algorithm distributes its activities around regions in the search space while distance-related measures describe the population's dynamics from a geometric viewpoint. These measures can be used to qualitatively characterise the exploration/exploitation behaviour of an algorithm. Also, a set of fitness-dependency measures are introduced that can be used as indicators of an algorithm's exploitation behaviour. These measures assess the dependency of different population activity on fitness values, in order to characterise the exploitation tendency of an algorithm. We used an entropy-based measure, namely the uncertainty coefficient, to quantify the dependency of different emergent features of the collective behaviour on the fitness.

The proposed model has been implemented using a modified version of a self-organising map. SOMs are used to capture the distribution of the initial population, and track its moves as the algorithm operates on solving the problem and creates new individuals in different areas of the search space. There is nothing special about this choice. Other topology-preserving algorithms might have been used. Also, for a number of features, one does not need a topology, and so a vector-quantisation algorithm would suffice. In the future, we will investigate the relative benefits and drawbacks of alternative techniques.

In this paper, we have studied the collective dynamic behaviour of EAs focusing on the effect of using different crossover methods, various parameter settings, and different population-handling techniques. The impact of combining consistent or conflictive biases on the emergent features of the collective behaviour has been studied. However, we believe our methods would also work well with other population-based algorithms. In the future, we will also explore this avenue of research.

The results of experimenting with different EAs have shown that the proposed measures are able to capture new aspects of the population dynamics and reveal important features of the collective behaviour. Experiments have also shown that our proposed measures of exploitation are useful in quantifying the exploitation behaviour of populations. The proposed measures enable a new level of algorithm comparison and provide tools to investigate the effect of combining different biases on the resulting collective behaviour of a population.

Our model is an effort to contribute to our understanding of the collective dynamic behaviour of population and the exploration/exploitation phenomena. It provides an empirical analysis tool that allows us to represent the dynamic behaviour at different levels of abstraction and to identify and measure different emergent features of the dynamic behaviour, thus enabling us to better describe and, potentially, tune that behaviour. By studying the dependence between the emergent features, we tried to define what constitutes exploitation (or exploration) and to understand the effect of operators with different biases on the collective dynamic behaviour of populations.

This model can also be of use by EA practitioners for analysing the effects of operators and for tuning their control parameters. It can help characterise the way contradictory/consistent biases combine together to produce the collective behaviour, or how populations behave differently according to the nature of the fitness landscape they are operating on. This can help identify what make certain algorithms fail/succeed in searching certain fitness landscapes. All these issues will be investigated in future work.

## References

## Appendix A: EA Parameter Settings, Naming Convention, Techniques, and Benchmark Problems

In this work, we used generational evolutionary algorithms with four population handling techniques: standard, fitness sharing, deterministic crowding, and random immigrants. Except for deterministic crowding, which uses its own selection mechanism, we used tournament selection with different tournament sizes to provide different levels of selection pressure. Table 2 summarises the common parameter settings among all the EAs that have been used in our experiments.

Parameter . | Value . |
---|---|

Population size | 100 |

Crossover rate | 0.75 |

Mutation rate | 0.2 |

Mutation step | [0, 1.0] |

Individual length | 5 |

Total number of evaluations per run | 10,000 |

Parameter . | Value . |
---|---|

Population size | 100 |

Crossover rate | 0.75 |

Mutation rate | 0.2 |

Mutation step | [0, 1.0] |

Individual length | 5 |

Total number of evaluations per run | 10,000 |

We used a naming standard for EAs to capture four facts about an algorithm: population-handling technique, crossover method, crossover control parameter, and tournament size (selection pressure). The algorithm name has the format **PPCCCC-r-t** where **PP** denotes the population-handling technique as follows: **St**: standard EA, **FS**: fitness sharing EA, **Cr**: deterministic crowding EA, and **RI**: random immigrant EA. **CCCC** refers to the crossover method, **Arth**: arithmetic crossover, **BLX**: blending crossover, and **Heur**: heuristic crossover. **r** is a real number representing the control parameter of the crossover method used, and **t** is an integer representing the tournament size. For example, FSArth-0.4-4 is an EA using fitness sharing and the arithmetic crossover with control parameter 0.4 and tournament size 4.

In addition to standard generational EA, we used three population-handling techniques. (1) Fitness sharing (Goldberg and Richardson, 1987) is a fitness scaling mechanism that alters only the fitness evaluation stage of an EA and can be used with any other techniques. In this method, similar individuals (located in the vicinity of each other on the fitness landscape) have to share the fitness among them. Thus, the number of individuals that can reside in an area is limited by the fitness of that area. A full description can be found in Back et al. (1997). (2) Crowding techniques (DeJong, 1975) insert new individuals into the population by replacing similar individuals. However, replacement errors may prevent crowding techniques from maintaining individuals on different but nearby peaks. Deterministic crowding (Mahfoud, 1992, 1995) is designed to minimise the replacement error. A full description of the algorithm and pseudocode can be found in Back et al. (1997). (3) The random immigrants mechanism (Cobb and Grefenstette, 1993) replaces part of the population with randomly created individuals. In this work, we replace 10% of the worst individuals in a population by random ones. We perform this replacement at the end of every EA round.

We tested the algorithms on 5-dimensional real-valued optimisation problems obtained by composing a random combination of 10 basic functions taken from the following family: Sphere, Ranstrigin, and Griewank. These functions were randomly generated, shifted, rotated, and combined. In our experiments, we conducted 100 runs for each algorithm setting (e.g., crossover method and its control parameter, tournament size,, etc.). Each run operated on a random fitness landscape. Full details on the benchmark functions and how to compose them can be found in Liang et al. (2005).

## Appendix B: Crossover Operators for Real-Coded Evolutionary Algorithms

Three crossover methods have been used in this paper. They all return two offspring and all have a control parameter that can be used to tune their exploration/exploitation behaviour. These methods are:

**Arithmetic Crossover.**(Arth-) (Michalewicz, 1996) Two offspring individuals,*C*_{1}and*C*_{2}, are produced from two parents,*P*_{1}and*P*_{2}, as follows: and .**Blending Crossover.**(BLX-; Eshelman and Schaffer, 1992) Two offspring individuals are produced by randomly (uniformly) generating values for their genes. Suppose that*p*^{i}_{1}and*p*^{i}_{2}are the*i*th parameter of parents*P*_{1}and*P*_{2}, respectively, the corresponding parameter,*c*, of offspring^{i}_{k}*C*is randomly chosen from the interval , where , , and_{k}*I*=*p*−_{max}*p*._{min}**Heuristic Crossover.**(Heur-; Herrera and Lozano, 1996; Wright, 1990) This method that creates one offspring individual around the parent with the highest fitness. Here we modified it slightly so as to produce two offspring. Suppose that we have two parent individuals,*P*_{1}and*P*_{2}, and*P*_{1}is the one with higher fitness, then the two offspring individuals,*C*_{1}and*C*_{2}, are created as: and .

## Notes

^{1}

An initial formulation and some results of our model were presented in Turkey and Poli (2012).

^{2}

SOMs have previously been used to improve the performance of evolutionary algorithms by enhancing the search strategy and avoiding genetic drift (Amor and Rettinger, 2005).

^{3}

In this work, we assume that the algorithm is operating on real-coded problems. However, a grid of nodes for different types of representations can be defined.

^{4}

Please note the difference between *s ^{t}*, which is a set of regions (set of sets), and

*S*, which is the set of all nodes in these regions. Formally, .

^{t}^{5}

We conducted our experiments on a grid of size 25 × 25.

^{6}

By standard evolutionary algorithm, we mean a generational evolutionary algorithm that uses no extra techniques. Full parameters settings are given in Appendix A.

^{7}

Of course, growth in those regions does not mean that an algorithm also improves the best-so-far fitness of the search.