**1. Introduction**

Memetic algorithms (MAs) are a class of hybrid evolutionary algorithms (EAs) and have been widely applied to many complex optimization problems [1], such as scheduling problems [2–4], combinatorial optimization problems [5–7], multi-objective optimization problems [8,9], and multimodal optimization problems [10,11]. Population-based evolutionary approaches, such as genetic algorithms (GAs) or particle swarm optimization (PSO) [12], are able to detect in a fast way the main regions of attraction, while their local search abilities represent a major drawback [13] from the point of view of solution accuracy and of the convergence behaviour, especially when applied to multimodal problems. MAs [1,14] have recently received attention as effective meta-heuristics to improve generic EA schemes by combining these latter with local search (LS) procedures [13,15]. In MAs, EA operations are employed for global rough exploration, and LS operators are used to execute further exploitation of single EA individuals. This allows one to enhance the EA's capability to solve complicated problems. In the context of PSO [12], the LS procedure can be regarded as an iterative local move of the particle. The scheme has been successfully applied to improve different meta-heuristics such as simulated annealing [16], GAs [17], and PSO [18]. A memetic algorithm that hybridizes PSO with LS, called MPSO, is proposed in [11] for locating multiple global and local optima in the fitness landscape of multimodal optimization problems (MMOPs). In this work, the local search is designed by means of two different operators, cognition-based local search (CBLS) and random walk with direction exploitation (RWDE), but no adaptation technique of the local search operators is employed.

Although many GA-based MAs and other EA approaches employing meme evolution have been proposed in the literature [13,19], the PSO memetic algorithms proposed so far fail to employ meme evolution, probably because of the difficulty of designing the descriptors of the local search operators, which are usually statically tailored to the problem, or the class of problems, at hand [7,8,10,11,20]. Indeed, they exhibit, at most, only some form of temporary adaptation that is not preserved or transmitted through the generations as meta-Lamarkian memes [21].

Since meme evolution is able to provide adaptivity of the local search operators, and many research results have shown that the advantages of the introduction of a local search in EAs, namely in PSO, are generalized, we have set as a research objective the investigation of the possibility of introducing memetic evolution of LS operators in a PSO context. The idea was to design a framework where memes are not selected in a preliminary evolutionary phase, but where they evolve and adapt during the search.

The research resulted in the coevolving memetic particle swarm optimization (CoMPSO) algorithm presented in this work. CoMPSO is an algorithm for numerical optimization in a hybrid sea and allows one to evolve local search operators in the form of memes associated with the PSO particles exploring a continuous and/or discrete hybrid search space. A special feature of local search operators, defined by the memes, is that they are applied in a probabilistic way. A meme is applied to the personal best position of the associated particle with a given probability and every given number of iterations. The basic idea is that memes are evolved themselves by means of the PSO-MEME evolution scheme, a classical PSO scheme modified to manage meme vectors composed of hybrid continuous/discrete components. Indeed, the PSO-MEME scheme distinguishes the evolution of continuous and discrete components where a probabilistic technique is applied.

The rest of the paper is structured as follows. Some general issues concerning MAs and PSO are recalled in Section 2.1, while the coevolving memetic PSO scheme is introduced in Section 2.2 by describing the proposed co-evolving particle-meme scheme and a probabilistic technique for PSO evolution of discrete integer components. Experimental results and comparisons with classical PSO and state-of-the-art static memetic PSO are presented and discussed in Section 3. Finally, in Section 4, a discussion of the advantages of the proposed approach and a description of future lines of research conclude the paper.

#### **2. Co-Evolution of Memes in PSO: Models and Methods**

In this section, we introduce the architecture of CoMPSO by first recalling some concepts of memetic algorithms and PSO, respectively in Sections 2.1.1 and 2.1.2. A detailed description of the CoMPSO scheme is then provided in Section 2.2, focusing on the algorithmic aspects introduced to guarantee diversity and on the representation for memes in the discrete, continuous, and hybrid space search, and finally describing and explaining the motivations for meme co-evolution strategies.

## *2.1. Memetic Algorithms and Particle Swarm Optimization*

## 2.1.1. Memetic Algorithms

MAs, first proposed by Moscato [14,22], represent a recent growing area of research in evolutionary computation. The MA meta-heuristics essentially combine an EA with local search techniques that can be seen as a learning procedure that makes individuals capable of performing local refinements. MAs take inspiration from the Dawkins notion of "meme" [23] representing a unit of cultural evolution that can exhibit refinement.

