Abstract

Quantum-behaved particle swarm optimization (QPSO), motivated by concepts from quantum mechanics and particle swarm optimization (PSO), is a probabilistic optimization algorithm belonging to the bare-bones PSO family. Although it has been shown to perform well in finding the optimal solutions for many optimization problems, there has so far been little analysis on how it works in detail. This paper presents a comprehensive analysis of the QPSO algorithm. In the theoretical analysis, we analyze the behavior of a single particle in QPSO in terms of probability measure. Since the particle's behavior is influenced by the contraction-expansion (CE) coefficient, which is the most important parameter of the algorithm, the goal of the theoretical analysis is to find out the upper bound of the CE coefficient, within which the value of the CE coefficient selected can guarantee the convergence or boundedness of the particle's position. In the experimental analysis, the theoretical results are first validated by stochastic simulations for the particle's behavior. Then, based on the derived upper bound of the CE coefficient, we perform empirical studies on a suite of well-known benchmark functions to show how to control and select the value of the CE coefficient, in order to obtain generally good algorithmic performance in real world applications. Finally, a further performance comparison between QPSO and other variants of PSO on the benchmarks is made to show the efficiency of the QPSO algorithm with the proposed parameter control and selection methods.

1.  Introduction

The particle swarm optimization (PSO) algorithm is a population-based optimization technique, originally developed by Kennedy and Eberhart in 1995. It was motivated by the social behavior of bird flocking or fish schooling and shares many similarities with evolutionary computation techniques. A PSO system is initialized with a population of random solutions and searches for optima by updating generations. However, unlike evolutionary algorithms, PSO has no evolution operators such as crossover and mutation. In PSO, the potential solutions, called particles, fly through the problem space by following their own experiences and the current best particles. It has been shown that the PSO algorithm is comparable in performance with and may be considered as an alternative method to evolutionary algorithms (Angeline, 1998a).

During the last decade, PSO gained increasing popularity since it can get better results for optimization problems in a faster and cheaper way compared with other methods, but has fewer parameters to be adjusted. In order to investigate in detail the mechanism of PSO, a lot of theoretical analyses have been done on the algorithm (Kennedy, 1998; Ozcan and Mohan, 1999; Clerc and Kennedy, 2002; van den Bergh, 2002; Shi and Eberhart, 1998b; Trelea, 2003; Emara and Fattah, 2004; Gavi and Passino, 2003; Kadirkamanathan et al., 2006; Jiang et al., 2007; Poli, 2008). These theoretical analyses were focused on the behavior of the individual particle, which is essential to the understanding of the search mechanism of the algorithm and to parameter selection. For example, Kennedy (1998) carried out an analysis of simplified particle behavior and showed the different trajectories of particles for a range of design choices aiming to gain some insights into the behavior of particles through simulations, and Clerc and Kennedy (2002) undertook the first formal analysis of the particle trajectory and of the stability properties of the algorithm.

Besides the theoretical analyses, there has been a considerable amount of work done in developing the original version of PSO through empirical simulations. In order to accelerate the convergence of the particle, Shi and Eberhart (1998a) introduced the concept of an inertia weight into the original PSO, and Clerc (1999) proposed an alternative version of PSO incorporating a parameter known as the constriction factor which should replace the restriction on velocities. Some researchers employed the operations of other evolutionary algorithms in PSO to enhance its performance, the most important work being done by Angeline (1998b). Another general form of particle swarm, referred to as the lbest model, was first proposed by Eberhart and Kennedy (1995), and then was studied in depth by many other researchers in order to find other topologies to improve the performance of PSO (Suganthan, 1999; Kennedy, 1999, 2002; Mendes et al., 2004; van den Bergh and Engelbrecht, 2004; Janson and Middendorf, 2005; Liang and Suganthan, 2005; Mohais et al., 2005; Parrott and Li, 2006). Bratton and Kennedy defined a standard PSO, which is an extension of the original PSO algorithm while taking into account the previous developments that can improve the performance of the algorithm (Bratton and Kennedy, 2007).

Some researchers have attempted to experiment with various ways to simulate the particle trajectory by directly sampling, using a random number generator, from a distribution of some theoretical interest, and thus have proposed many probabilistic PSO algorithms (Kennedy, 2003, 2004, 2006; Secrest and Lamont, 2003; Krohling, 2004; Sun, Feng, et al., 2004; Sun, Xu, et al., 2004; Sun et al., 2005; Krohling and Coelho, 2006; Richer and Blackwell, 2006), the most popular being the bare-bones PSO (BBPSO) family. In BBPSO, each particle has no velocity vector and its new position is sampled around a supposed good one according to a probability distribution. In the original BBPSO, this distribution is Gaussian (Kennedy, 2003). A moment analysis was proposed to determine the characteristics of the sampling distribution of BBPSO and other PSO algorithms (Poli, 2008).

The focus of this paper is on a probabilistic algorithm, quantum-behaved particle swarm optimization (QPSO), which was proposed by Sun, Feng et al. (2004). The inspiration of QPSO came from quantum mechanics and the trajectory analysis of PSO (Clerc and Kennedy, 2002). The trajectory analysis showed that each particle in PSO oscillates around and converges to its local attractor, or to put it in other words, each particle is in a bound state. In QPSO, the particle is assumed to have quantum behavior and to be in a bound state, and is further assumed to be attracted by a quantum potential well centered on its local attractor, thus having a new stochastic update equation for its position (Sun, Feng et al., 2004). Later, a global point known as the mean best position was introduced into the algorithm in order to enhance the global search ability of the QPSO algorithm (Sun, Xu et al., 2004; Sun et al., 2005).

The QPSO algorithm essentially belongs to the BBPSO family, but samples the new position with a double exponential distribution. Besides, its update equation uses an adaptive strategy and has fewer parameters to be adjusted, leading to a good performance of the algorithm as an overall result. The QPSO algorithm has aroused the interest of many researchers from different communities. It has been shown to successfully solve a wide range of continuous optimization problems. Among these applications, it has been used to tackle the problems of constrained optimization (Sun et al., 2007), multi-objective optimization (Omkara et al., 2009), neural network training (Li et al., 2007), electromagnetic design (Mikki and Kishk, 2006; Coelho and Alotto, 2008), semiconductor design (Sabata et al., 2009), clustering (Sun et al., 2006), system identification (Gao, 2008), engineering design (Coelho, 2008, 2010), image processing (Lei and Fu, 2008), power systems (Coelho and Mariani, 2008; Sun et al., 2009; Zhang, 2010; Sun and Lu, 2010), bioinformatics (Chen et al., 2008; Cai et al., 2008), to name only a few.

In addition to the applications, many efficient strategies have been proposed to improve the performance of QPSO. For example, Liu et al. (2005) introduced the mutation operation into QPSO to improve the search ability of the algorithm. Wang and Zhou (2007) proposed a local QPSO (LQPSO) as a generalized local search operator and incorporated LQPSO into a main QPSO algorithm, which leads to a hybrid QPSO scheme QPSO-LQPSO, with enhanced searching qualities. In Coelho (2008), it was shown that the chaotic mutation operation could diversify the population of QPSO and thus improve the performance of the algorithm. Pant et al. developed a new variant of QPSO, which used an interpolation-based recombination operator for generating a new solution vector in the search space (Pant et al., 2008). They also proposed a new mutation operator called the Sobal mutation to improve the performance of the QPSO algorithm (Pant et al., 2009). Xi et al. (2008) accelerated the convergence speed of QPSO by introducing a weighted mean best position into the algorithm. Huang et al. (2009) proposed an improved QPSO, employing a selection operation on the particles in order to filter the particle swarm and accelerate its convergence.

While empirical evidence has shown that the QPSO algorithm works well, thus far there has been little insight into how it works. In this paper, we make a comprehensive analysis of individual particle behavior for QPSO from the perspective of probability measure and show how to select the contraction-expansion (CE) coefficient, which is the most important algorithmic parameter. To achieve this goal, the work presented in the paper is mainly divided into the following two parts.

The first part includes the theoretical analysis of the algorithm. Since the behavior of the individual particle in QPSO exerts great influence on the convergence of the particle swarm and in turn on the convergence of the algorithm, we analyze the individual particle behavior for two different versions of QPSO. In QPSO, the CE coefficient controls the behavior of the individual particle, just as the inertia weight and acceleration coefficients influence the behavior of the particle in PSO. Therefore, this part of the work comes down to deriving the upper bound of the value of the CE coefficient that guarantees the convergence or boundedness of the particle.

The upper bound of the CE coefficient only provides the condition that leads to the convergence of the particle swarm, not the condition that ensures the efficiency of QPSO in practical applications. Thus, the second part of the paper involves empirical studies on how to select the value of the CE coefficient within the upper bound to lead the QPSO algorithm to good performance in general. With regard to this, we firstly perform stochastic simulations in order to visualize the influence of the value of the CE coefficient on particle convergence speed or the expected range of particle random vibration, as well as to verify the derived theoretical results of the parameter's upper bound. Secondly, we execute the QPSO algorithm with two parameter control methods on some well-known benchmark functions to find the parameters’ values resulting in good solutions in general. Finally, the QPSO algorithms with the parameters’ values found for both parameter control methods are tested on a suite of benchmark functions proposed by Suganthan at CEC 2005 (Suganthan et al., 2005), and the performances are compared with those of other PSO variants.

