Feedback controlled particle swarm optimization and its application in time-series prediction

Expert Systems with Applications 39 (2012) 8557–8572

Contents lists available at SciVerse ScienceDirect

Expert Systems with Applications journal homepage: www.elsevier.com/locate/eswa

Feedback controlled particle swarm optimization and its application in time-series prediction q W.K. Wong ⇑, S.Y.S. Leung, Z.X. Guo Institute of Textiles and Clothing, The Hong Kong Polytechnic University, Hong Kong, PR China

a r t i c l e

i n f o

Keywords: Particle swarm optimization Feedback control Time-series prediction

a b s t r a c t Particle swarm optimization (PSO) algorithm is an algorithmic technique for optimization by solving a wide range of optimization problems. This paper presents a new approach of extending PSO to solve optimization problems by using the feedback control mechanism (FCPSO). The proposed FCPSO consists of two major steps. First, by evaluating the fitness value of each particle, a simple particle evolutionary fitness function is designed to control parameters involving acceleration coefficient, refreshing gap, learning probabilities and number of the potential exemplars automatically. By such a simple particle evolutionary fitness function, each particle has its own search parameters in a search environment. Secondly, a local learning method using a competitive penalized method is developed to refine the solution. The FCPSO has been comprehensively evaluated on 18 unimodal, multimodal and composite benchmark functions with or without rotation. Compared with various state-of-the-art algorithms, including traditional PSO algorithms and representative variants of PSO algorithms, the performance of FCPSO is promising. The effects of parameter adaptation, parameter sensitivity and local search method are studied. Lastly, the proposed FCPSO is applied to constructing a radial basis neural network, together with the K-means method for time-series prediction. Ó 2012 Elsevier Ltd. All rights reserved.

1. Introduction Particle swarm optimization (PSO) is a population-based search optimization technique proposed by Kennedy and Eberhart (1995) in 1995. The PSO utilizes a simple mechanism that mimics swarm behavior in agents such as birds flocking and fish assembling to guide the particles to search for globally optimal solutions. In PSO, a swarm consists of a number of particles. Each particle stands for a potential solution of the optimization problem. All of the particles iteratively discover a probable solution. Each particle flies to a new position according to the new velocity and the previous positions of the particle. The PSO algorithm is easy to implement and has rapidly progressed in recent years (Andrews, 2006; Angeline, 1998; Bergh & Engelbrecht, 2004; Cai, Cui, Zeng, & Tana, 2008; Chen, Peng, & Jian, 2007; Clerc & Kennedy, 2002; Coello, Pulidos, & Lechuga, 2004; Eberhart & Shi, 2001; Kennedy, 1997, 1999; Kennedy & Eberhart, 1995; Kennedy & Mendes, 2002, 2006; Krohling & dos Santos Coelho, 2006; Liang, Suganthan, & Baskar, 2006; Lovbjerg, Rasmussen, & Krink, 2001; Mendes, Kennedy, &

q The work described in this paper was fully supported by a research Grant from The Hong Kong Polytechnic University (Project No. G-YH11). ⇑ Corresponding author. E-mail addresses: [email protected] (W.K. Wong), tcleungs@inet. polyu.edu.hk (S.Y.S. Leung), [email protected] (Z.X. Guo).

0957-4174/$ - see front matter Ó 2012 Elsevier Ltd. All rights reserved. doi:10.1016/j.eswa.2012.01.126

Neves, 2004; Ratnaweera, Halgamure, & Watson, 2004; Shi & Eberhart, 1998a, 1998b; Tang, Wang, & Fang, 2011; Yen & Leong, 2009; Zhan, Zhang, Li, & Chung, 2009). However, the conventional PSO algorithm can easily get trapped in the local optima when solving complex multimodal tasks. These weaknesses have restricted applications of the PSO in real-world optimization problems. Therefore, accelerating convergence rates and avoiding the local optima are two important and appealing aims in the study of PSO. A number of variants of PSO algorithms have been proposed to achieve these two aims (Andrews, 2006; Angeline, 1998; Bergh & Engelbrecht, 2004; Cai et al., 2008; Chen et al., 2007; Clerc & Kennedy, 2002; Eberhart & Shi, 2001; Kennedy, 1997, 1999; Kennedy & Eberhart, 1995; Kennedy & Mendes, 2002, 2006; Krohling & dos Santos Coelho, 2006; Liang et al., 2006; Lovbjerg et al., 2001; Mendes et al., 2004; Ratnaweera et al., 2004; Shi & Eberhart, 1998a, 1998b; Zhan et al., 2009). One of the most popular ways to achieve the goals are to use time-varying controlling strategies, which are based on the generation number in the PSO iterations using linearly decreasing rules (Shi & Eberhart, 1998a). Although these generation-number-based methods lead to some improvement on the PSO algorithm, they might give rise to the risk of inappropriately tuning the parameters since no information on the evolutionary state that reflects the population and fitness diversity is identified or utilized. In order to improve efficiency and to accelerate the search process, it is of great importance

8558

W.K. Wong et al. / Expert Systems with Applications 39 (2012) 8557–8572

to determine the most suitable values for the parameters including inertia weights and acceleration coefficients. To improve the performance of PSO on complex optimization problems, a particle swarm optimizer utilizing the feedback control method is developed in this paper. The PSO parameters including acceleration coefficients, refreshing gap, number of exemplars and learning probabilities are controlled by our designed particle’s evolutionary function. A leader local learning method is proposed in this paper to perform only on the globally best particle by utilizing a competitive penalized method. A series of experiments are to be carried out on various benchmark functions with different topological structures to verify the effectiveness of the FCPSO and to comprehensively compare FCPSO with other variants of PSO algorithms. On the other hand, as typical kind of neural networks (Wong, Zeng, & Au, 2009, 2009), although radial basis function neural networks (RBFNNs) have been widely used in time-series prediction (Buchtala, Klimek, & Sick, 2005; Du & Zhang, 2008; Harpham & Dawson, 2006; Makridakis, Wheelwright, & Hyndman, 1998; Zhang, 2003), it is still difficult to select an appropriate network structure (Buchtala et al., 2005; Chen, Wu, & Luk, 1999; Du & Zhang, 2008). In particular, to overcome the existing difficulties in developing an RBFNN, evolving RBFNN combining genetic algorithms (GAs) with other standard learning algorithms for constructing RBFNNs has been reported in the literature (for example see Buchtala et al., 2005; Chen et al., 1999; Du & Zhang, 2008). In addition, a number of input variables are available for model construction in many cases since the amount of collected data increases rapidly (Guyon & Elisseeff, 2003). In Du and Zhang (2008), an evolving RBFNN using a GA with binary codes has been proposed, where centers, widths and inputs are selected automatically. However, the encoding scheme is complex and the population size of GAs needs to be large to achieve high prediction accuracy. As a promising evolutionary algorithm, PSO is a competitive algorithm for constructing RBFNN. In this paper, we aim to utilize the proposed PSO to construct a RBFNN with a new simple encoding scheme. With the proposed encoding scheme, a novel RBFNN is applied to time-series prediction. In this paper, we aim to utilize the proposed PSO to construct a RBFNN with a simple encoding scheme. With the proposed encoding scheme, a novel RBFNN is applied to time-series prediction. This paper is organized as follows. Section 2 introduces the original PSO and its current variants. Section 3 describes the proposed PSO using feedback control method. Section 4 presents the test functions and the experimental setting for each algorithm. A time-series prediction method using the novel PSO is proposed in Section 5. Conclusions are given in Section 6. 2. PSO algorithms In the PSO, each particle moves to a new position according to the new velocity and the previous positions of the particle. This is compared with the best position generated by previous particles in the fitness function, and the best one is kept in the archive. Therefore, each particle is guided by both the local best solution and the global best position. If a particle moves and a better solution is found, other particles will move closer to it to explore the region in the process. 2.1. PSO algorithm In PSO, a swarm consists of S particles moving around in a D-dimensional search space. The position of the ith particle is represented by a vector xi(k) = (xi1(k), xi2(k), . . . , xiD(k)), where xij(k) 2 [xmin,j, xmax,j], 1 6 j 6 D, xmin,j and xmax,j are lower and upper

bounds for the jth dimension respectively. In the search process, the particle tunes its current position toward the global optimum according to two factors: the best position encountered by itself (pbest) denoted as pi = (pi1, pi2, . . . , piD) and the best position in the whole swarm (gbest) denoted as pg = (pg1, pg2, . . . , pgD). The velocity of the ith particle at the kth iteration is represented by vi(k) = (vi1(k), vi2(k), . . . , viD(k)), which is limited to a maximum velocity vmax = (vmax,1, vmax,2, . . . , vmax,D). Thus, if the magnitude of the updated velocity jvi,dj exceeds vmax,d, then vi,d is assigned the value sign(vi,d)vmax,d. In this paper, the maximum velocity vmax,d is set to 20% of the search range, as proposed in Eberhart and Shi (2001). The inertia weight w is a scaling factor controlling the influence of the old velocity on the new one which was proposed in Shi and Eberhart (1998a); r1,j and r2,j are two uniform random number samples from U(0, 1). The parameters c1 and c2 are called acceleration coefficients, called cognitive and social coefficients respectively. The velocity and position of the particle at next iteration are updated according to the following equations:

v i;j ðk þ 1Þ ¼ wðkÞv i ðkÞ þ c1 r1;j ðpi;j ðkÞ  xi;j ðkÞÞ þ c2 r2;j ðpg;j ðkÞ  xi;j ðkÞÞ; ð1Þ xi;j ðk þ 1Þ ¼ xi;j ðkÞ þ v i;j ðk þ 1Þ: 2.2. Some variants of PSO Since its introduction in 1995 by Kennedy and Eberhart (1995), PSO has attracted much research interest. A number of results have been developed to improve its performance in various ways, thereby deriving many interesting variants. One of the most famous variants (Shi & Eberhart, 1998a) introduces a parameter called inertia weight into the original PSO algorithms. The inertia weight w is utilized to balance the global and local search abilities. A large inertia weight is more appropriate for global search, and a small inertia weight facilitates local search. A linearly decreasing inertia weight over the course of search according to the generation number was proposed by Shi and Eberhart (1998a). Convergence analysis and stability studies have been reported by Clerc and Kennedy, and the constriction factor has been introduced into PSO for analyzing the convergence behavior (Clerc & Kennedy, 2002). Apart from the inertia weight and the constriction factor, the acceleration coefficients c1 and c2 are also vital to PSO. In two extreme cases conducted by Kennedy (1997), i.e., the social model and the cognitive model, experiments have shown that both acceleration coefficients are essential to PSO. It was suggested by Kennedy and Eberhart (1995) that the coefficients are fixed to 2.0, and this configuration has been adopted by many other researchers. Ratnaweera et al. (2004) proposed a PSO algorithm with linearly time-varying acceleration coefficients (PSO-TVAC), where a larger c1 and a smaller c2 were set at the beginning and were gradually reversed during the search. Similar to the inertia weight, the time-varying c1 and c2 can balance the global and local search abilities. In Zhan et al. (2009), an adaptive PSO (APSO) was presented by defining an evolutionary state function and an elite learning method was utilized in convergence rate for avoiding the particles trapping into the local optima. In addition to studying parameter control and auxiliary techniques, PSO topological structures are also widely investigated. It was found by Kennedy (1999, 2002, 2006) that PSO with a small neighborhood might perform better on complex problems, while PSO with a large neighborhood might perform better on simple problems. In the ‘‘fully informed particle swarm’’ (FIPS) algorithm (Mendes et al., 2004), the information of the entire neighborhood is employed to guide the particles. In Liang et al. (2006), a comprehensive learning PSO was proposed. Its learning strategy abandons the global best information, and all other particles’ past best positions might be used to update particles’ velocity.

W.K. Wong et al. / Expert Systems with Applications 39 (2012) 8557–8572

Hybridization by combining PSO with other search approaches is also a hot topic to improve the performance of the PSO. Evolutionary operators such as selection, crossover, and mutation have been introduced to PSO to keep the best particles (Andrews, 2006; Angeline, 1998), to increase the diversity of the population, and to improve the ability to escape local minima (Lovbjerg et al., 2001). In Lovbjerg et al. (2001), the swarm is divided into sub-populations, and a breeding operator is utilized within a subpopulation or between the sub-populations to increase the diversity of the population. A cooperative particle swarm optimizer (CPSO-H) (Bergh & Engelbrecht, 2004) was proposed to improve the performance of the original PSO on multimodal problems. 3. FCPSO Although PSO algorithms have been applied to a wide range of optimization problems and numerous variants of the PSO have been introduced, solving complex optimization problems with high efficiency and rapid convergence is still an active area of research. In order to allow better particles to achieve efficient search and rapid convergence speed of PSO, two techniques including feedback control scheme and local search method controlled by an efficient competitive penalized method are proposed. This section firstly introduces the positions and velocity updating equations of our method. Secondly, the control method of acceleration parameters using feedback mechanism is introduced. Thirdly, the control methods of learning probability, refreshing gap and number of exemplars are proposed. Finally, a leader local search is presented. 3.1. Positions and velocity updating equations of FCPSO The proposed FCPSO is given as follows:

v i;j ðk þ 1Þ ¼ wv i;j ðkÞ þ di ðkÞci ðkÞrj ðkÞðpr;j ðkÞ  xi;j ðkÞÞ þ ð1  di ðkÞÞci ðkÞr j ðkÞðpi;j ðkÞ  xi;j ðkÞÞ; xi;j ðk þ 1Þ ¼ xi;j ðkÞ þ v i;j ðk þ 1Þ;

ð2Þ ð3Þ

where pr,j(k) represents another particles’s history best fitness in the jth dimension (Liang et al., 2006); ci(k) is the acceleration coefficient which is time-varying according to the search circumstances. di(k) is stochastic variable that describes the following random events for the system (2):