MAs are usually distinguished into three main classes. First generation MAs [22] are characterized by a single fixed meme; i.e., a given local search operator is applied to candidate solutions of the main EA according to different selection and application criteria (see, for example, [20]). The criteria used in the selection of candidates and the frequency of the local search application represent the parameters of MAs. Second generation MAs, also called "meta-Lamarckian MAs" [21,24], reflect the principle of memetic transmission and selection; i.e., they are characterized by a pool of given local search operators (memes) that compete, based on their past merits, to be selected for future local refinements. The choice of the best memes can be based on different criteria, such as the absolute fitness of the candidate solutions associated with the meme or the fitness improvement due to past applications. Finally, the most recent generation of MA schemes are the co-evolving MAs [10,15,25,26], which introduce the idea of meme evolution. In this case, the pool of memes is not known a priori, but the memes are represented and evolved by means of EA methods similar to those applied for candidate solutions. The CoMPSO scheme proposed in this paper falls in this latter class.

## 2.1.2. Particle Swarm Optimization

Particle swarm optimization (PSO), originally introduced in [18], is a meta-heuristic approach to continuous optimization problems that is inspired by the collective behavior observed in flocks of birds. Further variants of PSO have been proposed in [12,27–33].

In PSO, a swarm of artificial entities—the so called *particles*—navigates the search space, aiming at optimizing a given objective/fitness function *f* : Θ → R, where Θ ⊆ R*<sup>d</sup>* is the feasible region of the space. The set *P* = {*p*1, ... , *pn*}, composed of *n* particles, is endowed with a neighborhood structure; i.e., a set *Ni* ⊆ *P* of neighbors is defined for each particle *pi*. Every particle has a position in the search space, which is iteratively evolved using its own search experience and the experience of its neighbors. Therefore, PSO adopts both cognitive and social strategies [34] in order to focus the search on the most promising areas.

At every iteration *t* ∈ N, any particle *pi* is composed of the following *d*-dimensional vectors: the *position* in the search space *xi*,*t*, the *velocity vi*,*t*, the *personal best* position *bi*,*<sup>t</sup>* visited so far by the particle, and the *neighborhood best position gi*,*<sup>t</sup>* visited so far among the particles in *Ni*.

Both positions and velocities of the PSO particles are randomly initialized. Then they are iteratively updated until a stop criterion is met (e.g., the allowed budget of fitness evaluations has been consumed), according to the following *move equations* (first studied in [18]):

$$
\omega\_{i,t+1} = \omega v\_{i,t} + q\_1 r\_{1,t} (b\_{i,t} - \mathbf{x}\_{i,t}) + q\_2 r\_{2,t} (g\_{i,t} - \mathbf{x}\_{i,t}), \tag{1}
$$

$$\mathbf{x}\_{i,t+1} = \mathbf{x}\_{i,t} + v\_{i,t+1} \,. \tag{2}$$

In Equation (1), the weight *ω* is the inertial coefficient, *ϕ*1 and *ϕ*2 are the acceleration factors, while the two random vectors *<sup>r</sup>*1,*t*,*r*2,*<sup>t</sup>* ∈ R*<sup>d</sup>* are uniformly distributed in [0, 1] *d*. The first term of Equation (1)—the *inertia* or momentum—is a memory of the trajectory followed by the particle so far. Its aim is to prevent drastic changes in the search performed by the particle. The second term is the *cognitive component* and represents the tendency of the particle to return to the best position ever visited by itself. Finally, the third term—the *social component*—models the contribution of the neighbors in the particle's trajectory.

After every position update, the fitness of *xi*,*<sup>t</sup>* is evaluated. If it is better than those of the previous personal and social best positions, these are updated accordingly.

Usually a fully connected topology is adopted as a neighborhood graph among the particles, so only one global social best *gt* is needed instead of *n* local copies of it. It should be noted that sometimes the particles can exceed the bounds of the feasible search space. A common practice to address this issue, adopted also in our work, is to restart the out-of-bounds components of a particle in a position randomly chosen between the old position and the exceeded bound [35,36].

## *2.2. Coevolving Memetic PSO*

Coevolving memetic PSO (CoMPSO) is a hybrid algorithm that combines PSO with the evolution of local search memetic operators that automatically adapt to the objective function landscape.