The rest of the paper is organized as follows. In Section 2, the basic principles of QPSO are introduced. Section 3 presents a theoretical analysis of individual particle behavior. Section 4 provides the experimental analysis of the algorithm. Some concluding remarks are given in the last section.

2.  Quantum-Behaved Particle Swarm Optimization

In the PSO with M individuals, each individual is treated as a volumeless particle in an N-dimensional space, with the current position vector and the velocity vector of particle i () at the nth iteration represented as and , respectively. The particle moves according to the following equations:
formula
1
formula
2
for , where c1 and c2 are known as the acceleration coefficients. Vector is the best previous position (the position giving the best objective function value or fitness value) of particle i, called the personal best (pbest) position, and vector is the position of the best particle among all the particles in the population, called the global best (gbest) position. Without loss of generality, we consider the following minimization problem:
formula
3
where f(X) is an objective function continuous almost everywhere and S is the feasible space. Accordingly, Pi,n can be updated by
formula
4
and Gn can be found by Gn=Pg,n, where . The parameters rji,n and Rji,n, varying with n for each i and j, are two different random numbers distributed uniformly on (0, 1), which is denoted by . Generally, the value of Vji,n is restricted in the interval [−Vmax, Vmax ].
Trajectory analysis (Clerc and Kennedy, 2002) demonstrated the fact that convergence of the PSO algorithm may be achieved if each particle converges to its local attractor, , defined at the coordinates
formula
5
or
formula
6
for , where , with regard to the random numbers rji,n and Rji,n in Equations (1) and (5). In PSO, the acceleration coefficients c1 and c2 are generally set to be equal, that is, c1=c2, and thus is a sequence of random numbers uniformly distributed on (0, 1). As a result, Equation (6) can be restated as
formula
7

The above equation indicates that Pi,n, the stochastic attractor of particle i, lies in the hyper-rectangle with Pi,n and Gn being the two ends of its diagonal so that it moves following Pi,n and Gn. In fact, as the particles are converging to their own local attractors, their current position, personal best positions, local attractors, and the global best positions are all converging to one point, leading the PSO algorithm to convergence. From the point view of Newtonian dynamics, in the process of convergence, the particle moves around and careens toward Pi,n with its kinetic energy (or velocity) declining to zero, like a returning satellite orbiting the earth. As such, the particle in PSO can be considered as the one flying in an attraction potential field centered at Pi,n in the Newtonian space. It has to be in a bound state for the sake of avoiding explosion and guaranteeing convergence. If these conditions are generalized to the case in which the particle in PSO has quantum behavior moving in an N-dimensional Hilbert space, it is also indispensable that the particle moves in a quantum potential field to ensure the bound state. From the perspective of quantum mechanics, the bound state in the quantum space, however, is entirely different from that in the Newtonian space, which may lead to a very different form of PSO. This is the motivation of the proposed QPSO algorithm which is described below (Sun, Feng et al., 2004).

In the quantum time-space framework, the quantum state of a particle is described by the wave function (Cohen-Tannoudji et al., 1997). In a three-dimensional space, the wave function of a particle satisfies the relation
formula
8
where Qdxdydz is the probability that the particle will appear in the infinitesimal element about the point (x, y, z). In other words, represents the probability density function satisfying
formula
9
Equation (8) or (9) gives the statistical interpretation for the wave function. Generally, varies in time according to the following equation
formula
10
where is Planck's constant and is the Hamiltonian operator defined by
formula
11
for a single particle of mass m in a potential field V(X). Equation (10) is known as the time-dependent Schrödinger equation.
We assume that each single particle in QPSO is treated as a spin-less particle moving in an N-dimensional Hilbert space with a given energy, and thus its state is characterized by a wave function which only depends on the position of the particle. Inspired by the convergence analysis of the particle in PSO (Clerc and Kennedy, 2002), we further assume that, at the nth iteration, particle i flies in the N-dimensional Hilbert space with a potential well centered at pji,n on the jth dimension (). In order to facilitate the description, we consider a particle in an one-dimensional space firstly and denote the position of the particle as X and Pi,n as p. With point p being the center of the potential well, the potential energy of the particle in the one-dimensional potential well is represented as
formula
12
where Y=Xp and is the intensity of the potential well. Hence, for this bound state problem, the Hamiltonian operator is
formula
13
The particle's state is subject to the following stationary Schrödinger equation:
formula
14
where E is the energy of the particle and is the wave function of the particle which only depends on its position. As presented in Theorem 1 below, we can obtain the normalized wave function of the particle in the bound state by solving the above equation with the bound condition: , as .
Theorem 1: 
For the particle that moves in the one-dimensional potential well formulated by Equation (12), its normalized wave function in the bound state is given by
formula
15
where .
Proof: 
Integrating Equation (14) with respect to Y from to and taking leads to
formula
16
For , Equation (14) can be written as
formula
17
where (E<0). To satisfy the bound condition
formula
18
the solution of Equation (17) must be
formula
19
Note that it can be proved that only an even wave function satisfies the bound condition in Equation (18), and then the solution of Equation (17) can be written as
formula
20
where A is the normalization constant. According to Equation (16), we have
formula
Thus
formula
21
and
formula
22
The function satisfies normalization condition
formula
23
leading to . is the characteristic length of the potential well. Inserting and into Equation (20), the normalized wave function can then be written as
formula
24

This completes the proof of the theorem.

In terms of the statistical interpretation of the wave function, the probability density function of Y is given by
formula
25
and the corresponding probability distribution function is
formula
26

With the given probability distribution function, we can measure the position of the particle using Monte Carlo inverse transformation. This technique is described in the proof of Theorem 2. Such a process of measuring the particle's position in quantum mechanics is essentially achieved by collapsing the quantum state to the classical state.

Theorem 2: 
If a particle moves in the bound state in the one-dimensional potential well as described by Equation (12), its position can be determined by using the stochastic equation
formula
27
where u is a random number uniformly distributed on (0, 1), that is, .
Proof: 
Let v be a random number uniformly distributed on (0, 1), that is
formula
Substituting v for F(Y) in Equation (26) and following Monte Carlo inverse transformation, we obtain
formula
28
Since , putting u=1−v leads to . Thus Equation (28) can be written as
formula
29
from which we can immediately obtain
formula
Because Y=Xp, we have
formula
30

This completes the proof of the theorem.

Now we generalize Equation (27) to the case in the N-dimensional Hilbert space where each dimension of the particle's position is bounded in a potential well and updated independently. Considering that the particle's position, its local attractor, the characteristic length of the potential well and the random variable u develop with the iteration number n, we can use the following equation to measure the jth () component of the position of particle i () at the (n+1)th iteration.
formula
31
where uji,n+1 is a sequence of random numbers uniformly distributed on (0, 1), varying with n for each i and j.
The value of Lji,n in Equation (31) can be determined by either of the following two equations (Sun, Feng et al., 2004; Sun, Xu et al., 2004):
formula
32
or
formula
33
where is known as the mean best (mbest) position which is defined by the average of the pbest positions of all particles, that is, . Therefore, the position of the particle is updated by using either of the following two equations:
formula
34
or
formula
35
The parameter in Equations (32)– (35) is a positive real number, the CE coefficient, which can be adjusted to balance the local and global search of the algorithm during the search process. The PSO with Equation (34) or (35) is the QPSO. To distinguish them, we denote the QPSO with Equation (34) as QPSO-Type 1 and that with Equation (35) as QPSO-Type 2.

The QPSO algorithm starts with the initialization of the particles’ current positions and their pbest positions (setting Pi,0=Xi,0), followed by the iterative update of the particle swarm. At each iteration of the procedure, the mbest position of the particle swarm is computed (for QPSO-Type 2) and the current position of each particle is updated according to Equation (34) or (35) with the coordinates of its local attractor evaluated by Equation (7). Before each particle updates its current position, its fitness value is evaluated together with an update of its pbest position and the current gbest position. In Equation (34) or (35), the probability of using either the + operation or the − operation is equal to 0.5. The iterative process continues until the termination condition is met.

The procedure of the QPSO algorithm is outlined in Algorithm 1. Note that , is used to denote the random numbers that are separately generated and uniformly distributed on (0, 1).

3.  Theoretical Analysis

The weighting of the CE coefficient in the QPSO algorithm may result in a kind of explosion as the position coordinates careen toward infinity. This section demonstrates that properly selected can prevent explosion, and further, this coefficient can induce the particle to converge to its local attractor (in QPSO-Type 1) or to be probabilistic bounded (in QPSO-Type 2).

As in the PSO algorithm, an important source of the swarm's search capability is the interactions among particles as they react to one another's findings. However, the theoretical analysis in this paper does not involve the analysis of interparticle effects but is focused on the stochastic movements of single particles.

3.1.  Preliminaries

This section provides a brief introduction to mathematical preliminaries on probability measure and sequences of random variables, which are essential to the theoretical analysis of a single particle's behavior. To study this knowledge in depth, one may refer to the related references such as Dudley (2003) and Shiryayev (1984).
formula
Definition 1: 

The set containing all possible outcomes of a random experiment is called the space of elementary events or the sample space. Each outcome or a point in is known as an elementary event or a sample point.