(

Event 1 : system ð2Þ experiences pr;j ðkÞ; Event 2 : system ð2Þ experiences pi;j ðkÞ:

ð4Þ

Let di(k) be Bernoulli distributed sequences defined by

 di ðkÞ ¼

1; if Event 1 occurs; 0; if Event 2 occurs;

ð5Þ

where di(k) satisfies Prob{di(k) = 1} = di0, Prob{di(k) = 0} = 1  di0. In the above equations, di0 denotes the learning probability rate. The velocity update equation is written in the form of (2) to highlight the importance of di0 in the PSO searching process. The stochastic variable di(k) is introduced into velocities updating equations to describe different learning strategies occurring in a probabilistic way. Using binary sequence switching to illustrate stochastic events such as miss measurement, random packet losses and stochastic delay is common in some works (Tang, Fang, Xia, & Yu, 2009; Wang, Ho, & Liu, 2003, 2007). On the other hand, the stochastic variables used in our paper can describe the loss information conveyed in a flock of birds. The randomly occurring information losses are known in the process of communication between the particles in the swarm. The information losses are

8559

considered as a binary switching sequence that is specified by a conditional probability distribution. In our technique, the learning probability is time-varying according to the search environment. It was found that different learning-probability values will affect the search results for the same problem if the same value of learning probability was used for all the particles in the swarm (Liang et al., 2006). Thus, the learning probability of each particle is assigned differently. Similar to the method in Liang et al. (2006), for each dimension of particle i, a random number is generated. If the random number is larger than di0, the corresponding dimension will learn from its own pi,j(k); otherwise it will learn from another particle’s pi,j(k). To guarantee that a particle learns from good exemplars and to minimize the time wasted on poor ways, the particle is allowed to learn from the exemplars until the particle ceases improving for a certain number of generations called the refreshing gap mi(k); then the value of refreshing gap for the particle is reset. 3.2. Control of acceleration parameters using feedback mechanism The acceleration coefficient c is usually used to balance the global and local search abilities of PSO, which is vital to the search performance of PSO (Kennedy, 1997). In Ratnaweera et al. (2004), PSO with the acceleration coefficients was time-varying with the generation number is proposed. However, only time information is taken into account in the proposed scheme. The search information at each generation is not considered although such information is of great importance for enhancing the performance of PSO. From Eq. (2), it can be observed that a large acceleration coefficient will enable the particle to learn from other particles or its best position found so far quickly, while a small acceleration coefficient will make the particle retain in the current position and learn from other particles slowly to search for the particle’s nearby position. In this subsection, an efficient and simple method is proposed to automatically assign the acceleration parameters to each particle. It should be noted that the feedback mechanism was proposed in Cai et al. (2008) for controlling acceleration coefficient. However, different from the work in Cai et al. (2008), the feedback method is further extended to control other parameters for PSO in this paper. The fitness information is time-varying in a PSO evolutionary process and it characterizes the different levels of each particle. Each particle in the swarm plays a different role in the searching process. For instance, the worst particle (which has the worst fitness value) might become the best particle (which has best fitness) in the next generation according to the velocity updating equation. Therefore, how to employ the search information such as fitness information to control the acceleration parameters would be an important and promising research topic in PSO. Generally, if the worst particle has a large velocity to move toward the exemplar, it will enable itself to achieve global search rapidly. Meanwhile, if a small speed is assigned to the best particle, it would make the particle search around itself to refine the current search solution. At each generation, we can easily sort the current fitness value and obtain the rank of each particle. According to the discussion above, a novel, simple and efficient control method of acceleration coefficients is proposed by designing an evolutionary fitness function (EFF) in this paper. The acceleration coefficients can be calculated according to EFF easily. At each step, the best fitness of each particle i found so far is calculated and then a set of fitness values Fi (i = 1, 2, . . . , S) are obtained. We compare each Fi, and obtain the maximum and minimum value of Fmax and Fmin among Fi. The EFF can be measured using the set of fitness values:

Gi ¼

F i  F min ; F max  F min

ði ¼ 1; 2; . . . ; SÞ 2 ½0; 1:

ð6Þ

8560

W.K. Wong et al. / Expert Systems with Applications 39 (2012) 8557–8572

The EFF can be viewed as normalized fitness function at each iteration. From the above equation, it is clear to see that each particle might have different values of EFF at each generation. In this paper, the EFF can be used to determine the acceleration coefficient of each particle, that is, the worst particle has the largest speed to enable global search and the best particle has the smallest speed to refine the current solution. Hence, it would be beneficial to allow c to follow the EFF using the following equation cðGi Þ : Rþ ! Rþ :

ci ¼ ðc2  c1 Þ  Gi þ c1 ;

ði ¼ 1; 2; . . . ; SÞ;

ð7Þ

where c1 is the lower bound and c2 is the upper bound of the acceleration coefficient. The method is similar to the time-varying acceleration coefficient (Ratnaweera et al., 2004). However, the proposed method adjusts acceleration coefficient according to the population fitness. By using this linear transformation, the swarm has the acceleration coefficient belonging to the range [c1, c2] in the search process. Particles in the swarm play a different role in the evolutionary process. The worse particle has better global search abilities and the better particle has better local search abilities. Thus, the swarm can benefit from this method to have global and local search capabilities at each generation simultaneously. The acceleration coefficient of each particle will change its value according to its fitness instead of generation so as to adjust the value automatically. Under our design, since ci is not only monotonic with time, but also monotonic with Gi, ci will be adapted to the search circumstance characterized by Gi. For example, if a better solution is discovered outside the swarm, almost all the particles in the swarm will have large Gi and thus generate large ci, which will benefit the global search. Meanwhile, the best particle will stay comparatively statically to wait for other particles fly to it. If the swarm assembles for searching solutions more accurately, nearly all the particles will have small Gi and ci, which will lead the swarm to have the local search. It is easy to implement EFF in PSO since the calculation of Gi and ci does not need extra computation. In order to show the control scheme clearly, we utilize the method for solving the Griewank function with 20 particles and 30 dimensions. The change in the acceleration coefficient of one particle and mean acceleration coefficient of the swarm are illustrated in Fig. 1. It can be seen from Fig. 1 that the technique can dynamically tune the acceleration coefficient of the particle according to its fitness. It can be seen that, in the early phase, the acceleration coefficient is large so as to make velocity large and ensure the swarm a global search. In the latter search stage, the acceleration coefficient decreases and thus the swarm has the ability of searching locally.

3.3. Control of learning probability using feedback mechanism In this subsection, control scheme of learning probability is proposed using fitness feedback mechanism, which is similar to the method in controlling acceleration coefficients. From Liang et al. (2006), it is found that different values of learning probability yielded different results on the same problem if the same value was used for all the particles in the population. The authors proposed to set so that each particle has a different value. As a result, particles have different levels of exploration and exploitation ability in the swarm and are able to solve various kinds of problems. The learning probabilities for each particle were empirically developed the following expression (Liang et al., 2006):

exp di0 ¼ 0:05 þ 0:45 



10ði1Þ S1



1 ;

expð10Þ  1

ði ¼ 1; 2; . . . ; SÞ:

ð8Þ

In this paper, we define that the best particle has the minimum fitness value. Different from the fixed learning probability assigned to each particle proposed in Liang et al. (2006), the best fitness value Fi of each particle found so far is sorted in an ascending way. and remember the former position as Ii. Thus, the learning probability of each particle is assigned as follows:

exp dIi 0 ðGi Þ ¼ 0:05 þ 0:45 





10ði1Þ S1

1

expð10Þ  1

;

ði ¼ 1; 2; . . . ; SÞ:

ð9Þ

From the above learning methods, it can be found that the worst particle will have the largest learning probability to learn from other particles, while the best particle will have the smallest learning probability to retain the current position. By using this method, each particle has diverse learning probability along the search process and thus will produce different effects. An example is used to show the adaptive control of learning probability of one particle, as seen in Fig. 2. The example is for solving the Griewank function with 20 particles and 30 dimensions. The adjustment in the learning probability of one particle is depicted in Fig. 2. It can be observed from Fig. 2 that the technique can automatically alter the learning probability of the particle according to its fitness. One can see that the particle learns from other particles frequently in the early stage and retains the position to refine the solution accurately.

1.75 c1(k)

1.7

0.55

Mean ci(k)

δ10(k)

0.5

1.65 0.45 1.6 0.4 1.55 0.35 1.5 0.3 1.45 0.25 1.4 0.2 1.35 0.15 1.3 0.1 1.25

0

500

1000

1500

2000

generation

0.05

0

500

1000

1500

2000

generation Fig. 1. The parameter adaption process for acceleration coefficient of one particle and mean acceleration coefficient of the swarm.

Fig. 2. The parameter adaption process for learning probability of one particle.

8561

W.K. Wong et al. / Expert Systems with Applications 39 (2012) 8557–8572

3.4. Control of refreshing gap using feedback mechanism

3 Mean ni(k) of the swarm

In this subsection, a simple control technique of refreshing gap is proposed using the above feedback mechanism. The refreshing gap is another important factor in CLPSO to influence the results. The refreshing gap is set at seven for all test functions in Liang et al. (2006). However, this value cannot be applied to solve different optimization functions. Usually, it is much easier for the worst particle is much easier to update its best fitness than the best particle. Thus, a small refreshing gap is given to the worse particles and a large refreshing gap is assigned to the better particles. The mechanism means that the worse particles will learn from different directions rapidly and the better particles will be slower in learning from other directions. The worse particles tend to abandon the current search direction after fewer generations if the current direction cannot help it to find a better solution. Conversely, better particles have a tendency to keep the current search direction for more generations although the current direction cannot help it to find a better solution. The refreshing gap in this paper is updated using the following equation:

mi ðkÞ ¼ m  floorðk  Gi Þ;

ði ¼ 1; 2; . . . ; SÞ;

ð10Þ

where m is the maximum refreshing gap and k is an integer standing for the refreshing gap coefficient; floor(x) is the function for calculating the integer not larger than x. From (10), we can see that the refreshing gap ranges from [m  k, m]. In this paper, m is set to seven. By utilizing the above method, each particle possesses different refreshing gaps in the search process. An example is also presented to show the adaptive control of learning probability of one particle, as seen in Fig. 3. The adjustment in the refreshing gap of one particle for solving the Griewank function with 20 particles and 30 dimensions is depicted in Fig. 3. Fig. 3 shows that the technique can alter the refreshing gap of the swarm efficiently (see Fig. 4). 3.5. Control of number of exemplar using feedback mechanism In the original CLPSO, if a random number is larger than di0, the corresponding dimension will learn from its own pi,j(k); otherwise it will learn from another particle’s pi,j(k). A tournament selection procedure when the particle’s dimension learns from another particle’s pi,j(k) is proposed. Two particles are randomly chosen from the population and their fitness values are compared to select a better one to serve as an exemplar. However, we define the number of potential exemplars using fitness feedback mechanism as follows: 6.8 6.6

2.9 2.8 2.7 2.6 2.5 2.4 2.3 2.2 2.1

0

500

1000

1500

2000

generation Fig. 4. The parameter adaption process for mean number of potential exemplar.

ni ðkÞ ¼ 2 þ floorðl  Gi Þ;

ði ¼ 1; 2; . . . ; SÞ;

ð11Þ

where l is an integer coefficient to control the number of exemplars. From the above equation, each particle in the population has its own number of exemplars to follow. It can be observed that the number of potential exemplars belongs to [2, 2 + l]. Different from the method in Liang et al. (2006), if the random number is smaller than di0(k), ni(k) particles in the population are used and the fitness values of these ni(k) are compared in order to select a better/the best one to guide particle i. The worst particle will have the largest values of the number of potential exemplars, while the best particle will have smallest values of the number of potential exemplar. Benefiting from the control method, the bad particles will learn from the good particle rapidly while the good particle will maintain the diversity and slowly learn from other particles. We can see that the swarm can achieve global search and local search simultaneously. The swarm will have the abilities of quick convergence rate and will keep the diversity of the population at the same time. An example is given in Fig. 3 showing the adaptive control of learning probability of one particle. The adjustment in the number of potential exemplars for solving the Griewank function with a population of 20 particles and 30 dimensions is depicted in Fig. 3. The technique can alter the refreshing gap of the swarm efficiently. It can be observed that, in the early stage, the swarm needs more potential exemplars to achieve fast convergence speed and it needs fewer potential exemplars to have the ability of diversity in the latter stage. 3.6. Leader local search

Mean mi(k) of the swarm 6.4

In this subsection, a leader local search is used to refine the solution of the optimization problem. This approach is controlled automatically in the evolutionary process by a competitive penalized method. First, let s(k)(k P 0) be a Markov chain taking values in a finite state S ¼ f1; 2; . . . ; Ng with probability transition  space  ðkÞ matrix CðkÞ ¼ cij given by

6.2 6 5.8

NN

5.6

ðkÞ

P ¼ fsðk þ 1Þ ¼ jjsðkÞ ¼ ig ¼ cij ;

5.4

ðkÞ

5.2 5

0

500

1000

1500

2000

generation Fig. 3. The parameter adaption process for mean refreshing gap of the swarm.

i; j ¼ 1; 2; . . . ; N;

ð12Þ

where c P 0ði; j 2 SÞ is the transition rate from i to j and PN ðkÞ ij j¼1 cij ¼ 1. The leader local search randomly selects one dimension of the best particle, which is written as pg,j(k) for the jth dimension. Only one dimension is chosen since the local optima are likely to have better solution in one dimension. The leader local search is carried out by the following two search strategies

8562

W.K. Wong et al. / Expert Systems with Applications 39 (2012) 8557–8572

pg;j ðkÞ ¼ pg;j ðkÞ þ ððk1;j ðkÞ  k2;j ðkÞÞ  r1 ðkÞ þ k2;j ðkÞÞ  R1 ðkÞ;

ð13Þ

pg;j ðkÞ ¼ pg;j ðkÞ þ ððg1 ðkÞ  g2 ðkÞÞ  r2 ðkÞ þ g2 ðkÞÞ  R2 ðkÞ;

ð14Þ

where r1(k) and r2(k) are uniform random numbers samples from U(0, 1). k1,j(k) and k2,j(k) are the upper and lower bounds of the chosen particle in the swarm in jth dimension at each step. g1(k) and g2(k) are the maximal and minimal value of all dimensions of the globally best particle at each generation. R1(k) and R2(k) are the search range of the method which are updated by a probability matrix C(k). For choosing the R1(k) and R2(k), consider the following   cðkÞ 1  cðkÞ ; where c(k) is probability transition matrix CðkÞ ¼ ðkÞ ðkÞ c 1c time-varying to choose the radius. s(k) = 1 denotes that the Ri(k)(i = 1, 2) is set to 1 and s(k) = 2 indicates that the Ri(k) is chosen as 0.1. Here, we suggested that c(k) is linearly decreased with time like linearly decreased inertia weight, which is given by

cðkÞ ¼ ðc1  c2 Þ 

kmax  k þ c2 ; kmax

ð15Þ

where c1 and c2 are the upper and lower bounds of c(k). c1 = 1 and c2 = 0 are set in all experiments in this paper, which reveals that the Ri(k)(i = 1, 2) = 1 occurs frequently and Ri(k)(i = 1, 2) = 0.1 seldom happens in the early stage. In leader local search, the new position can be accepted if the fitness is better than the current globally best position. Otherwise, the new position is abandoned. A competitive penalized learning method is proposed to deal with the problem of these two search methods are not always useful in the search process. Define Q1 and Q2 as (13) activation variable and (14) activation variable, respectively. The activation probability variable Qi is updated as follows:

8 Q i þ j1 ; if the search strategy i > > > < succeed to update the globally best particle; Qi ¼ > Q i  j2 ; if the search strategy i > > : fails to improve the globally best particle; ð16Þ where i = 1, 2 represents the search strategy (13) and (14) respectively. j1 is the learning rate and j2 is the penalized rate. The j1 = 0.1 and j2 = 0.001. If the magnitude of the updated Qi(i = 1, 2) exceeds the value of 1, then Qi is assigned the value of 1. If Qi decreases under 0.1, then Qi is set to the value 0.1 for reoccurring.

3.7. Differences between CLPSO and the proposed PSO In addition to the leader local search, there exist four major differences in the learning strategies between CLPSO (Liang et al., 2006) and the proposed learning methods. The main differences arise from the usage of available fitness information. In this paper, the fitness information of each particle is used to control the parameters in PSO automatically. The four major differences are listed as follows: (1) Instead of using the fixed acceleration coefficient c, we employ the EFF value for feedback control acceleration coefficient c. Each particle has its own ci(k) and thus will lead to rapid convergence speed. (2) Instead of using unchangeable learning probability for each particle, each particle is assigned unique learning probability using fitness feedback mechanism. (3) Instead of assigning unalterable refreshing gap for all particles, each particle is given a refreshing gap by the idea of feedback control. (4) Instead of using two potential exemplars for the particles, the number of potential exemplars is calculated for each particle in the swarm using fitness feedback control. In summary, benefiting from the feedback control method of parameters, the particles can achieve global search and local search simultaneously. Thus, the swarm will have the capabilities of fast convergence rate and keep the diversity of the population at the same time. The following experimental results show the good performance of the proposed method. 4. Experimental results 4.1. Experiments setup Eighteen benchmark functions are listed in Table 1 and (17)–(32) are used to test the performance of PSOs (Yao, Liu, & Lin, 1999; Liang, Suganthan, & Deb, 2005; Suganthan et al., 2005). Biased initializations were used for the functions whose global optimum is at the center of the search range. All the functions were tested on 30 dimensions. The population size of all the PSOs was set to 20. f1(x), f2(x) and f9(x) represented unimodal optimization problem. f1(x) was used to test the convergence speeds of PSOs. f9(x) was unimodal function with rotation. f3(x) to f8(x) were multimodal problems which are hard to optimize. Note that some functions are separable and can be solved by using D one-dimensional searches. Hence, eight

Table 1 Benchmark configurations. Functions

Name

Dimension

Search space

Minimum

Threshold

Initial range

f1(x) f2(x) f3(x) f4(x) f5(x) f6(x) f7(x) f8(x) f9(x) f10(x) f11(x) f12(x) f13(x) f14(x) f15(x) f16(x) f17(x) f18(x)

Sphere Rosenbrock Rastrigin Noncontinuous Rastrigin Ackley Griewank Penalized Weierstrass Rotated Quadric Rotated Rosenbrock Rotated Rastrigin Rotated Noncontinuous Rastrigin Rotated Ackley Rotated Griewank Rotated Penalized Rotated Weierstrass CF3 Rotated CF3

30 30 30 30 30 30 30 30 30 30 30 30 30 30 30 30 30 30

[100, 100]D [10, 10]D [5, 5]D [5, 5]D [32, 32] [600, 600]D [50, 50]D [0.5, 0.5]D [100, 100]D [10, 10]D [5, 5]D [5, 5]D [32, 32]D [600, 600]D [10, 10]D [0.5, 0.5]D [500, 500]D [500, 500]D

0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

0.01 100 50 50 0.01 0.01 0.01 0.01 100 100 100 100 0.01 100 0.01 10 100 100

[100, 50]D [10, 5]D [5, 2]D [5, 2]D [32, 16]D [600, 200]D [50, 20]D [0.5, 0.2]D [100, 50]D [10, 5]D [5, 2]D [5, 2]D [32, 16]D [600, 200]D [10, 5]D [0.5, 0.2]D [5, 5]D [5, 5]D

8563

W.K. Wong et al. / Expert Systems with Applications 39 (2012) 8557–8572 Table 2 PSO algorithms for comparison.

Rotated Rastrigin : f 11 ðxÞ ¼

Parameters

Reference

LDIW TVAC CLPSO APSO FCPSO

w:0.9  0.4, c1 = c2 = 2 w:0.9  0.4, c1:2.5  0.5, c2:0.5  2.5 w:0.9  0.4, c = 1.49, m = 7 Automatically chosen w:0.65

(Shi & Eberhart, 1998a) (Ratnaweera et al., 2004) (Liang et al., 2006) (Zhan et al., 2009) this paper

y ¼ M  x;

x2i ;

ð17Þ

i¼1

Rosenbrock : f 2 ðxÞ ¼

D1 X ð100ðxiþ1  xi Þ2 þ ðxi  1Þ2 Þ;

ð18Þ

i¼1

Rastrigin : f 3 ðxÞ ¼

D X

ðx2i  10 cosð2pxi Þ þ 10Þ;

i¼1

Noncontinuous Rastrigin : ( where yi ¼

D X 

f 4 ðxÞ ¼

ð19Þ

 y2i  10 cosð2pyi Þ þ 10 ;

i¼1

xi ;

jxi j < 0:5;

roundð2xi Þ ; 2

jxi j P 0:5; rffiffiffiffiffiffiffiffiffiffiffi D P 1 2

0:2

xi

D

Ackley : f 5 ðxÞ ¼ 20e

i¼1

ð20Þ

1 D

e

D P

f12 ðxÞ ¼

D X 

where zi ¼

i¼1

þ 20 þ e;

Griewank : f 6 ðxÞ ¼

1 4000

Generalized Penalized : 2

f7 ðxÞ¼ Dp f10 sin ðpy1 Þþ þ

D P

x2i 

D Y

i¼1

i¼1

D1 P

2

  xi cos pffi þ 1; i

f1

f2

f3

f4

f5

f8

2

ðyi 1Þ ½1þ10 sin ðpyiþ1 ÞþðyD 1Þ g

i¼1

f10

uðxi ; 10; 100; 4Þ;

i¼1 8 m xi >a; < kðxi aÞ ;   1 where yi ¼ 1þ 4 ðxi þ1Þ ; uðxi ; a; k; mÞ¼ 0; a6xi 6 a; : kðxi  aÞm ; xi < a: ð23Þ

f11

f12

f13

Weierstrass :

! kmax kmax D X X X k k k k ½a cosð2pb ðxi þ 0:5ÞÞ  D ½a cosð2pb  0:5Þ ; f 8 ðxÞ ¼ i¼1

k¼0

Rotated Quadric : f 9 ðxÞ ¼

ð24Þ

N i X X i¼1

!2 xj

f15

f16

;

y ¼ M  x;

ð25Þ

j¼1

f17

D1 X Rotated Rosenbrock : f 10 ðxÞ ¼ ð100ðyiþ1  yi Þ2 þ ðyi  1Þ2 Þ; i¼1

y ¼ M  x;

f14

k¼0

a ¼ 0:5; b ¼ 3; kmax ¼ 20;

ð26Þ

jyi j < 0:5;

roundð2yi Þ ; 2

jyi j P 0:5; rffiffiffiffiffiffiffiffiffiffiffi D P 1 2 yi

D

ð28Þ

y ¼ M  x;

i¼1

e

1 D

D P

cos 2pyi

i¼1

ð29Þ ð30Þ

Table 3 Search result comparisons among five PSOs on 18 test functions.

f9 2

yi ;

þ 20 þ e;   D D Y 1 X y Rotated Griewank : f 14 ðxÞ ¼ y2  cos piffi þ 1; y ¼ M  x; 4000 i¼1 i i¼1 i

ð21Þ ð22Þ

(

Rotated Ackley : f 13 ðxÞ ¼ 20e

f7 D X

 z2i  10 cosð2pzi Þ þ 10 ;

0:2

f6 cos 2pxi

ð27Þ

Rotated noncontinuous Rastrigin :

i¼1

rotated unimodal and multimodal problems are used to test the performance of the PSOs. To rotate a function, first an orthogonal matrix M is generated according to the Salomon’s method (Salomon, 1996). The original variable x is left multiplied by the orthogonal matrix M to get the new rotated variable y = M ⁄ x. This variable y is used to compute the fitness value f. Clearly, when one dimension in x is changed, all dimensions in y will be influenced. Thus, the rotated function cannot be solved by D one-dimensional searches. f17 and f18 are composite functions with or without rotation. The composition functions are asymmetrical multimodal problems, with different properties in different areas (Liang et al., 2005). D X

ðy2i  10 cosð2pyi Þ þ 10Þ;

i¼1

Algorithm

Sphere : f 1 ðxÞ ¼

D X

f18

Mean Std. Dev. Mean Std. Dev. Mean Std. Dev. Mean Std. Dev. Mean Std. Dev. Mean Std. Dev. Mean Std. Dev. Mean Std. Dev. Mean Std. Dev. Mean Std. Dev. Mean Std. Dev. Mean Std. Dev. Mean Std. Dev. Mean Std. Dev. Mean Std. Dev. Mean Std. Dev. Mean Std. Dev. Mean Std. Dev.

LDIW

TVAC

CLPSO

APSO

FCPSO

3.89e11 1.94e10

5.4e8 2.5e7

8.91e29 8.37e29

3.55e43 4.81e45

1.10e052 3.46e052

837.37 2771.8

35.75 29.52

18.93 5.56

27.4 12

19.25 11.76

43.48 16.32

35.66 11.13

2.02e14 4.53e14

2.4 1.1

0 0

19.32 10.011

28.08 10.90

0 0

0.002 2.46e04

0 0

0.79 3.97

0.14 0.42

2.21e14 4.38e15

1.25e5 1.67e4

7.69e015 1.24e015

0.016 0.018

0.04 0.11

3.65e14 1.74e13

0.02 0.01

0 0

3.54e31 4.24e59

1.66 1.57e29

6.20e31 1.48e59

0.0003 0.0009

1.5e32 1.15e33

0.089 0.321

0.17 0.13

7.25e11 1.07e10

0.08 0.06

7.13e11 1.01e010

9668.18 9464.89

3050.92 5580.1

845.55 387.08

2.9e4 3.2e5

55.20 24.8

197.16 388.71

62.60 49.04

29.96 6.21

39.5 36.6

40.5 25.1

81.48 30.56

57.65 17.1

95.18 21.98

77.1 18.9

47.49 16.9

88.35 29.67

64.40 22.34

80.88 12.7

91.0 19.3

38.8 20.5

2.18 0.54

2.11 0.69

0.001 0.005

2.7 0.9

7.55e15 0

4.10 20.49

0.02 0.035

1.78e5 5.17e5

0.02 0.02

2.7e7 1.55e6

0.06 0.16

1.66 9.6e16

1.63e7 3.41e6

63.5 52.3

4.0e15 5.27e15

8.42 3.18

9.30 2.02

10.67 1.26

13.0 2.7

7.4 4.7

1819.35 139.54

1784.58 191.56

13.32 45.35

676.36 541.40

3.9 19.6

1837.93 128.34

1812.28 184.83

35.83 55.79

737.21 611.45

18.2 49.4

8564

W.K. Wong et al. / Expert Systems with Applications 39 (2012) 8557–8572

1050 FCPSO PSOLDIW PSOTVAC APSO CLPSO

100

b

Mean fitness

Mean fitness

a

FCPSO PSOLDIW PSOTVAC APSO CLPSO

103

102

10−50 101

10−100 0

50000

100000

100 150000

200000

0

50000

FEs

c

FCPSO PSOLDIW PSOTVAC APSO CLPSO

d

10−10

1010

FCPSO PSOLDIW PSOTVAC APSO CLPSO

0

50000

100000

150000

10−10

10−20

200000

0

50000

FEs

FCPSO PSOLDIW PSOTVAC APSO CLPSO

f

100

150000

200000

1010

FCPSO PSOLDIW PSOTVAC APSO CLPSO

100

Mean fitness

Mean fitness

100000

FEs

1010

10−10

10−20

200000

100

Mean fitness

Mean fitness

100

e

150000

FEs

1010

10−20

100000

0

50000

100000

150000

200000

10−10

10−20

0

50000

FEs

100000

150000

200000

FEs

Fig. 5. Performance of the algorithms for twelve 30-dimensional benchmark functions. (a) f1. (b) f2. (c) f3. (d) f4. (e) f5. (f) f6. (g) f7. (h) f8. (i) f9. (j) f10. (k) f11. (l) f12. (m) f13. (n) f14. (o) f15. (p) f16. (q) f17. (r) f18.

Rotated Generalized Penalized : D1 P 2 2 ðyi  1Þ2 ½1 þ 10 sin ðpyiþ1 Þ f15 ðxÞ ¼ Dp f10 sin ðpy1 Þ þ

Rotated Weierstrass : f 16 ðxÞ

! kmax kmax D X X X k k k k ½a cosð2pb ðyi þ 0:5ÞÞ  D ½a cosð2pb  0:5Þ ; ¼

i¼1

D P þðyD  1Þ2 g þ uðxi ; 10; 100; 4Þ; i¼1   where yi ¼ 1 þ 14 ðxi þ 1Þ ; uðzi ; a; k; mÞ 8 m zi > a; < kðzi  aÞ ; ¼ 0; a 6 zi 6 a; : kðzi  aÞm ; zi < a:

i¼1

ð31Þ

z¼Mx

a ¼ 0:5;

k¼0

b ¼ 3;

k¼0

kmax ¼ 20;

y ¼ M  x;

ð32Þ

Composite Function 3 in Liang et al. (2005) f17: (CF3) is composed using ten Griewank functions. Rotated Composite Function 3 in Liang et al. (2005) f18: (Rotated CF3) is composed using ten Griewank functions with rotation.

8565

W.K. Wong et al. / Expert Systems with Applications 39 (2012) 8557–8572

g

h

1010

100

105

FCPSO PSOLDIW PSOTVAC APSO CLPSO

100 PSOLDIW

Mean fitness

Mean fitness

FCPSO PSOTVAC

10−10

APSO CLPSO

10−5

10−20 10−10 10−30 10−15 0

50000

100000

150000

200000

0

50000

FEs

i

100000

150000

j

105

FCPSO PSOLDIW PSOTVAC APSO CLPSO

Mean fitness

Mean fitness

103

101 FCPSO PSOLDIW PSOTVAC APSO CLPSO

10−1

200000

FEs

102

101

100 10−3 0

50000

100000

150000

10-1

200000

0

50000

FEs

k

l

FCPSO PSOLDIW PSOTVAC APSO CLPSO

300

200000

350 FCPSO PSOLDIW PSOTVAC APSO CLPSO

300 250

Mean fitness

250

Mean fitness

150000

FEs

350

200 150

200 150

100

100

50

50

0

100000

0

0.5

1

FEs

1.5

0

2 5

0

0.5

1

FEs

x 10

1.5

2 5

x 10

Fig. 5 (continued)

Experiments were carried out to compare five PSO algorithms including the proposed FCPSO on the 18 test problems with 30 dimensions. The population size of five PSO was 20. Four existing PSO algorithms are shown in Table 2. The first PSO was LDIW (Shi & Eberhart, 1998a) with linearly decreasing inertia weight. TVAC (Ratnaweera et al., 2004) was a PSO with time-varying acceleration parameters and incorporating a self-organizing method. CLPSO delivered a comprehensive-learning strategy, which is used to yield better performance for multimodal functions (Liang et al., 2006). APSO was an adaptive PSO, which can adjust the acceleration coefficients and inertia weight adaptively (Zhan et al., 2009) according to an evolutionary factor by calculating average distance. The parameters for these PSOs were provided in Table 2.

In all the experiments, the algorithm configuration of the FCPSO was listed as follows. The inertia weight was set to 0.65 since it is the mean value of upper bound and lower bound of LDIW. In (7), c1 and c2 were set to 1.25 and 1.75 so that the mean value of c1 and c2 was 1.5, which was close to the value of c = 1.494 in CLPSO. In the following experiment, it was shown that the bounds of c2  c1 were vital for the search performance. The other coefficients in (10) and (11) were assigned as k = 3 and l = 2. Learning rates j1 and penalizing rate j2 were set to 0.1 and 0.001, respectively. All the algorithms used the same number of 2  105 fitness evaluations (FEs) for each test function. All the experiments were performed on the same machine with a Core 2 2.26-GHz CPU, 2GB memory, and Windows XP operating system. Each algorithm

8566

W.K. Wong et al. / Expert Systems with Applications 39 (2012) 8557–8572

m

n

105

FCPSO PSOLDIW PSOTVAC APSO CLPSO

5

10

10-0

Mean fitness

Mean fitness

100

10−5

FCPSO PSOLDIW

10−10

10−5

10−10

PSOTVAC APSO CLPSO

10−15 0

50000

100000

150000

200000

0

50000

FEs

150000

200000

FEs

1010

FCPSO PSOLDIW PSOTVAC APSO CLPSO

p

50 FCPSO PSOLDIW PSOTVAC APSO CLPSO

45 40

100

35

Mean fitness

Mean fitness

o

100000

30 25

10−10

20 15 10

10−20

5 0

50000

100000

150000

200000

0

0.5

FEs

r

104

FCPSO PSOLDIW PSOTVAC APSO CLPSO

102

100

10−2

Mean fitness

Mean fitness

q

1

1.5

FEs

2 x 105

104

FCPSO PSOLDIW PSOTVAC APSO CLPSO

102

100

0

50000

100000

150000

200000

10−2

FEs

0

50000

100000

150000

200000

FEs Fig. 5 (continued)

repeated 25 times independently for eliminating random discrepancy. 4.2. Comparisons on the solution accuracy The mean solutions, the best solution and standard deviation (Std. Dev.) of the solutions are listed in Table 3. The best result among those PSOs is indicated by boldface in the table. Fig. 5 depicts the comparisons in terms of convergence, mean solutions and evolution processes in solving 18 benchmark functions. It can be seen from the Table 3 and Fig. 5(a) that the FCPSO provides the best performance on the Sphere function f1, which is used

to test the convergent rates. However, on f2 and f10, CLPSO outperforms FCPSO a little and FCPSO ranks second on f2 and third on f10. On Rotated Quadric function f9, APSO outperforms other variants of PSO significantly. Although FCPSO could not deliver good performance on f9, it can be seen that FCPSO has better search performance than CLPSO. Table 3 and Fig. 5(c)-(h) illustrate the comparisons on the multimodal functions without rotation that are hard to optimize. FCPSO achieves the global optimum on the optimization of complex functions f3, f4, and f6. FCPSO performs much better than other PSOs on f3-f8, which is usually difficult to optimize. From Table 3 and Fig. 5(c)- (h), it can be found that FCPSO can offer good search performance on multimodal functions. When

8567

W.K. Wong et al. / Expert Systems with Applications 39 (2012) 8557–8572 Table 4 Convergence speed and algorithm reliability comparisons; ‘’ representing no runs reached an acceptable solution.

f1 f2 f3 f4 f5 f6 f7 f8 f9 f10 f11 f12 f13 f14 f15 f16 f17 f18

Mean FEs Right Percentage Mean FEs Right Percentage Mean FEs Right Percentage Mean FEs Right Percentage Mean FEs Right Percentage Mean FEs Right Percentage Mean FEs Right Percentage Mean FEs Right Percentage Mean FEs Right Percentage Mean FEs Right Percentage Mean FEs Right Percentage Mean FEs Right Percentage Mean FEs Right Percentage Mean FEs Right Percentage Mean FEs Right Percentage Mean FEs Right Percentage Mean FEs Right Percentage Mean FEs Right Percentage Mean Reliability

(%) (%) (%) (%) (%) (%) (%) (%) (%) (%) (%) (%) (%) (%) (%) (%) (%) (%)

LDIW

TVAC

CLPSO

APSO

FCPSO

107617 100 96049 92 104540 68 122168 100 114125 96 112570 44 94477 100 118712 72 130517 20 103001 68 108145 72 118632 80 – 0 112592 48 140204 52 92092 72 – 0 – 0 54.67

45669 100 42884.3 96 46573 84 52573.4 96 50955.1 88 47290 32 – 0 – 0 62405 64 38136 76 63237.2 100 68153.8 96 53927 4 49028 32 – 0 44725 68 – 0 – 0 52

56502 100 58801 100 56953.8 100 44623.7 100 63211.3 100 63947.3 100 25983 100 78425.5 100 – 0 47542.8 100 119620 64 92260 92 106367 96 71641 100 30075 100 139420 28 53530 96 66404 84 86.7

8953 100 10376.56 100 9379 100 4319 100 56420 100 17763.82 44 28566 100 116900 4 23066.80 100 4199.80 100 2489 92 2481 72 – 0 9318 48 11016.9 12 3010 16 8008 12 7925 12 61.7

20902.1 100 23474.4 100 7869.52 100 3962.4 100 25787.2 100 32041.9 100 4107 100 35914.6 100 101754 60 31925.3 100 37976.8 100 18336.7 100 31496.1 100 34000 100 2334 100 47280 68 16961.1 96 34037.2 88 95.1

dealing with rotated functions, FCPSO surpasses other PSOs on six multimodal rotated functions f11-f16, as seen from Table 3 and Fig. 5(k)-(p). One can observe that the proposed method can help the PSO to search for the optimum and maintain a high convergence speed when dealing with multimodal rotated functions. On composite functions f17 and f18, FCPSO can deliver the best performance among these PSOs and performs a little better than CLPSO. The capability of avoiding local optima and finding global optimum of multimodal functions indicates the superiority of FCPSO. Comparing the results and the convergence graphs, among these five PSO algorithms, CLPSO has good global search ability and slow convergence speed. APSO can converge to the best solution found so far very quickly though it is easy for APSO to stuck in the local optima. FCPSO has good local search ability and global search ability at the same time. 4.3. Comparisons on convergent rate The convergent rate for achieving the global optimum is another key point for testing algorithm performance. In solving real-world optimization problems, the ‘‘FE’’ overwhelms the algorithm overhead. Usually, the FE accounts for the most time of running algorithms as the PSO is highly computation efficient. Table 4 shows that FCPSO needs least FEs to achieve the acceptable solution on f2, f3, f4, f5, f7, f8, f9 and f10, which reveals that FCPSO has a higher convergent rate than other algorithms. Though APSO outperforms FCPSO on the other functions, FCPSO ranks second on these functions. In addition, APSO has much worse successful ratio and accuracy than FCPSO on the tested functions.

Table 5 Computational complexity of algorithms. Algorithm

T0

T1

b2 T

b 2  T 1 Þ=T 0 ðT

LDIW TVAC CLPSO APSO FCPSO

1.89 1.89 1.89 1.89 1.89

17.42 17.42 17.42 17.42 17.42

40.67 41.39 19.96 41.15 21.68

12.3 12.7 1.35 12.6 2.25

4.4. Comparisons on successful percentage Table 4 also shows that FCPSO yields the highest percentage for achieving acceptable solutions in 25 runs. According to the no free lunch theorem (Wolpert & Macready, 1997), any elevated performance over one class of problems is offset by performance over another class. Hence, one algorithm cannot perform better than all the others on every optimization problem. In summary, the FCPSO performs best on both unimodal and multimodal functions. Owing to feedback mechanism and leader local search, the FCPSO processes capabilities of fast convergence, highest successful ratio, good FEs results and best search accuracy among these PSOs. 4.5. Analysis on the computational complexity of algorithms Following the methods in Suganthan et al. (2005), the computational complexity of the algorithms discussed in this paper is comb 2 can be puted as shown in Table 5. Computations of T0, T1, and T referred to Suganthan et al. (2005) and thus one can obtain

8568

W.K. Wong et al. / Expert Systems with Applications 39 (2012) 8557–8572

Table 6 Advantages of feedback control mechanism and leader local search. Algorithms

Functions

f1

f4

f5

f6

f11

f13

FCPSO without leader local search and feedback control FCPSO without leader local search

Average Std. Dev. FEs Average Std. Dev. FEs Average Std. Dev. FEs

8.91e29 8.37e42 56502 5.25e42 4.67e42 28535 1.1e52 3.46e52 20902

0 0 44623.7 0 0 7622 0 0 3962

2.21e14 4.38e15 63211 7.54e15 3.99e15 34111 7.69e15 1.24e15 25787

3.65e14 1.74e13 63947 0 0 34098 0 0 32041

95.18 21.98 119620 48.2 11.0 28996 47.4 16.9 37976

7.55e15 0 106367 7.55e15 0 45718 7.55e15 0 31496

FCPSO

b 2  T 1 Þ=T 0 . The results of LDIW, TVAC, APSO, CLPSO and FCPSO ðT are shown in Table 5. All of the time values were measured in CPU seconds. The results of LDIW, TVAC, APSO, CLPSO and FCPSO were obtained on a PC with CPU Dual Core 2.26 GHz, RAM 2G, Platform: MATLAB 2009a. The computational complexity of FCPSO was modest, with a complexity value smaller than 3. 4.6. Performance of feedback control and leader local search Two techniques, i.e., feedback control mechanism and leader local search were also used to test on the search performance of FCPSO. Results of 25 independent runs are shown in Table 6. The results show that FCPSO without leader local search can outperform FCPSO without both feedback control mechanism and leader local search on accuracy and convergence rate. This indicates that feedback control mechanism is a useful tool for achieving good performance of PSO. After introducing leader local search, FCPSO can further refine solutions on unimodal functions and can make the swarm converge more quickly. FCPSO can deliver not only the highest accuracy on unimodal functions, but also offer a good global search performance on multimodal functions. To summarize, the full FCPSO is the most powerful for the tested functions. The results verify that feedback control method can accelerate the convergence and enhance accuracy and leader local search can make the swarm to have a better search performance. 4.7. Sensitivity of parameters setting In this subsection, the effects of c1 and c2 are investigated on FCPSO here. Suitable c1 and c2 affect the search performance of FCPSO greatly. For the sake of simplicity, in this paper, the mean value of c1 and c2 was set to 1.5, which was close to the value 1.494 adopted in CLPSO. After fixing the mean value, the range of c2  c1 was changed gradually. Appropriate c1 and c2 can enhance global and local search capabilities and reduce inefficient learning strategy leading to consumption saving of FEs. The results of mean

Table 7 Effects of c1 and c2 on search accuracy and convergence rate (The best value is in Bold font.).

f1 f3 f4 f6 f13 f14

Average Mean FEs Average Mean FEs Average Mean FEs Average Mean FEs Average Mean FEs Average Mean FEs

c1 = 1.05

c1 = 1.15

c1 = 1.25

c1 = 1.35

c1 = 1.45

5.9e50 22174 0.11 7907 0 3876 0.00002 34697 7.55e15 34381 1.7e7 36522

1.1e51 20802 0.11 8502 0 4106 0 33521 7.83e15 34685 1.7e7 33422

1.1e52 20902 0 7869 0 3962 0 32041 7.55e15 31496 2.7e7 34000

2.1e53 21062 0.04 7887 0 3955 0 31907 7.12e15 31450 1.7e7 35965

8.6e56 20316 0.07 7895 0 3825 0.0002 32349 7.55e15 32892 6.3e5 39106

values and standard deviations of the solutions are provided in Table 7. c1 is altered from 1.05 to 1.45 with an interval 0.1, i.e., the range of c2  c1 was changed from 0.9 to 0.1 with an interval 0.2. Results showed that c1 = 1.25 and c1 = 1.35 reveal the better performance than other cases. A larger c1 can make the FCPSO perform better on Sphere function f1. c1 = 1.25 performs best on f3 and this is similar to the best results on other functions. Therefore, in this paper, c1 = 1.25 is adopted, that is, the range is set to c2  c1 = 0.5 and c2 is chosen as 1.75. 5. Application to time-series prediction Particle swarm optimization is widely applied in many fields. One of the potential area is to design the suitable structure of RBFNN. A fast and efficient PSO is essential to construct RBFNN since it is a multimodal optimization. In this section, a novel RBFNN using an evolving clustering method is presented, which combines the proposed FCPSO and the K-means method. A new encoding scheme is proposed for time series prediction. By using the novel RBFNN, the hidden centers and controlling parameters are selected automatically, and the input variables are chosen at the same time. 5.1. Evolving RBFNN When using PSO to design a RBFNN, one of the most important issues is to encode the RBF network into the particle with an efficient method. When the structure (number of hidden nodes and number of inputs) and the network parameters (centers, widths, and weights) associated are specified, the RBF network is specified. The proposed encoding scheme uses a particle that consists of three parts, as shown in Fig. 6. These three parts deal with the optimization of the number of hidden nodes, optimizing the controlling parameters, and the selection of hidden nodes. The fixed length for a particle makes the evolution procedure easier, and the exact number of nodes and inputs cannot be given a priori, while the maximum number can be estimated or given. The maximum number of nodes and inputs can thus be defined in the particle such that the dimensions of the particle can be fixed. In the first part, the first dimension of the particle stands for the number of hidden nodes. Note that we define the maximum number of hidden nodes to ensure Ni 6 Ni,max. The value of Ni is used for the initial number of cluster centers for the K-means clustering method. Once the number of hidden nodes Ni is fixed, in the second part of the particle, the Ni number of controlling parameters is activated. In the third part of the particle, Mmax stands for the maximum number of input variables. The Mmax entries are positive floating-point numbers in (0, 1), each of which determines whether the corresponding input variable is to be activated (i.e. to be used for prediction). If a specified dimension in the third part is smaller than 0.5 (Ai,j < 0.5), then the corresponding input variable is not selected,

8569

W.K. Wong et al. / Expert Systems with Applications 39 (2012) 8557–8572

Fig. 6. The novel structure for optimizing RBFNN.

and vice versa (Ai,j P 0.5). Thus Ai,j behaves like control dimension in the particle governing the selection of the active input variable. From the above description, we can see that the information of particle in the second part represents whether or not input is used for the current architecture of an individual. For example, in the first part, we use L1 = 1 dimension to represent the network hidden layer nodes. If the maximum number of nodes is set as Nmax, the dimension of the second part is L2 = Nmax. If the maximum number of input variable is chosen as Mmax, the dimension of third part is L3 = Mmax. Thus, the length of a chromosome to represent an RBF network becomes D = L1 + L2 + L3 = 1 + Nmax + Mmax. The initialization of the newly encoded chromosome using PSO can be done with the standard PSOs using randomly generated numbers. In the first part, the initialization of Ni belongs to Ni 2 (1, Nmax). In the evolution process, the Ni is a real number belong to Ni 2 (0, Nmax). We set the exact hidden number as Ni,hidden = floor (Ni) + 1. If floor (Ni) + 1 > Nmax, we set Ni,hidden = Nmax. For the third part, in the extreme case that all the dimensions are

generated as zeros, some dimensions should be enforced to be a real number belonging to [0.5, 1] to avoid the situation where no inputs can be selected and the algorithm cannot be executed. The encoded particle will be evolved by the conventional evolutionary operations. The summation of the Ai,j P 0.5 in the third part gives the actual number of inputs to be used, and the position of Ai,j P 0.5 indicates the corresponding input to be selected. As the evolutionary operations alter the value of each dimension in the particle, with the proposed encoding scheme, the flexible architecture has been presented. Since the architecture and the network parameters evolve at the same time, the structure and parameters selection problem for the RBFNN can also be performed simultaneously. To summarize, the mechanism make it possible to find a RBFNN that is composed of a compact number of nodes and inputs with optimal accuracy. Combining the newly proposed encoding scheme into constructing RBFNN, the algorithm for evolving RBFNN based on a PSO, as shown in Fig. 7, is described in the following steps: 1. Encode all the parameters into the particle using the proposed encoding scheme. 2. Initialize parameters of each particle. Randomize an initial population of NP individual particle. Each particle represents a set of parameters. 3. Preprocess data using the inputs selected by the third part of each particle. 4. Use the K-means method to select hidden nodes by the first part of the each particle. 5. Obtain the weight matrix W using OLSA and the controlling parameters obtained by the second part of the particle. 6. Calculate objective functions. After the parameters are generated, the objective function of each particle is calculated. The mean square error (MSE) for the training data or the testing data is regarded as the objective function of each particle accounting for the accuracy requirement. In order to obtain the compact network, i.e., the network with fewer network parameters, the number of nodes and inputs should also be included in the objective function as well. The overall objective function for each particle is calculated as

Ji ¼

Q 1 X ^i;j Þ2 ; ðy  y Q j¼1 i;j

i ¼ 1; . . . ; S;

j ¼ 1; . . . ; Q;

ð33Þ

^i;j is the prewhere Q is the data length; yi,j is the target output and y dicted output; 7. Update the velocity and positions of particles. Steps 3–6 correspond to one generation. The evolution process repeats for a fixed number of generations or will cease when the search process converges with a given accuracy. The best particle in the swarm will be used to determine the optimal network parameters.

Fig. 7. The novel structure for optimizing RBFNN.

One of the most important concerns in constructing RBFNN is to select the appropriate number and the width of hidden nodes. In this paper, we use one dimension of the particle to represent the number of hidden nodes, while the K-means method is employed to select width of the centers. In the evolution process, the number and width of hidden nodes varies. In the end of the evolution, we

8570

W.K. Wong et al. / Expert Systems with Applications 39 (2012) 8557–8572

Table 8 Comparison results for Box–Jenkins time series prediction with input selection. Refs.

Inputs

Number of rules or nodes

MSE of testing

Kim et al. (2006) Sugeno and Yasukawa (1993) Lin and Cunningham (1995)

y(t), u(t  3), y(t  1), u(t  4), u(t  3), y(t  2), y(t  1), u(t  6), u(t  5), u(t  3), y(t  1), u(t  5), u(t  4), u(t  3) y(t  3), y(t  2), y(t  1), u(t  6), u(t  5), u(t  4) y(t  2), y(t  1), u(t  6), u(t  5), u(t  2), u(t  1)

2 4 25

0.129 0.19 0.261

45

0.169

94

0.065

30

0.0582

Nie (1995) Du and Zhang (2008) our method

can find the appropriate number and width of the hidden nodes. On the other hand, one limitation of the K-means method is the initial clustering number. In this paper, the initial clustering number is evolved by FCPSO, which can be automatically selected. Thus, the centers of hidden nodes can be automatically selected by the Kmeans method. The dimension of this encoding scheme is cut down by nearly half, compared with the encoding method in Du and Zhang (2008). Recently, some evolving K-means methods have been proposed in the literature. It is worth mentioning that the clustering method can easily be extended using the first part of the particle depicted in Fig. 6 and setting an appropriate fitness function.It is worth mentioning that the range of the first part of the particle is in [1,Nmax], which is different from other parts of the particle. In order to tune the number of hidden nodes rapidly,

a

62 prediction data source data

60

the maximum velocity in the first part is different from other parts. In the experiments, the Vmax in the first part is set to 20% of the search range, i.e. Vmax = 0.2 ⁄ Nmax. 5.2. Application to Box–Jenkins gas furnace data In this subsection, the presented evolving RBF networks are applied to a set of time-series prediction problems: Box–Jenkins (Hyndman, xxxx). The presented results are compared with the results reported in the recent studies. In the experiment, the population size is 10, the generation step for PSO is 50 and the generation step for K-means is 50. The parameters setting for PSO are the same as in Section 4. The Box–Jenkins furnace data was recorded from a combustion process of a methane-mixture. The data set

b

1.5

1 58 0.5

54

Error

x (t)

56

52 50

0

−0.5

48 −1 46 44

−1.5 0

50

100

150

200

250

0

50

100

t

150

200

250

t

c 25

times

20

15

10

5

0 y(t−1) y(t−2) y(t−3) y(t−4) u(t−1) u(t−2) u(t−3) u(t−4) u(t−5) u(t−6)

ith input variable Fig. 8. Prediction results for Box–Jenkins time series with input selection. (a) prediction data and real data; (b) prediction error; (c) input selection in 25 runs.

W.K. Wong et al. / Expert Systems with Applications 39 (2012) 8557–8572

originally consisted of 296 data points [y(t), u(t)]. The input u(t) was the gas flow rate into the furnace and the output y(t) was carbondioxide (CO2) concentration in the outlet gas. In this experiment, we are trying to predict y(t) based on

½yðt  1Þ; yðt  2Þ; yðt  3Þ; yðt  4Þ; uðt  1Þ; uðt  2Þ; uðt  3Þ; uðt  4Þ; uðt  5Þ; uðt  6Þ:

ð34Þ

Consequently, the effective number of data points is reduced to 290, providing 145 for training and 145 for testing. The maximum number of inputs was set to 10, and the input candidates are the same as shown in (34); the maximum number of nodes is 40. The best result within 25 runs is listed in Table 8 and the prediction results are shown in Fig. 8. The average number of hidden centers is 30. It can be seen from Table 8 and Fig. 8 that the best prediction accuracy can be obtained while keeping the numbers of nodes and inputs in an acceptable range. Compared with the recent results presented in Table 8, it can be seen that with appropriate number of nodes and inputs, the prediction error is less. Although the prediction results in Du and Zhang (2008) are close to the results of our method, the hidden nodes are nearly three times than ours and the input variables required are much more than the method used in this study. 6. Conclusion An algorithm termed FCPSO is proposed for solving optimization problems in this paper. A novel PSO employing feedback control mechanism is presented where the fitness information is used to define an evolutionary fitness function and thus determine the parameters in PSO. In such a way, each particle has been assigned unique parameters such as the acceleration coefficient, learning probability, number of potential exemplar and refreshing gap. The new strategy results in the particles having parameters adaptively according to the search environment. A leader local search method is also proposed to refine the search solution. From experiments, it can be observed that this learning strategy enables the FCPSO to make use of the information in swarm more effectively to generate better quality solutions frequently, in comparison with four PSO variants. Based on the results of the five algorithms on the 18 test problems belonging to four classes, it is confirmed that FCPSO improves the PSO’s performance significantly and offers the best performance on these problems irrespective of whether they are rotated or not. It is worth mentioning that another attractive property of the FCPSO is that it does not introduce any complex calculations to the original simple PSO. Hence, the FCPSO is simple and easy to implement like the original PSO. The proposed FCPSO is applied to constructing a RBFNN, together with K-means method for time-series prediction, where the inputs selection is considered. The experiment performed on Box–Jenkins gas furnace data shows the superiority of the proposed scheme. The proposed fitness feedback method can be used not only in CLPSO, but also other kinds of PSOs. References Andrews, P. S. (2006). An investigation into mutation operators for particle swarm optimization. In Proceedings of IEEE congress on evolutionary computations, Vancouver, BC, Canada (pp. 1044–1051). Angeline, P. J. (1998). Using selection to improve particle swarm optimization. In Proceedings of IEEE congress on evolutionary computation, Anchorage, AK (pp. 84– 89). Bergh, F. VD., & Engelbrecht, A. P. (2004). A cooperative approach to particle swarm optimization. IEEE Transactions on Evolutionary Computation, 8, 225–239. Buchtala, O., Klimek, M., & Sick, B. (2005). Evolutionary optimization of radial basis function classifiers for data mining applications. IEEE Transactions on System Man and Cybernetics, 35, 928–947. Cai, X., Cui, Z., Zeng, J., & Tana, Y. (2008). Dispersed particle swarm optimization. Information Processing Letters, 105, 231–235.

8571

Chen, C., Wu, Y., & Luk, B. L. (1999). Combined genetic algorithm optimization and regularized orthogonal least squares learning for radial basis function network. IEEE Transactions on Neural Networks, 10, 1239–1243. Chen, Y. P., Peng, W. C., & Jian, M. C. (2007). Particle swarm optimization with recombination and dynamic linkage discovery. IEEE Transactions on System Man and Cybernetics B, 37, 1460–1470. Clerc, M., & Kennedy, J. (2002). The particle swarm: Explosion, stability, and convergence in a multi-dimensional complex space. IEEE Transactions on Evolutionary Computation, 6, 58–73. Coello, C. A. C., Pulidos, G. T., & Lechuga, M. S. (2004). Handling multiple objectives with particle swarm optimization. IEEE Transactions on Evolutionary Computation, 8, 256–279. Du, H., & Zhang, N. (2008). Time series prediction using evolving radial basis function networks with new encoding scheme. Neurocomputing, 71, 1388–1400. Eberhart, R.C., and Shi, Y.H., 2001. Particle swarm optimization: Developments, applications and resources. Proc. IEEE Congr. Evol. Comput. Seoul, Korea, pages 204–210. Guyon, I., & Elisseeff, A. (2003). An introduction to variable and feature selection. Journal of Machine Learning Research, 3, 1157–1182. Harpham, C., & Dawson, C. W. (2006). The effect of different basis functions on a radial basis function network for time series prediction: A comparative study. Neurocomputing, 69, 2161–2170. Hyndman. R.J., Time series data library. http://www-personal.buseco.monash edu. au/hyndman/TSDL/. Kennedy, J., 1997. Dynamic multiple swarms in multiobjective particle swarm optimization. Proc. IEEE Int. Conf. Evol. Comput., Indianapolis, IN, pages 303– 308. Kennedy, J., 1999. Small worlds and mega-minds: Effects of neighborhood topology on particle swarm. Proceedings of the 1999 Congress on Evolutionary Computation, pages 1931–1938. Kennedy, J., and Eberhart, R., 1995. Particle swarm optimization. Proceedings of IEEE International Conference On Neural Network, pages 1942–1948. Kennedy, J., and Mendes, R., 2002. Population structure and particle swarm performance. Proc. IEEE Congr. Evol. Comput., Honolulu, HI, pages 1671–1676. Kennedy, J., & Mendes, R. (2006). Neighborhood topologies in fully informed and best-of-neighborhood particle swarms. IEEE Transactions on System Man and Cybernetics C Applied Revision, 36, 515–519. Kim, M. S., Kim, C. H., & Lee, J. J. (2006). Evolving compact and interpretable Takagi– Sugeno fuzzy models with a new encoding scheme. IEEE Transactions on System Man and Cybernetics, 36, 1006–1023. Liang, J. J., Suganthan, P. N., & Deb, K. (2005). Novel composition test functions for numerical global optimization. Proceedings of Swarm Intellegince Symposium, 82–102. Liang, J. J., Suganthan, P. N., & Baskar, S. (2006). Comprehensive learning particle swarm optimizer for global optimization of multimodal functions. IEEE Transactions on Evolutionary Computation, 10, 281–295. Lin, Y., & Cunningham, G. A. (1995). A new approach to fuzzy-neural system modeling. IEEE Transactions on Fuzzy Systems, 3, 190–197. Lovbjerg, M., Rasmussen, T. K., & Krink, T. (2001). Hybrid particle swarm optimizer with breeding and subpopulations. In Proceedings of Genetic Evolutionary Computation Conference (pp. 469–476). Makridakis, S., Wheelwright, S., & Hyndman, R. (1998). Forecasting: Methods and Applications. New York: Wiley. Mendes, R., Kennedy, J., & Neves, J. (2004). The fully informed particle swarm: Simpler, maybe better. IEEE Transactions on Evolutionary Computation, 8, 204–210. Nie, J. (1995). Constructing fuzzy model by self-organizing counter propagation network. IEEE Transactions on System Man and Cybernetics, 25, 963–970. Krohling, R. A., & dos Santos Coelho, L. (2006). Coevolutionary particle swarm optimization using gaussian distribution for solving constrained optimization problems. IEEE Transactions on System Man and Cybernetics B, 36, 1407–1416. Ratnaweera, A., Halgamure, S., & Watson, H. (2004). Self-organizing hierarchical particle swarm optimizer with time-varying acceleration coefficients. IEEE Transactions on Evolutionary Computation, 8, 240–255. Salomon, R. (1996). Reevaluating genetic algorithm performance under coordinate rotation of benchmark functions. BioSystems, 39, 263–278. Shi, Y., & Eberhart, R. C. (1998). A modified particle swarm optimizer. In Proceedings of IEEE congress evolutionary computation (pp 69–73). Shi, Y., & Eberhart, R. C. (1998b). Parameter selection in particle swarm optimization. In Proceedings of the 7th international conference on evolutionary programming VII. LNCS (1447, pp. 591–600). New York: Springer-Verlag. Suganthan, P. N., Hansen, N., Liang, J. J., Deb, K., Chen, Y.-P., Auger, A., et al. (2005). Problem definitions and evaluation criteria for the cec2005 special session on real-parameter optimization. Proceedings IEEE Congress Evolutionary Computation, 1–50. Sugeno, M., & Yasukawa, T. (1993). A fuzzy-logic-based approach to qualitative modeling. IEEE Transactions on Fuzzy Systems, 1, 7–31. Tang, Y., Fang, J., Xia, M., & Yu, D. (2009). Delay-distribution-dependent stability of stochastic discrete-time neural networks with randomly mixed time-varying delays. Neurocomputing, 72, 3830–3838. Tang, Y., Wang, Z., & Fang, J. (2011). Controller design for synchronization of an array of delayed neural networks using a controllable probabilistic pso. Information Sciences, 181, 4715–4732. Wang, Z., Ho, D. W. C., & Liu, X. (2003). Variance-constrained filtering for uncertain stochastic systems with missing measurements. IEEE Transactions on Automatic Control, 48, 1254–1258.

8572

W.K. Wong et al. / Expert Systems with Applications 39 (2012) 8557–8572

Wang, Z., Yang, F., Ho, Daniel W. C., & Liu, X. (2007). Robust h-infinity control for networked systems with random packet losses. IEEE Transactions on Systems Man and Cybernetics – Part B, 37, 916–924. Wolpert, D. H., & Macready, W. G. (1997). No free lunch theorems for optimization. IEEE Congress on Evolutionary Computation, 1, 67–82. Wong, W. K., Yuen, C. W. M., Fan, D. D., Chan, L. K., & Fung, E. H. K. (2009). Stitching defect detection and classification using wavelet transform and bp neural network. Expert Systems with Applications, 36, 3845–3856. Wong, W. K., Zeng, X. H., & Au, W. M. R. (2009). A decision support tool for apparel coordination through integrating the knowledge-based attribute evaluation expert system and the T–S fuzzy neural network. Expert Systems with Applications, 36, 2377–2390.

Yao, X., Liu, Y., & Lin, G. M. (1999). Evolutionary programming made faster. IEEE Transactions on Evolutionary Computation, 3, 82–102. Yen, G. G., & Leong, W. F. (2009). Dynamic multiple swarms in multiobjective particle swarm optimization. IEEE Transactions on System Man Cybernetics – Part A: Systems and Humans, 39, 890–911. Zhan, Z., Zhang, J., Li, Y., & Chung, H. S. H. (2009). Adaptive particle swarm optimization. IEEE Transactions on System Man and Cybernetics-B, 39, 1362–1381. Zhang, G. P. (2003). Time series forecasting using a hybrid arima and neural network model. Neurocomputing, 50, 159–175.