CoMPSO evolution can be seen as the combination of two co-evolving populations: the particles and the memes. The PSO particles represent the candidate solutions and evolve by means of the usual PSO laws (see Section 2.1.2). The memes represent the local search operators applicable to particles and are evolved by a modified PSO scheme that employs an innovative technique designed to deal with the discrete domains present in the meme representation. Furthermore, another difference with respect to the standard PSO scheme is the introduction of a diversity control mechanism that prevents the premature convergence of the particle population to a local optimum. Other example of the hybridization of PSO can be found in the area of pattern recognition of strings [37–40].

In the following, the general CoMPSO scheme (Section 2.2.1), the diversity control mechanism of particles evolution (Section 2.2.2), the meme representation (Section 2.2.3), meme evolution (Section 2.2.4), and the novel technique for managing discrete components in meme PSO (Section 2.2.5) are introduced.

## 2.2.1. The CoMPSO Scheme

In the CoMPSO general scheme, a population of *n* particles, i.e., *P* = {*p*1, ... , *pn*}, navigates the search space following PSO dynamics while trying to find the optimum of the given fitness function *f* . At each iteration, a local search is possibly applied to a subset of particles in order to improve their fitness. This local search phase is realized by a meme population, i.e., *M* = {*<sup>m</sup>*1, ... , *mn*}, of local search operators. Each *mi* is associated with particle *pi* and is evolved using a PSO-like approach on the meme space.

Particles in *P* are encoded using the usual *d*-dimensional vectors of the PSO scheme (see Section 2.1.2), i.e., the position vector *x*(*p*) *i*,*t* , the velocity vector *v*(*p*) *i*,*t* , and the personal best vector *b*(*p*) *i*,*t* , other than the common global best *g*(*p*) *t* . In this work, a complete neighborhood graph has been adopted. Similarly, a generic meme *mj* ∈ *M* is encoded by the hybrid vectors *x*(*m*) *j*,*t* , *v*(*m*) *j*,*t* , *b*(*m*) *j*,*t* , and *g*(*m*) *t* , i.e., vectors that combine both discrete and continuous components as described in Section 2.2.5.

The general CoMPSO scheme is described by the pseudo-code reported in Algorithm 1.

The populations of particles *P* and memes *M* are randomly initialized in their respective feasible regions. During main-loop iterations, each particle *pi* is evolved following the scheme described in Section 2.1.2. *LOCAL\_SEARCH* activates the memetic part of CoMPSO. The local search operators, defined by the memes in *M*, are applied in a probabilistic way: *mi* is applied to the personal best position of particle *pi* with a probability *γ* at every *φ* iteration, and the local search application is denoted by *mi*(*pi*) or equivalently by *x*(*m*) *i* (*b*(*p*) *i* ). Furthermore, for every iteration, the local search is also applied to the global best particle in position *g*(*p*) *t* . In this case, the global best particle and its respective meme are indicated by *pg* and *mg*. Before being applied, the memes are evolved by means of the *PSO\_MEME* evolution scheme described in Section 2.2.4. If the meme application leads to a fitness improvement, the new candidate solution obtained replaces both the personal best that the particle current position. Finally, the global best is updated accordingly.

**Algorithm 1** CoMPSO Pseudo-Code.

1: **procedure** COMPSO 2: *t* ← 0 3: *P* ← INITIALIZE\_PARTICLES() 4: *f*(*P*) ← EVALUATE\_FITNESS() 5: *M* ← INITIALIZE\_MEMES() 6: **while** termination criterion is not met **do** 7: *t* ← *t* + 1 8: **if** DIV\_CONTROL() **then** 9: DIV\_RESTORE() 10: **end if** 11: **for all** *pi* ∈ *P* **do** 12: *v*(*p*) *i*,*t* ← UPDATE\_VELOCITY(*v*(*p*) *i*,*t*−1, *x*(*p*) *i*,*t*−1, *b*(*p*) *i*,*t*−1, *g*(*p*) *<sup>t</sup>*−<sup>1</sup>) 13: *x*(*p*) *i*,*t* ← UPDATE\_POSITION(*x*(*p*) *i*,*t*−1, *v*(*p*) *i*,*t* ) 14: *f*(*x*(*p*) *i*,*t* ) ← EVALUATE\_FITNESS(*x*(*p*) *i*,*t* ) 15: *b*(*p*) *i*,*t* ← UPDATE\_PERSONAL\_BEST(*b*(*p*) *i*,*t*−1, *x*(*p*) *i*,*t* ) 16: **end for** 17: *g*(*p*) *t* ← UPDATE\_GLOBAL\_BEST() 18: **if** LOCAL\_SEARCH(*<sup>γ</sup>*, *φ*) **then** 19: **for all** *pi* ∈ *P* **do** 20: PSO\_MEME(*mi*) 21: *mi*(*pi*) 22: **end for** 23: **end if** 24: PSO\_MEME(*mg*) 25: *mg*(*pg*) 26: *g*(*p*) *t* ← UPDATE\_GLOBAL\_BEST() 27: **end while** 28: **return** getBestSolution 29: **end procedure**