Definition 2: 

The space together with a -algebra of its subsets is a measurable space which is denoted by .

Definition 3: 

An ordered triple where

  • is a set of points ,

  • is a -algebra of subsets of ,

  • P is a probability measure on ,

is called a probabilistic model or a probability space. Here is the sample space or the space of elementary events, each set A in is known as an event, and P(A) is the probability of the event A.

Let be the real line. ( or is a measure space. or is called the Borel algebra of subsets.

Definition 4: 

A real function defined on is an -measurable function or a random variable, if for every ; or equivalently, if the inverse image is a measurable set in .

Definition 5: 

Let be random variables defined on a probability space . The sequence is called a sequence of independent identically distributed (i.i.d.) random variables, if each random variable has the same probability distribution as the others and all are mutually independent.

Just as in analysis, there are various kinds of convergence of random variables in probability theory. Four of these are particularly important: with probability one, in probability, in distribution, and in mean of order r (Shiryayev, 1984).

Definition 6: 
Let and be random variables defined on a probability space . The sequence converges almost surely (with probability one, almost everywhere) to the random variable if
formula
36
that is, if the set of sample points for which does not converge to has probability zero. This convergence is denoted by (P-a.s.), or or .
Definition 7: 
Let and be random variables defined on a probability space . The sequence converges in probability to the random variable (notation: ) if for every ,
formula
37
or
formula
38
Definition 8: 
Let and be random variables defined on a probability space with distribution function Fn(x) and f(X) respectively. The sequence converges to a distribution to the random variable (notation: ) if for every continuous point x,
formula
39
is satisfied.
Definition 9: 
Let and be random variables defined on a probability space with , . The sequence converges to a mean of order r to the random variable (notation: ) if
formula
40

The following theorem indicates the relationships between the four kinds of convergence. It should be noted that the converses of Equations (41), (42), and (43) are false in general.

Theorem 3: 
We have the following implications:
formula
41
formula
42
formula
43

The following two theorems present the strong and weak laws of large numbers, respectively.

Theorem 4 (Kolmogorov's Strong Law of Large Numbers): 
Let be a sequence of independent identically distributed random variables with . Then
formula
44
Theorem 5 (Khintchine's Weak Law of Large Numbers): 
Let be a sequence of independent identically distributed random variables with . Then for every ,
formula
45
or
formula
46

3.2.  Simplification of the Iterative Equations

It is obvious that the convergence or boundedness of the position of an individual particle is consistent with the convergence or boundedness of its component in each dimension. In detail, for QPSO-Type 1, the necessary and sufficient condition for Xi,n to converge to Pi,n in any kind of convergence is that Xji,n converges to pji,n for each in that kind of convergence; for QPSO-Type 2, the necessary and sufficient condition for Xi,n to be probabilistic bounded is that Xji,n is probabilistic bounded for each .

Furthermore, it appears from Equation (34) or (35) that each dimension of the particle's position is updated independently following the same equation. The only link between the dimensions of the problem space relies on the objective function, and in turn, through the locations of the personal and global best positions found so far, or the mean best position among particles. As a result, to investigate the convergence or boundedness of an individual particle, we only need to analyze the convergence or boundedness of the component in any dimension. Without loss of universality, the issue of the convergence or boundedness for a particle in the N-dimensional space can be reduced to the one for a single particle in the one-dimensional space, and we only need to investigate the behavior of the particle in the one-dimensional space using the iterative equation given by
formula
47
or
formula
48

In the above equations, the local attractor of the particle and the mean best position are replaced by p and C, which are treated as probabilistic bounded random variables, instead of constants as in Clerc and Kennedy (2002). Here, the probabilistic boundedness of p and C means that and . The position sequence is a sequence of random variables and is a sequence of independent random variables with for all n>0. For convenience, we denote a particle moving according to Equation (34) or (47) as a Type-1 particle, and a particle moving according to Equation (35) or (48) as a Type-2 particle. The remainder of this section focuses on the behaviors of the two types of particles.

Moreover, it is obvious that the sequence in Equation (47) or (48) is a sequence of independent identically distributed random variables, since each random variable ui is generated independently and has the same probability distribution, that is, the uniform distribution on (0, 1), as the others.

3.3.  Convergence of the Type-1 Particle

Rewriting Equation (47) as
formula
49
where , we obtain
formula
50
where X0 is the initial position of the particle. As such, the convergence analysis of the sequence of random variables can be reduced to that of the infinite product
formula
51

The following theorem gives an integral which is of significance to the succedent analysis of the convergence of and .

Theorem 6: 
The improper integral
formula
52
holds with the Euler-Mascheroni constant .
Proof: 
Let , and thus we have
formula
Since
formula
where is the Gamma function,
formula
From Courant (1989) we have that
formula
which implies that
formula

This completes the proof of the theorem.

With the above preliminaries, in the following part of this section, we derive the sufficient and necessary condition for to converge to random variable p in each kind of convergence. The detailed proof of the almost sure convergence is provided in the text, while those of other kinds of convergence are presented in the Appendix.

3.3.1.  Almost Sure Convergence

Lemma 1: 
If is a sequence of independent identically distributed random variables with for all n>0 and , then
formula
53
Proof: 
Since is a sequence of independent identically distributed (i.i.d.) random variables, is also a sequence of i.i.d. random variables. Theorem 6 implies that
formula
Thus, by Theorem 4 (Kolmogorov's strong law of large numbers), we have
formula
This completes the proof of the lemma.
Theorem 7: 

The necessary and sufficient condition for the position sequence of a Type-1 particle to converge almost surely to p is that .

Proof: 

Before the proof, we provide the following two groups of equivalent propositions at first.

(a) If , then
formula
which is equivalent to
formula
54
where Z+ is the set of all positive integers.
(b) Lemma 1 has the following equivalences
formula
Proof of Necessity:

If , the propositions in (a) hold. Since the propositions in (b) also hold due to Lemma 1, for any positive integer m, there exists K=max(K1, K2) such that whenever ,

formula
and
formula
55
holds simultaneously, and thus it follows that . That is, , , such that whenever ,
formula
This proposition is equivalent to
formula
This ends the proof of necessity.

Proof of Sufficiency

(i) From the equivalences of Lemma 1, we have that , such that whenever
formula
56
Since , , we have
formula
57
and therefore
formula
From Equation (56), we have
formula
and thus
formula
58
Since , we find that
formula
resulting in the fact that
formula
59
(ii) , such that whenever , , from which we have
formula
and thus have
formula
60
From Equations (59) and (60), we have that , such that whenever ,
formula
which is equivalent to
formula
61
Since , we obtain
formula
The condition that implies that , so
formula
62
From the inequality in Equation (61) and from Equation (62), we obtain
formula
and thus
formula
Considering Equation (54), we find that
formula
63

This ends the proof of sufficiency.

This completes the proof of the theorem.

3.3.2.  Convergence in Probability

Theorem 8: 

The necessary and sufficient condition for the position sequence of a Type-1 particle to converge to p in probability is that .

Proof: 

See the proof in the Appendix.

3.3.3.  Convergence in Distribution

Theorem 9: 

The necessary and sufficient condition for the position sequence of a Type-1 particle to converge to p in distribution is that .

Proof: 

See the proof in the Appendix.

3.3.4.  Convergence in Mean of Order r

Theorem 10: 

The necessary and sufficient condition for the position sequence of a Type-1 particle to converge to p in mean of order r is that .

Proof: 

See the proof in the Appendix.

3.3.5.  Some Discussion

The above analysis indicates that the sufficient and necessary condition for the position sequence of a Type-1 particle to converge to p in any of four kinds of convergence is that . In other words, the four kinds of convergence in this case are equivalent conditionally. This conclusion does not conflict with Theorem 3. It is shown by Equation (41) that almost sure convergence implies convergence in probability, while the converse is in general not true. Under some special circumstances, the converse of Equation (41) may hold. Kolmogorov's strong law and Khintchine's weak law of large numbers reveal that a convergence in probability may imply almost sure convergence for certain cases, leading to the equivalence between the two kinds of convergence of the Type-1 particle.

In any case, it is due to the probabilistic boundedness of p and its uniform distribution on (P, G) or (G, P) that the converse of Equation (42) holds for the Type-1 particle, as shown in the proof of Theorem 9. Owing to the equivalence between the convergence in probability and the convergence in distribution, the point p can be treated as a constant as we do in the stochastic simulations in the next section. Furthermore, as shown by the proof of Theorem 10 in the Appendix, since is uniformly integrable, the sufficient and necessary condition for its convergence in the mean of order r is the same as for its convergence in probability.

It is shown in the proof of Theorem 9 that when , the position of the Type-1 particle can be any real number but infinity, which means that the position is probabilistically bounded as n tends to infinity. Therefore, we may conclude that if , the position of the particle is probabilistically bounded; otherwise, it is divergent.

3.4.  Boundedness of the Type-2 Particle

Now we turn our attention to the behavior of a Type-2 particle, beginning by rewriting Equation (48) as
formula
64
where and . From the above equation, we have the following two inequalities:
formula
65
and
formula
66

Based on the inequalities in Equations (65) and (66), the proofs of the following theorems derive the necessary and sufficient condition for the probabilistic boundedness of the Type-2 particle's position, which is also related to the probabilistic boundedness of .

Theorem 11: 

The necessary and sufficient condition for to be probabilistically bounded (i.e., ) is that .

Proof: 

The proof is essentially given by the first two parts of Theorem 10 in the Appendix.

Theorem 12: 

The necessary and sufficient condition for the position sequence of a Type-2 particle to be probabilistically bounded (i.e., ), is that .

Proof: 
Since is a continuous random variable, it is evident that . Therefore, implies that
formula
Denoting where , we have that for every n>0, , namely
formula
67
Proof of Sufficiency: By replacing |Cp| in Equation (65) by that in the inequality in Equation (67), we obtain the following inequality:
formula
from which we find that
formula
Thus the following inequality holds:
formula
68
Since ,
formula
By Theorem 11, we have that, whenever ,
formula
Considering that , we have . Therefore
formula
from which we immediately have that , namely, |XnC| is probabilistically bounded. This implies that Xn is also probabilistically bounded, that is, . This ends the proof of sufficiency.
Proof of Necessity: Replacing |Cp| in Equation (66) by that in Equation (67), we find that
formula
69
from which we obtain
formula
70
Thus we have that
formula
As a result, it follows that
formula
71
If Xn is probabilistically bounded, |XnC| is also probabilistically bounded, that is, . Since , |XnC|+r is probabilistically bounded, that is,
formula
Considering the inequality in Equation (71), we have
formula
which immediately results in
formula
72
Due to the probabilistic boundedness of |X0C|+r, we have
formula
73
From Equations (72) and (73), we obtain which means as shown by Theorem 11. This implies that is the necessary condition for the probabilistic boundedness of Xn. This ends the proof of necessity.

This completes the proof of the theorem.

Theorem 12 reveals that the sufficient and necessary condition for the probabilistic boundedness of a Type 2 particle is the same as that for the probabilistic boundedness of a Type 1 particle, since the behaviors of both types of particles are related to the probabilistic boundedness of . Besides, the behavior of the Type 2 particle is also influenced by point C. In practice, when the QPSO algorithm is running, the personal best positions of all the particles converge to the same point. This implies that |Cp| converges to zero almost surely, and thus also converges to zero almost surely. Hence, if and only if , . According to the inequality in Equation (68), we can find that or .

4.  Experimental Analysis

4.1.  Stochastic Simulations on the Particle's Behavior

It has been demonstrated that the behavior of an individual particle in QPSO is subject to how the value of is selected. In order to verify the correctness of the previous analysis, we carried out stochastic simulations for the behavior of the particle. Two groups of simulations were implemented, one for a Type-1 particle and the other for a Type-2 particle.

The value of for either group of simulations was selected from a series of numbers in [0.5, 2.0] and the corresponding maximum number of iterations (nmax) executed was sufficiently large for the purpose of convergence, boundedness, and divergence. When the stochastic simulation was executed, the logarithmic value of the distance between the current position Xn and the point p was recorded as the ordinate, and the number of iterations was the abscissa. Without loss of generality, in the simulations for a Type-1 particle, p was fixed at the origin, and the initial position of the particle was set as X0=1000; in the simulations for a Type-2 particle, C was fixed at X=0.001, and p and the initial position were the same as those in the simulations for the Type 1 particle. The simulations were performed on Matlab 7.0 and the results are shown in Figures 1 and 2.

Figure 1:

Simulation results for a Type-1 particle. When , (i.e., ) as n increased, and the smaller resulted in a faster convergence speed. When , (i.e., ) as n increased, and the larger led to a faster divergence speed.

Figure 1:

Simulation results for a Type-1 particle. When , (i.e., ) as n increased, and the smaller resulted in a faster convergence speed. When , (i.e., ) as n increased, and the larger led to a faster divergence speed.

Figure 2:

Simulation results for a Type-2 particle. When , ln |Xnp| was bounded as n increased, and the smaller resulted in a narrower oscillation range. When , (i.e., ) as n increased, and the larger led to a faster divergence speed.

Figure 2:

Simulation results for a Type-2 particle. When , ln |Xnp| was bounded as n increased, and the smaller resulted in a narrower oscillation range. When , (i.e., ) as n increased, and the larger led to a faster divergence speed.

Figure 1 shows the simulation results for the position of the Type-1 particle. When ln |Xnp| is smaller than −700 or is larger than 700, |Xnp| reaches the minimum or maximum positive value that the computer can identify, and thus we can consider that |Xnp| converges to zero (Xn converges to p) or diverges to infinity. It can be seen from Figure 1 that the particle's position Xn converged to p when , and diverged when . When Xn converged to p, the smaller the value of , the faster Xn converged. On the other hand, when Xn diverged, the larger value of resulted in a faster speed at which the position careened to infinity.

As for a Type-2 particle, according to Figure 2, its position oscillated around p and C when , while the position exploded when . The results also show that when Xn was probabilistically bounded, the smaller value of resulted in a narrower oscillation range of the particle's position. When Xn diverged, the larger value of led to a faster divergence speed.

It can be concluded that there exists an within (1.775, 1.785) such that whenever , the position of the Type-1 particle converges (when ) or is probabilistic bounded (when ), and that of the Type-2 particle is probabilistic bounded; otherwise the position of either type of particles diverges. Therefore, it is obvious that the simulation results are consistent with the theoretical results in the previous section that setting can prevent explosion of the particle's position.

The fact that the smaller value of results in the faster convergence speed of a Type-1 particle and the narrower oscillation range of a Type-2 particle can be illustrated as follows. For a Type-1 particle, according to Equation (47), we define the convergence speed at the nth iteration by
formula
74
Thus we have
formula
75
and note that the velocity at which |Xnp| declines, increases with the smaller cn, namely, the smaller . For a Type-2 particle, according to Equation (48), we find that
formula
76
implying that a smaller value of results in a smaller conditional expected value of |Xnp| and thus a narrower oscillation range of Xn. It is obvious that the slower convergence of a Type-1 particle or the wider oscillation range of a Type-2 particle means a stronger global search ability of the particle, while the faster convergence or the narrower oscillating scope implies a stronger local search ability. Excessive global search may result in slow convergence of the QPSO algorithm, and on the other hand, excessive local search may cause premature convergence of the algorithm. It is significant, therefore, to balance the global search (exploration) and local search (exploitation) to obtain a generally good performance of the algorithm. The parameter , whose value selection is studied empirically in the rest of this section, plays a major role in balancing the exploration and exploitation of the particles in QPSO.

4.2.  Methods of Parameter Control

When an algorithm is employed to solve a problem at hand, one of the most important issues is how to select its parameters. It is evident that is the most influential parameter on the convergence properties of QPSO, except for the population size. Although it was demonstrated in Section 3 that it is necessary and sufficient to set to prevent the individual from explosion and guarantee the convergence of the particle swarm, this does not mean that any value of less than or equal to can lead to a satisfactory performance of QPSO in practical applications. In the remainder of this section, we aim at finding out through empirical study how to control and select the value of so that QPSO may yield good performance in general. We also provide a performance comparison with other forms of PSO. To do so, we tested the algorithms on 12 benchmark functions F1 to F12 presented by Suganthan et al. (2005), which are not listed here due to space limitations.

When QPSO is applied to practical problems, there are several control methods for the parameter . A simple approach is to set to be a fixed value when the algorithm is running. Another efficient method is to decrease the value of linearly in the course of search. That is, is determined by
formula
where and () are the initial and final values of , respectively, n is the current iteration number, and nmax is the maximum number of allowable iterations.

We also propose a variant of QPSO-Type 2 in which the mean best position Cn in Equation (15) is replaced by the pbest position of a randomly selected particle in the swarm at each iteration. To distinguish between them, we denote the QPSO-Type 2 with Cn as QPSO-Type 2-I and that with a randomly selected pbest position as QPSO-Type 2-II.

4.3.  Empirical Studies on Parameter Selection

When either control strategy for is used in QPSO, the value of in the fixed-value method or and in the linear time-varying approach must be determined. To select the values that can yield generally good algorithmic performance, we tested QPSO-Type 1, QPSO-Type 2-I and QPSO-Type 2-II on three frequently used functions in Suganthan's benchmark suite: Shifted Rosenbrock Function (F6), Shifted Rotated Griewank's Function (F7), and Shifted Rastrigin's Function (F9), using the two methods of controlling . Rastrigin and Griewank are two difficult multimodal problems and Rosenbrock is a unimodal problem. For each of these functions, N=30. The expressions and bounds of the functions are shown in Suganthan et al. (2005). For each setting, each algorithm using 20 particles was tested for 100 runs on each problem. In each trial run, the initial positions of the particles were determined randomly within the search bounds. To determine the effectiveness of each algorithm for the setting of each control method on each problem, the best objective function value (i.e., the best fitness value) found after 3,000 iterations was averaged over 100 runs of the algorithm for the same setting and the same benchmark function. The results were compared by settings across the three benchmarks. The best setting of each control method was selected by ranking the averaged best objective function values for each problem, adding the ranks, and taking the value that had the lowest summed rank, provided that the performance is acceptable (in the top half of the rankings) in all the tests for a particular setting of .

The results for QPSO-Type 1 are presented in Table 1. When the fixed-value method was used, was set to different values smaller than in each case. The results obtained for outside the range [1.2, 0.85] are very poor and are not listed in the table. It can be observed from Table 1 that if the value of was fixed during the search process, QPSO-Type 1 generated quite different objective function values for different values of , particularly for the Shifted Rosenbrock Function, although the difference between two adjacent values of is not significant. The results indicate that the performance of QPSO-Type 1 was somewhat sensitive to the value of . The best results were obtained when . When time-varying was used, and () were selected from a series of different values less than . Only acceptable results are listed in the table. It can be seen that the algorithmic performance was sensitive to both and . It is also found that decreasing linearly from 1.0 to 0.9 led to the best performance in general.

Table 1:
Mean best fitness values obtained by QPSO-Type 1 on three benchmark functions.
Fixed
RosenbrockRastriginGriewankRosenbrockRastriginGriewank
1.2 1.9761e+03 73.2189 75.2146 0.97 140.7474 65.2704 0.0216 
1.10 207.2907 51.3072 0.7604 0.96 128.861 63.8414 0.0229 
1.05 179.4073 48.2471 0.0676 0.95 173.8644 69.3123 0.0211 
1.02 178.5613 53.0128 0.0259 0.94 157.0963 71.4457 0.0254 
1.01 151.0301 55.9349 0.0246 0.93 168.8516 70.5280 0.0208 
1.00 138.0746 56.4232 0.0218 0.92 193.6964 75.5213 0.0281 
0.99 188.3713 52.7332 0.0217 0.90 152.7539 84.4062 0.1007 
0.98 170.1586 60.7105 0.0196 0.85 282.2530 96.1226 51.0799 
Linearly-Decreasing        
 Rosenbrock Rastrigin Griewank  Rosenbrock Rastrigin Griewank 
 195.5114 60.3980 0.0688  174.0003 59.4232 0.0390 
 180.0561 64.5759 0.2481  200.9611 65.3739 0.2538 
 178.0053 63.8211 0.4468  204.2796 61.2373 0.6978 
 170.1478 63.3740 0.7066  261.7130 71.5227 0.9640 
 143.5103 50.7536 0.0276  208.1228 81.9640 2.1316 
 181.5607 57.9190 0.0346  262.8041 81.9872 7.3829 
 145.4834 55.5282 0.1123  311.9669 87.7742 12.0388 
 165.9128 61.5130 0.3720  506.0865 87.1547 21.0660 
 218.6484 58.8480 0.5683  7.6355e+06 119.1037 259.0206 
 139.9815 54.4278 0.0209  3.3683e+06 111.3555 282.3997 
Fixed
RosenbrockRastriginGriewankRosenbrockRastriginGriewank
1.2 1.9761e+03 73.2189 75.2146 0.97 140.7474 65.2704 0.0216 
1.10 207.2907 51.3072 0.7604 0.96 128.861 63.8414 0.0229 
1.05 179.4073 48.2471 0.0676 0.95 173.8644 69.3123 0.0211 
1.02 178.5613 53.0128 0.0259 0.94 157.0963 71.4457 0.0254 
1.01 151.0301 55.9349 0.0246 0.93 168.8516 70.5280 0.0208 
1.00 138.0746 56.4232 0.0218 0.92 193.6964 75.5213 0.0281 
0.99 188.3713 52.7332 0.0217 0.90 152.7539 84.4062 0.1007 
0.98 170.1586 60.7105 0.0196 0.85 282.2530 96.1226 51.0799 
Linearly-Decreasing        
 Rosenbrock Rastrigin Griewank  Rosenbrock Rastrigin Griewank 
 195.5114 60.3980 0.0688  174.0003 59.4232 0.0390 
 180.0561 64.5759 0.2481  200.9611 65.3739 0.2538 
 178.0053 63.8211 0.4468  204.2796 61.2373 0.6978 
 170.1478 63.3740 0.7066  261.7130 71.5227 0.9640 
 143.5103 50.7536 0.0276  208.1228 81.9640 2.1316 
 181.5607 57.9190 0.0346  262.8041 81.9872 7.3829 
 145.4834 55.5282 0.1123  311.9669 87.7742 12.0388 
 165.9128 61.5130 0.3720  506.0865 87.1547 21.0660 
 218.6484 58.8480 0.5683  7.6355e+06 119.1037 259.0206 
 139.9815 54.4278 0.0209  3.3683e+06 111.3555 282.3997 

Table 2 records the results for QPSO-Type 2-I. For the case in which the fixed-value method was used, the results of the algorithm for outside the range [0.6, 1.2] are not listed because of poor quality. It can be observed from the results that the performance of QPSO-Type 2-I was less sensitive to than that of QPSO-Type 1. The results seemed fairly stable when was valued in the interval [0.8, 0.7]. The value of that generated the best results in this case was 0.75. When decreased linearly in the course of search, the objective function values obtained by QPSO-Type 2-I with different settings were not so different from each other as those obtained by QPSO-Type 1. It has been identified that varying linearly from 1.0 to 0.5 could yield the best quality results.

Table 2:
Mean best fitness values obtained by QPSO-Type 2-I on three benchmark functions.
Fixed
RosenbrockRastriginGriewankRosenbrockRastriginGriewank
1.20 5.0457e+007 209.4031 35.3916 0.75 82.9908 39.0991 0.0203 
1.00 2.9384e+004 164.7224 1.5544 0.74 119.3931 41.9011 0.0205 
0.95 1.5836e+003 150.1191 0.9811 0.73 167.4110 45.8569 0.0196 
0.90 153.7730 127.1872 0.1448 0.72 162.2084 49.3278 0.0251 
0.85 129.7591 97.7650 0.0193 0.71 275.5061 49.9125 0.0745 
0.80 115.6558 50.7660 0.0216 0.70 217.7180 58.0127 1.8027 
0.78 138.0215 39.3672 0.0184 0.65 1.0587e+005 73.1374 322.2796 
0.76 108.5730 35.5272 0.0228 0.60 1.8762e+007 102.6281 877.2902 
Linearly-Decreasing        
 Rosenbrock Rastrigin Griewank  Rosenbrock Rastrigin Griewank 
 199.7309 39.4630 0.0243  113.5787 30.6210 0.0391 
 157.9285 33.7100 0.0422  138.9362 31.6721 0.0667 
 192.3246 31.8872 0.0709  106.3916 30.2193 0.0186 
 159.6439 34.1352 0.1024  124.7271 31.0468 0.0246 
 210.3639 35.8151 0.0233  149.4981 31.4059 0.0442 
 165.3690 30.8869 0.0343  177.1897 35.1173 0.0885 
 157.2430 32.1224 0.0491  179.4460 32.9276 0.0226 
 184.4184 33.5440 0.0697  169.4544 36.3398 0.0471 
 74.5490 33.7418 0.0190  206.0009 38.5983 0.1958 
 88.0494 29.9218 0.0208  400.7041 37.8844 0.5162 
Fixed
RosenbrockRastriginGriewankRosenbrockRastriginGriewank
1.20 5.0457e+007 209.4031 35.3916 0.75 82.9908 39.0991 0.0203 
1.00 2.9384e+004 164.7224 1.5544 0.74 119.3931 41.9011 0.0205 
0.95 1.5836e+003 150.1191 0.9811 0.73 167.4110 45.8569 0.0196 
0.90 153.7730 127.1872 0.1448 0.72 162.2084 49.3278 0.0251 
0.85 129.7591 97.7650 0.0193 0.71 275.5061 49.9125 0.0745 
0.80 115.6558 50.7660 0.0216 0.70 217.7180 58.0127 1.8027 
0.78 138.0215 39.3672 0.0184 0.65 1.0587e+005 73.1374 322.2796 
0.76 108.5730 35.5272 0.0228 0.60 1.8762e+007 102.6281 877.2902 
Linearly-Decreasing        
 Rosenbrock Rastrigin Griewank  Rosenbrock Rastrigin Griewank 
 199.7309 39.4630 0.0243  113.5787 30.6210 0.0391 
 157.9285 33.7100 0.0422  138.9362 31.6721 0.0667 
 192.3246 31.8872 0.0709  106.3916 30.2193 0.0186 
 159.6439 34.1352 0.1024  124.7271 31.0468 0.0246 
 210.3639 35.8151 0.0233  149.4981 31.4059 0.0442 
 165.3690 30.8869 0.0343  177.1897 35.1173 0.0885 
 157.2430 32.1224 0.0491  179.4460 32.9276 0.0226 
 184.4184 33.5440 0.0697  169.4544 36.3398 0.0471 
 74.5490 33.7418 0.0190  206.0009 38.5983 0.1958 
 88.0494 29.9218 0.0208  400.7041 37.8844 0.5162 

The results for QPSO-Type 2-II are summarized in Table 3. It is clear from the results that the value of , whether it was fixed or time-varying, should be set to be relatively smaller so that the algorithm was comparable in performance with the other types of QPSO. The results for outside the range [0.4, 0.8] were of poor quality and are not listed in the table. As shown in Table 3, were set in the range of 0.5 to 0.6 when the fixed-value method was used. The best results were obtained by setting . On the other hand, the algorithm exhibited the best performance when decreased linearly from 0.6 to 0.5 for the linear time-varying method.

Table 3:
Mean best fitness values obtained by QPSO-Type 2-II on three benchmark functions.
Fixed
RosenbrockRastriginGriewankRosenbrockRastriginGriewank
0.80 7.7628e+05 195.3315 2.8068 0.55 167.7840 45.7357 0.0193 
0.70 176.5589 165.4325 0.2060 0.54 105.7474 42.4817 0.0163 
0.65 93.9793 132.4897 0.0162 0.53 111.5862 47.1062 0.0198 
0.60 107.5629 71.4177 0.0206 0.52 182.5722 53.4841 0.0225 
0.59 120.8616 64.2002 0.0195 0.51 176.2885 51.7469 0.0296 
0.58 145.5502 47.8130 0.0196 0.50 173.9438 56.2267 0.0217 
0.57 145.3949 45.4865 0.0177 0.45 1.1344e+05 85.6050 169.9241 
0.56 126.0142 46.9918 0.0232 0.40 9.3422e+07 113.9682 1.2831e+03 
Linearly-Decreasing        
 Rosenbrock Rastrigin Griewank  Rosenbrock Rastrigin Griewank 
 270.1094 123.8036 0.0366  81.5586 44.3686 0.0234 
 205.8190 54.4520 0.0152  293.5975 47.4461 0.0202 
 285.8976 48.9965 0.0549  131.1352 69.5491 0.0166 
 554.4213 46.6750 0.0341  176.7990 42.3857 0.0193 
 185.8477 109.5310 0.0173  131.4735 45.0470 0.0164 
 206.1129 49.5677 0.0165  234.9833 43.7979 0.0343 
 142.1582 49.8265 0.0315  89.6543 43.8327 0.0150 
 264.8597 44.7424 0.0579  157.4971 41.0110 0.0208 
 187.7583 92.9885 0.0118  209.7387 42.6532 0.0190 
 127.4620 44.5099 0.0143  306.1035 46.4197 0.0394 
Fixed
RosenbrockRastriginGriewankRosenbrockRastriginGriewank
0.80 7.7628e+05 195.3315 2.8068 0.55 167.7840 45.7357 0.0193 
0.70 176.5589 165.4325 0.2060 0.54 105.7474 42.4817 0.0163 
0.65 93.9793 132.4897 0.0162 0.53 111.5862 47.1062 0.0198 
0.60 107.5629 71.4177 0.0206 0.52 182.5722 53.4841 0.0225 
0.59 120.8616 64.2002 0.0195 0.51 176.2885 51.7469 0.0296 
0.58 145.5502 47.8130 0.0196 0.50 173.9438 56.2267 0.0217 
0.57 145.3949 45.4865 0.0177 0.45 1.1344e+05 85.6050 169.9241 
0.56 126.0142 46.9918 0.0232 0.40 9.3422e+07 113.9682 1.2831e+03 
Linearly-Decreasing        
 Rosenbrock Rastrigin Griewank  Rosenbrock Rastrigin Griewank 
 270.1094 123.8036 0.0366  81.5586 44.3686 0.0234 
 205.8190 54.4520 0.0152  293.5975 47.4461 0.0202 
 285.8976 48.9965 0.0549  131.1352 69.5491 0.0166 
 554.4213 46.6750 0.0341  176.7990 42.3857 0.0193 
 185.8477 109.5310 0.0173  131.4735 45.0470 0.0164 
 206.1129 49.5677 0.0165  234.9833 43.7979 0.0343 
 142.1582 49.8265 0.0315  89.6543 43.8327 0.0150 
 264.8597 44.7424 0.0579  157.4971 41.0110 0.0208 
 187.7583 92.9885 0.0118  209.7387 42.6532 0.0190 
 127.4620 44.5099 0.0143  306.1035 46.4197 0.0394 

4.4.  Performance Comparison

To determine whether QPSO can be as effective as other variants of PSO, PSO with inertia weight (PSO-In; Shi and Eberhart, 1998a, 1998b, 1999), PSO with constriction factor (PSO-Co; Clerc, 1999; Clerc and Kennedy, 2002), the standard PSO (Bratton and Kennedy, 2007), Gaussian PSO (Secrest and Lamont, 2003), Gaussian bare-bones PSO (Gaussian BBPSO; Kennedy, 2003, 2004), PSO with the exponential distribution (PSO-E; Krohling and Coelho, 2006), Lévy PSO (Richer and Blackwell, 2006), and QPSO were all compared by running a series of experiments on the first 12 functions from the CEC2005 benchmark suite (Suganthan et al., 2005). Functions F1 to F6 are unimodal, while functions F7 to F12 are multimodal. Each algorithm was run 100 times on each problem using 20 particles to search the global best fitness value. In each run, the particles in the algorithm started in new and randomly-generated positions, which were uniformly distributed within the search bounds. Each run of each algorithm lasted 3,000 iterations and the best fitness value (objective function value) for each run was recorded.

For the three types of QPSO, both methods of controlling were used and the parameters for each case were set as those values that generated the best results in the previous experiments. The parameter configurations of other PSO variants were the same as those recommended by the existing publications. For PSO-In, Shi and Eberhart (1998b) showed that the PSO-In with linearly decreasing inertia weight performs better than the one with fixed inertia weight. They varied the inertia weight linearly from 0.9 to 0.4 over the course of the run and fixed the acceleration coefficients (c1 and c2) at 2.0 in their empirical study (Shi and Eberhart, 1999). For PSO-Co, Clerc and Kennedy found that the values of the constriction factor and acceleration coefficients (c1 and c2) need to satisfy some constraints in order for the particle's trajectory to converge without a restriction on velocity (Clerc, 1999; Clerc and Kennedy, 2002). They recommended using a value of 4.1 for the sum of c1 and c2, which results in a value of the constriction factor and c1=c2=2.05. Eberhart and Shi (2000) also used these values of the parameters when comparing the performance of PSO-Co with that of PSO-In. In our experiments, we also employed these parameter values for PSO-In and PSO-Co, although they may not be optimal. For the standard PSO, the LBEST ring topology was used with other parameters set as those in PSO-Co (Bratton and Kennedy, 2007). For PSO-E, except for the population size and the maximum number of iterations, all the other parameters were configured as those in Krohling and Coelho (2006). The parameter configurations for Gaussian PSO, Gaussian BBPSO, and Lévy PSO were the same as those in Secrest and Lamont (2003), in Kennedy (2003), and in Richer and Blackwell (2006), respectively.

The mean best fitness value and standard deviation out of 100 runs of each algorithm on each problem is presented in Table 4. To investigate whether the differences in mean best fitness values between algorithms were significant, the mean values for each problem were analyzed using a multiple comparison procedure, an ANOVA with 0.05 as the level of significance. Unlike Tukey's honestly significant (THS) difference test used in Richer and Blackwell (2006), the procedure employed in this work is called a stepdown procedure which takes into account that all but one of the comparisons are less extreme than the range. When doing all pairwise comparisons, this approach is the best available if confidence intervals are not needed and sample sizes are equal (Day and Quinn, 1989).

Table 4:
Mean (SD) of the best fitness values over 100 trial runs of different algorithms.
AlgorithmsF1F2F3F4F5F6
PSO-In 3.8773e-13 785.0932 3.9733e+7 1.1249e+4 6.0547e+3 263.7252 
 (1.6083e-12) (661.2154) (4.6433e+7) (5.4394e+3) (2.0346e+3) (437.4145) 
PSO-Co 1.5713e-26 0.1267 8.6472e+6 1.3219e+4 7.6892e+3 123.0243 
 (1.4427e-25) (0.3796) (9.1219e+6) (6.0874e+3) (2.3917e+3) (266.2520) 
Standard PSO 8.2929e-26 78.2831 6.6185e+6 1.3312e+4 6.2884e+3 153.5178 
 (1.2289e-25) (52.3272) (3.0124e+6) (4.1076e+3) (1.4318e+3) (246.1049) 
Gaussian PSO 7.3661e-26 0.0988 1.1669e+7 2.3982e+4 8.0279e+3 150.7872 
 (5.9181e-25) (0.3362) (2.5153e+7) (1.2512e+4) (2.3704e+3) (303.3368) 
Gaussian 1.7869e-25 16.8751 7.7940e+6 1.1405e+4 9.5814e+3 144.1377 
BBPSO (8.4585e-25) (16.2021) (4.3240e+6) (6.7712e+3) (3.0227e+3) (165.2616) 
PSO-E 5.2531e-24 20.2750 6.2852e+6 8.2706e+3 7.2562e+3 189.8292 
 (2.2395e-23) (15.2414) (2.8036e+6) (3.6254e+3) (1.8666e+3) (375.8636) 
Lévy PSO 1.1880e-24 36.9986 1.7366e+07 7.4842e+3 8.2543e+3 133.9526 
 (1.1455e-23) (29.1360) (1.9001e+7) (6.6588e+3) (2.2297e+3) (293.8460) 
QPSO-Type 1 3.5936e-28 40.2282 4.8847e+6 6.2397e+3 8.0749e+3 138.0746 
( 1.00) (1.5180e-28) (23.3222) (2.1489e+6) (2.4129e+3) (1.7099e+3) (209.1735) 
QPSO-Type 1 5.0866e-29 4.5003 3.2820e+6 6.4303e+3 7.8471e+3 139.9815 
((4.4076e-29) (2.9147) (1.9953e+6) (2.9744e+3) (1.7878e+3) (206.8138) 
QPSO-Type 2-I 1.9838e-27 0.1771 1.6559e+6 3.1321e+3 5.7853e+3 82.9908 
(0.75) (5.2716e-28) (0.1137) (7.1264e+5) (2.0222e+3) (1.2483e+3) (119.836) 
QPSO-Type 2-I 1.2672e-27 120.6051 4.4257e+6 4.0049e+3 3.3684e+3 88.0494 
((3.7147e-28) (62.2340) (2.3302e+6) (2.7218e+3) (975.6551) (159.7481) 
QPSO-Type 2-II 3.1554e-36 0.0715 1.8544e+6 3.1443e+3 5.7144e+3 105.7474 
( 0.54) (2.3913e-36) (0.0530) (6.4710e+5) (3.8785e+3) (1.4898e+3) (155.4583) 
QPSO-Type 2-II 2.6728e-35 1.4099 2.1737e+6 2.1835e+3 4.3398e+3 89.6543 
((6.5932e-35) (7.8582) (1.0089e+6) (2.8487e+3) (1.4313e+3) (151.6908) 
Algorithms F7 F8 F9 F10 F11 F12 
PSO-In 0.9907 0.0414 39.5528 239.5814 41.0529 3.6785e+4 
 (4.7802) (0.2393) (16.1654) (72.2521) (6.0318) (4.0943e+4) 
PSO-Co 0.0255 5.1120 96.7296 171.6488 36.0339 9.9648e+3 
 (0.0327) (4.5667) (28.0712) (58.5713) (7.2659) (1.6158e+4) 
Standard PSO 0.0218 0.2744 79.1219 128.9865 30.3424 1.8178e+4 
 (0.0165) (0.6795) (20.2619) (32.3662) (2.7409) (1.4866e+4) 
Gaussian PSO 0.0224 2.7722 103.6245 184.2657 33.5448 6.8875e+4 
 (0.0178) (1.4603) (28.6113) (57.3675) (6.5823) (6.5610e+4) 
Gaussian 0.0205 3.5460 80.9496 164.2914 29.8088 3.4327e+4 
BBPSO (0.0208) (6.1929) (22.0621) (72.8542) (3.2671) (6.2435e+4) 
PSO-E 0.0493 3.5881 66.5112 163.7187 29.2666 1.7161e+4 
 (0.0538) (5.5286) (20.9853) (55.0921) (3.2083) (1.0862e+4) 
Lévy PSO 0.0446 2.2168 74.0446 154.3838 28.9923 1.6282e+4 
 (0.1182) (1.3575) (21.6913) (76.3070) (5.0212) (2.5184e+4) 
QPSO-Type 1 0.0218 0.1217 56.4232 137.0334 28.2096 1.2145e+4 
( 1.00) (0.0204) (0.4504) (16.7090) (38.5269) (3.0216) (9.7844e+3) 
QPSO-Type 1 0.0209 0.0916 54.4278 126.1298 29.4137 1.0576e+4 
((0.0203) (0.3166) (16.6044) (44.9531) (2.8907) (9.0572e+3) 
QPSO-Type 2-I 0.0203 0.0683 39.0991 128.5351 19.8616 7.2794e+3 
( 0.75) (0.0164) (0.3080) (12.4904) (57.6255) (7.0620) (8.2210e+3) 
QPSO-Type 2-I 0.0208 2.0961e-14 29.9218 118.4549 28.1887 1.2938e+4 
((0.0130) (1.9099e-14) (10.5736) (53.0216) (6.2233) (1.3787e+4) 
QPSO-Type 2-II 0.0163 0.0762 42.4817 185.6351 19.0976 4.6519e+3 
( 0.54) (0.0134) (0.3075) (12.1384) (46.6356) (6.7920) (4.3177e+3) 
QPSO-Type 2-II 0.0150 7.5318e-15 43.8327 207.0548 23.0303 5.6134e+3 
((0.0119) (1.7046e-15) (17.881) (14.4658) (8.7874) (4.6511e+3) 
AlgorithmsF1F2F3F4F5F6
PSO-In 3.8773e-13 785.0932 3.9733e+7 1.1249e+4 6.0547e+3 263.7252 
 (1.6083e-12) (661.2154) (4.6433e+7) (5.4394e+3) (2.0346e+3) (437.4145) 
PSO-Co 1.5713e-26 0.1267 8.6472e+6 1.3219e+4 7.6892e+3 123.0243 
 (1.4427e-25) (0.3796) (9.1219e+6) (6.0874e+3) (2.3917e+3) (266.2520) 
Standard PSO 8.2929e-26 78.2831 6.6185e+6 1.3312e+4 6.2884e+3 153.5178 
 (1.2289e-25) (52.3272) (3.0124e+6) (4.1076e+3) (1.4318e+3) (246.1049) 
Gaussian PSO 7.3661e-26 0.0988 1.1669e+7 2.3982e+4 8.0279e+3 150.7872 
 (5.9181e-25) (0.3362) (2.5153e+7) (1.2512e+4) (2.3704e+3) (303.3368) 
Gaussian 1.7869e-25 16.8751 7.7940e+6 1.1405e+4 9.5814e+3 144.1377 
BBPSO (8.4585e-25) (16.2021) (4.3240e+6) (6.7712e+3) (3.0227e+3) (165.2616) 
PSO-E 5.2531e-24 20.2750 6.2852e+6 8.2706e+3 7.2562e+3 189.8292 
 (2.2395e-23) (15.2414) (2.8036e+6) (3.6254e+3) (1.8666e+3) (375.8636) 
Lévy PSO 1.1880e-24 36.9986 1.7366e+07 7.4842e+3 8.2543e+3 133.9526 
 (1.1455e-23) (29.1360) (1.9001e+7) (6.6588e+3) (2.2297e+3) (293.8460) 
QPSO-Type 1 3.5936e-28 40.2282 4.8847e+6 6.2397e+3 8.0749e+3 138.0746 
( 1.00) (1.5180e-28) (23.3222) (2.1489e+6) (2.4129e+3) (1.7099e+3) (209.1735) 
QPSO-Type 1 5.0866e-29 4.5003 3.2820e+6 6.4303e+3 7.8471e+3 139.9815 
((4.4076e-29) (2.9147) (1.9953e+6) (2.9744e+3) (1.7878e+3) (206.8138) 
QPSO-Type 2-I 1.9838e-27 0.1771 1.6559e+6 3.1321e+3 5.7853e+3 82.9908 
(0.75) (5.2716e-28) (0.1137) (7.1264e+5) (2.0222e+3) (1.2483e+3) (119.836) 
QPSO-Type 2-I 1.2672e-27 120.6051 4.4257e+6 4.0049e+3 3.3684e+3 88.0494 
((3.7147e-28) (62.2340) (2.3302e+6) (2.7218e+3) (975.6551) (159.7481) 
QPSO-Type 2-II 3.1554e-36 0.0715 1.8544e+6 3.1443e+3 5.7144e+3 105.7474 
( 0.54) (2.3913e-36) (0.0530) (6.4710e+5) (3.8785e+3) (1.4898e+3) (155.4583) 
QPSO-Type 2-II 2.6728e-35 1.4099 2.1737e+6 2.1835e+3 4.3398e+3 89.6543 
((6.5932e-35) (7.8582) (1.0089e+6) (2.8487e+3) (1.4313e+3) (151.6908) 
Algorithms F7 F8 F9 F10 F11 F12 
PSO-In 0.9907 0.0414 39.5528 239.5814 41.0529 3.6785e+4 
 (4.7802) (0.2393) (16.1654) (72.2521) (6.0318) (4.0943e+4) 
PSO-Co 0.0255 5.1120 96.7296 171.6488 36.0339 9.9648e+3 
 (0.0327) (4.5667) (28.0712) (58.5713) (7.2659) (1.6158e+4) 
Standard PSO 0.0218 0.2744 79.1219 128.9865 30.3424 1.8178e+4 
 (0.0165) (0.6795) (20.2619) (32.3662) (2.7409) (1.4866e+4) 
Gaussian PSO 0.0224 2.7722 103.6245 184.2657 33.5448 6.8875e+4 
 (0.0178) (1.4603) (28.6113) (57.3675) (6.5823) (6.5610e+4) 
Gaussian 0.0205 3.5460 80.9496 164.2914 29.8088 3.4327e+4 
BBPSO (0.0208) (6.1929) (22.0621) (72.8542) (3.2671) (6.2435e+4) 
PSO-E 0.0493 3.5881 66.5112 163.7187 29.2666 1.7161e+4 
 (0.0538) (5.5286) (20.9853) (55.0921) (3.2083) (1.0862e+4) 
Lévy PSO 0.0446 2.2168 74.0446 154.3838 28.9923 1.6282e+4 
 (0.1182) (1.3575) (21.6913) (76.3070) (5.0212) (2.5184e+4) 
QPSO-Type 1 0.0218 0.1217 56.4232 137.0334 28.2096 1.2145e+4 
( 1.00) (0.0204) (0.4504) (16.7090) (38.5269) (3.0216) (9.7844e+3) 
QPSO-Type 1 0.0209 0.0916 54.4278 126.1298 29.4137 1.0576e+4 
((0.0203) (0.3166) (16.6044) (44.9531) (2.8907) (9.0572e+3) 
QPSO-Type 2-I 0.0203 0.0683 39.0991 128.5351 19.8616 7.2794e+3 
( 0.75) (0.0164) (0.3080) (12.4904) (57.6255) (7.0620) (8.2210e+3) 
QPSO-Type 2-I 0.0208 2.0961e-14 29.9218 118.4549 28.1887 1.2938e+4 
((0.0130) (1.9099e-14) (10.5736) (53.0216) (6.2233) (1.3787e+4) 
QPSO-Type 2-II 0.0163 0.0762 42.4817 185.6351 19.0976 4.6519e+3 
( 0.54) (0.0134) (0.3075) (12.1384) (46.6356) (6.7920) (4.3177e+3) 
QPSO-Type 2-II 0.0150 7.5318e-15 43.8327 207.0548 23.0303 5.6134e+3 
((0.0119) (1.7046e-15) (17.881) (14.4658) (8.7874) (4.6511e+3) 

The algorithms were ranked to determine which algorithm could reliably be said to be the most effective for each problem. The algorithms that were not statistically different from each other were given the same rank; those that were not statistically different from more than one other group of algorithms were ranked with the best-performing of these groups. For each algorithm, the resulting rank for each problem and the total rank are shown in Table 5.

Table 5:
Ranking by algorithms and problems.
Total
AlgorithmsF1F2F3F4F5F6F7F8F9F10F11F12rank
PSO-In 13 13 13 =9 =3 =12 13 =1 =2 13 13 =11 116 
PSO-Co =6 =1 =9 =11 =7 =1 =3 13 =12 =9 12 =3 87 
Standard PSO =6 11 =7 =11 =3 =1 =3 =7 =9 =1 10 =9 78 
Gaussian PSO =6 =1 =9 13 10 =1 =3 =9 =12 =9 11 13 97 
Gaussian BBPSO =10 =7 =9 =9 13 =1 =3 =9 =9 =6 =4 =11 91 
PSO-E 12 =7 =7 =7 =12 =11 =9 =6 =4 =9 100 
Lévy PSO =10 =9 12 =5 =10 =1 =11 =9 =9 =6 =4 =5 91 
QPSO-Type 1 =9 =5 =5 =10 =1 =3 =7 =6 =4 =5 64 
( 1.00)              
QPSO-Type 1 =5 =7 =1 =3 =1 =6 =1 =4 =5 46 
(             
QPSO-Type 2-I =6 =4 =2 =3 =1 =3 =1 =2 =1 =1 =3 28 
( 0.75)              
QPSO-Type 2-I 12 =5 =2 =1 =3 =1 =1 =4 =5 41 
(             
QPSO-Type 2-II =1 =2 =3 =1 =1 =1 =2 =9 =1 =1 25 
( 0.54)              
QPSO-Type 2-II =4 =1 =1 =1 =2 12 =1 33 
(             
Total
AlgorithmsF1F2F3F4F5F6F7F8F9F10F11F12rank
PSO-In 13 13 13 =9 =3 =12 13 =1 =2 13 13 =11 116 
PSO-Co =6 =1 =9 =11 =7 =1 =3 13 =12 =9 12 =3 87 
Standard PSO =6 11 =7 =11 =3 =1 =3 =7 =9 =1 10 =9 78 
Gaussian PSO =6 =1 =9 13 10 =1 =3 =9 =12 =9 11 13 97 
Gaussian BBPSO =10 =7 =9 =9 13 =1 =3 =9 =9 =6 =4 =11 91 
PSO-E 12 =7 =7 =7 =12 =11 =9 =6 =4 =9 100 
Lévy PSO =10 =9 12 =5 =10 =1 =11 =9 =9 =6 =4 =5 91 
QPSO-Type 1 =9 =5 =5 =10 =1 =3 =7 =6 =4 =5 64 
( 1.00)              
QPSO-Type 1 =5 =7 =1 =3 =1 =6 =1 =4 =5 46 
(             
QPSO-Type 2-I =6 =4 =2 =3 =1 =3 =1 =2 =1 =1 =3 28 
( 0.75)              
QPSO-Type 2-I 12 =5 =2 =1 =3 =1 =1 =4 =5 41 
(             
QPSO-Type 2-II =1 =2 =3 =1 =1 =1 =2 =9 =1 =1 25 
( 0.54)              
QPSO-Type 2-II =4 =1 =1 =1 =2 12 =1 33 
(             

For the Shifted Sphere Function (F1), QPSO-Type 2-II with either fixed or time-varying generated better results than other methods. The results for Shifted Schwefel's Problem 1.2 (F2) show that QPSO-Type 2-II with fixed and PSO-Co got the best results but the performances of PSO-In and QPSO-Type 2-I with linearly decreasing were inferior to those of the other competitors. For the Shifted Rotated High Conditioned Elliptic Function (F3), the QPSO-Type 2-I with fixed seemed to outperform the other methods at a level of statistical significance. QPSO-Type 2-II showed to be the winner among all the tested algorithms for Shifted Schwefel's Problem 1.2 with Noise in Fitness (F4). F5 is Schwefel's Problem 2.6 with Global Optimum on the Bounds, and for this benchmark, QPSO-Type 2-I with linearly decreasing yielded the best results. For benchmark F6, Shifted Rosenbrock Function, the performance of PSO-In and PSO-E was inferior to those of the other algorithms, among which there was no statistically significant difference. The results for the Shifted Rotated Griewank's Function without Bounds (F7) suggest that QPSO-Type 2-II, either with fixed or with time-varying , was able to find the solution for the function with the best quality among all the algorithms. Benchmark F8 is Shifted Rotated Ackley's Function with Global Optimum on the Bounds. All the QPSO methods except QPSO-Type 1 with time-varying , along with PSO-Co, showed a better performance for this problem than the others. As far as Shifted Rastrigin's Function (F9) is concerned, the QPSO-Type 2-I with linearly decreasing yielded the best result, and the results obtained by PSO-In and all other QPSO-based methods rank tied for the second best. F10 is Shifted Rotated Rastrigin's Function which appears to be a more difficult problem than F9. For this benchmark, QPSO-Type 2-II, QPSO-Type 1 with time-varying and the standard PSO outperformed the other competitors in a statistically significant manner. The best results for Shifted Rotated Weierstrass Function (F11) were obtained by QPSO-Type 2 when was fixed during the search process. When searching the optima of Schewefel's Problem 2.13 (F12), QPSO-Type 2-II, whether it employed fixed or time-varying , was found to have the best performance.

As shown by the total ranks listed in Table 5, the methods based on QPSO-Type 2 obtained a better overall performance than all the other tested algorithms. For each of the benchmark functions, their performance were as good as or significantly better than that of the other algorithms. It is also revealed by the total ranks that QPSO-Type 2-II performed slightly better than QPSO-Type 2-I. However, more detailed comparisons reveal that for F10, the results obtained by QPSO-Type 2-II had poorer quality than those by QPSO-Type 2-I even than those by other PSO variants, and that the QPSO-Type 2-I with fixed had the most stable performance across all of the benchmark functions with the worst rank being 6 for F6. It can also be observed that for either version of QPSO-Type 2, time-varying led the algorithm to an overall performance no better than for fixed , particularly for QPSO-Type 2-I.

The second best performing algorithm was the QPSO-Type 1 algorithm as indicated by the total ranks. It can be seen that there is a great deal of difference between the performance of the two parameter control methods. However, in contrast with QPSO-Type 2, QPSO-Type 1 appeared to work better by using the time-varying control method. The standard PSO was the next best algorithm, and it enhances the search ability of PSO by incorporating the lbest ring topology into PSO-Co, which was ranked fourth in overall performance. It is evident from the ranking list that the standard PSO and PSO-Co were both a great improvement over the PSO-In algorithm. The other four probabilistic algorithms did not work so effectively as the QPSO algorithms and the standard PSO. Among the four methods, Gaussian BBPSO and Lévy PSO showed better search ability than PSO-E and Gaussian PSO.

5.  Conclusion

In this paper, after describing the background of the QPSO algorithm, which belongs to the bare-bones PSO family, we theoretically analyzed the behavior of the individual particle in QPSO. Then, through empirical studies on a well-known benchmark suite, we provided guidelines for parameter selection for the algorithm, as well as the performance comparison between QPSO and some other forms of PSO.

Two types of particles, corresponding to two search strategies of the QPSO algorithm, were analyzed by using the theory of probability measure. We derived the sufficient and necessary condition for the position of a single particle to be convergent or probabilistically bounded. Since the behavior of the particle is determined by the CE coefficient , the derivation of the condition is reduced to find out how to select the value of to guarantee the convergence or probabilistic boundedness of the particle's position. For a Type-1 particle, if , the position sequence of the particle converges to its local attractor in probability (or almost surely, in distribution, in mean of order r); if , the position is probabilistically bounded; otherwise, the position diverges. For a Type-2 particle, if , the position is probabilistically bounded, or else it diverges. Therefore, provides an upper bound for the value of