2.2.2. Diversity Control

A diversity control mechanism has been introduced in order to avoid stagnation or premature convergence to local optima. When premature convergence is detected, a subset of particles whose positions will be reinitialized is selected. The method consists of two main components: a *diversity measure δ* and a *diversity restore mechanism* invoked when the population diversity becomes too low according to measure *δ*.

The diversity measure is computed on the fitness values of the particle population according to

$$\delta(P\_t) \leftarrow \text{std}(f(P\_t))\tag{3}$$

where std(*f*(*Pt*)) represents the standard deviation of the fitness values of particle population *P* at time *t*. It should be noted that, although genotypic distances between particles seem to be in principle more appropriate than fitness values, it has been experimentally observed that the use of these latter indicators provides a good diversity measure with the additional property of a more efficient computation (this is due to the fact that the computation is performed on R and not on R*<sup>d</sup>*).

The diversity *<sup>δ</sup>*(*<sup>P</sup>*0) of the initial random population *P*0 is employed to compute the threshold value *τ* = 0.2 · *<sup>δ</sup>*(*<sup>P</sup>*0) used in the routine *DIV\_CONTROL* in order to recognize population convergence or stagnation. This threshold value is computed basing on the unbiased particle population at time 0 since it is the only one generated in a pure random way.

Therefore, at each iteration, in the case that *δ*(*Pt*) < *τ*, the *DIV\_RESTORE* procedure is invoked. This routine randomly restarts the positions of the worst fit 50% of the population. Only the current positions of particles are restarted, while velocities and personal best values are kept unchanged.

## 2.2.3. Meme Representation

The local search operators adopted in CoMPSO are a generalization of the random walk (RW) technique [41–43]. RW is an iterative and stochastic local optimization method. Let *yt* be the candidate (local) minimum at the *t*-th iteration. Then the new value *yt*+<sup>1</sup> is computed according to

$$y\_{t+1} \leftarrow \begin{cases} y\_t + w\_t z\_t & \text{if } f(\mathbf{x}\_t + w\_t z\_t) < f(\mathbf{x}\_t) \\ y\_t & \text{otherwise} \end{cases} \tag{4}$$

where *wt* is a scalar step-length, and *zt* is a unit-length random vector. The step length *wt* is initialized to a given value *w*0 and is halved at the end of every iteration in the case that the fitness has not been improved. The process is repeated until a given number *q* of iterations is reached.

In our generalization, other than the previously described parameters *w*0 and *q*, two other parameters are introduced: the ramification factor *b* and the number of chosen successors *k* (with *k* ≤ *b*). Briefly, the idea is to expand RW in a (*b*-*k*-bounded) breadth search conversely from the pure depth-first style of the original RW. Initially, *k* copies of the starting (or seed) point of RW are generated. Then *b* new points are generated from the current *k* ones following Rule (4). Therefore, from these *b* intermediate points, the *k*-fittest ones are chosen as next-iteration current points, and the process is iterated until the deepness parameter *q* is reached.

Using this local search scheme, a generic meme is represented by means of the following four parameters: a real value *w*0 and the three integers *b*, *k*, and *q*. Furthermore, for each parameter, a range of admissible values has been experimentally established: *w*0 ∈ [0.5, <sup>4</sup>], *b* ∈ [1, 8] ∩ N, *k* ∈ [1, *b*] ∩ N, and *q* ∈ [4, 16] ∩ N.

## 2.2.4. Meme Evolution

As introduced in the previous section, a meme is represented in the hybrid (continuous/discrete) meme space by the 4-ple *m* = (*<sup>w</sup>*0, *b*, *k*, *q*). Meme PSO evolution proceeds asynchronously with respect to particle evolution, this is due to the fact that memes are evolved only before they are applied and according to probability *γ* and frequency *φ* described above.

Memes, like particles, at each iteration *t*, are characterized by a position *x*(*m*) *i*,*t* , i.e., the representation of the meme in the meme space, a velocity *v*(*m*) *i*,*t* (only for the continuous component), a personal best position *b*(*m*) *i*,*t* , and a global best meme *g*(*m*) *t* . Two alternative functions have been considered as meme fitness *f* (*m*): (1) the "absolute" fitness *f* of the particle *pi* associated to meme *mi*, and (2) the fitness improvement Δ*f* realized with the application of the meme *mi* to particle *pi*. Some preliminary experiments have shown that the two approaches do not significantly differ.

Memes are evolved by a classical PSO scheme modified to manage the meme vectors composed by hybrid continuous/discrete components. The *PSO\_MEME* scheme distinguishes the evolution of continuous components (according to Update Rules (1) and (2)) and that of discrete components where a probabilistic technique, presented in the next subsection, is applied. For the sake of completeness, *PSO\_MEME* pseudo-code is reported in Algorithm 2.


1: **procedure** PSO\_MEME(*mi*) 2: **for** each continuous dimension *j* **do** 3: *v*(*m*) *i*,*j* ← UPDATE\_CONTINUOUS\_VELOCITY(*v*(*m*) *i*,*j* , *x*(*m*) *i*,*j* , *b*(*m*) *i*,*j* , *g*(*m*) *j* ) 4: *x*(*m*) *i*,*j* ← UPDATE\_CONTINUOUS\_POSITION(*x*(*m*) *i*,*j* , *v*(*m*) *i*,*j* ) 5: **end for** 6: **for** each discrete dimension *j* **do** 7: *x*(*m*) *i*,*j* ← UPDATE\_DISCRETE\_POSITION(*x*(*m*) *i*,*j* , *b*(*m*) *i*,*j* , *g*(*m*) *j* ) 8: **end for** 9: *f* (*m*)(*mi*) ← EVALUATE\_MEME\_FITNESS(*mi*, *pi*) 10: *b*(*m*) *i* ← UPDATE\_MEME\_PERSONAL\_BEST(*b*(*m*) *i* , *x*(*m*) *i* ) 11: *g*(*m*) ← UPDATE\_MEME\_GLOBAL\_BEST(*g*(*m*), *x*(*m*) *i* ) 12: **end procedure**

2.2.5. Meme PSO for Discrete Domains

The meme discrete components are managed using a probabilistic technique that, by exploiting the total order of the integer domains here considered, simulates the classical PSO dynamics for continuous domains.

For each meme and for each discrete component domain *Dj* = {*dj*,1, ... , *dj*,*<sup>r</sup>*}, an appropriate probability distribution *PDj* over the values of *Dj* is built, then the new position *xj* for component *j* is computed according to a randomized roulette wheel tournament among values in *Dj* and by using probabilities *PDj* .

As described in the following, the probability distribution is set in a way that simulates the properties of the classical continuous PSO [44].

**Probability Distribution on Discrete Components** Let *dx*, *db*, *dg* ∈ *Dj* represents the values: the current meme position, the meme personal best position, and the global best meme position for component *j* of the meme encoding, then the distribution *Dj* is defined as a mixture of four discrete probability distributions:


Each triangular distribution is determined by the center value *dk*, an *amplification* factor *α* > 1 and a *width* 2 · *λ* + 1. The probability of *dk* is amplified by the factor *α*, i.e., the previous probability *PDj*(*dk*) is multiplied by *α*, while the 2 · *λ* values of domain *Dj* located around *dk*, i.e., *dk*−*λ*, ... , *dk*−<sup>1</sup> and *dk*+1, ... , *dk*+*λ*, are amplified by a *smoothing* factor *β* = *α λ*+1 · |*λ* + 1 − *s*|, where *s* is the distance of the value from the center *dk*. Moreover, it should be noted that the probabilities of each triangular distribution are normalized in order to sum up to 1.

The amplification factor *α* are set by using the PSO parameters; that is 1 + *ω* for the center *dx*, 1 + *ϕ*1 for *db*, and 1 + *ϕ*2 for *dg*.

**Discrete position update.** The new position of a discrete component *xj* is computed by a random roulette wheel tournament that uses the probability distribution above.

The technique for discrete components based on probability distribution values can be interpreted as initially considering all the values as equally probable, except for *dx*, *db*, and *dg*, whereas the amplification factors based on classical PSO parameters *ω*, *ϕ*1, *ϕ*2 confer a greater probability. Moreover, the smoothness factors, *β*s, also amplify the probabilities of centers neighbors, i.e., values close to the centers.

In other words, amplification and smoothing factors implement, for the discrete components, the probabilistic counterparts of the typical behavior of PSO velocity, i.e., the tendency to remain in the current position (inertia) and the tendency to move toward a personal and global best (cognitive and social dynamics). It should be noted that a similar idea is present in [45].
