2.1. PSO with mutation

In the PSO, each single solution is a “bird” in the search space. The particles fly through the problem space by following the current optimum particles. In every iteration, each particle is updated by following two “best” values. The first one is the best solution it has achieved so far, it is a personal best position denote by pib . Another “best” value, which is tracked by the particle swarm optimizer, is the best value obtained so far by any particle in the population, it is a global best position and defined as pgb . All of particles have fitness values which are evaluated by the fitness function to be optimized. Each particle updates according to those two best values, and then a new generation of population is created.

Suppose that the searching space is D dimensional with m randomly initialized particles in it, the particle swarm can be indicated by following parameters: xi=(xi1,xi2,...,xiD) stands for the location of particle i in the D dimensional space and it is also regarded as a potential solution, vi=(vi1,vi2,...,viD) stands for the flight velocity of particle i in the D dimensional space, pi=(pi1,pi2,...,piD) stands for the personal best position of particle i, pi=(pg1,pg2,...,pgD) stands for the global best position in the whole swarm. The particle updates its velocity and position with the following rules:

in which, i = 1,2 …,m; d = 1,2,…,D; k is iteration number; c₁, c₂ are called acceleration coefficients, which are used to adjust the maximum flight step of personal best value and global best value, rand() returns a random number between (0,1); ω is inertia weight which affects the balance of global search ability and local search ability.

In Zhong et al. in 2007 a large number of experiments proved that once ω decreases linearly with the iteration, the convergence of algorithm would be significantly improved. Therefore here we let ω=ωmax−(ωmax−ωmin)*kK , where k is the current iteration number, K is the maximum iteration number, ωmax is the maximum inertia weight, ωmin is the minimum inertia weight. The basic principles of Eq. (1) is that the velocity achieves information from the original velocity, personal best value and global best value, the number of information depends on ω , c₁ and c₂. The first part of Eq. (1) is called memory term, which denotes the impact of velocity and direction of previous iteration. The second part (the distance between current position of particle i and the personal best position) is called self-awareness term, which denotes the information that comes from its own experience. The third part (the distance between current position of particle i and the global best position) is called population-awareness term, which denotes the information that comes from another particles of the whole swarm, which reflects the knowledge sharing and cooperation. The PSO algorithm can be implemented in the following 6 steps:

Initialize generation and all particles, viz. set the initial position X of each particle and the initial velocity V randomly.

For each particle, calculate the fitness value.

For each particle, if the fitness value is better than the best fitness value f(pib) in history, set current value as the new pib

Choose the particle with the best fitness value of all the particles as pgb .

For each particle, calculate the particle velocity according to Eq. (1), and update particle position according to Eq. (2).

If the maximum iteration or the minimum error criteria is not attained, return to Step 2; otherwise end the iteration.

The PSO has been successfully applied in many continuous optimization problems. The Traveling Salesman Problem (TSP) is a typical discrete combinatorial problem. If one wants to solve the TSP with PSO, some improvements of basic PSO must be done. In Huang et al., in 2003 the concept of swap operator and swap sequence were introduced for solving the TSP. Suppose that the solution sequence of the TSP with n nodes is S=(ai),i=1,...,n. The definition of swap operator SO(i1,i2) is the points ai1 and ai2 in the solution sequence S. Swap sequence is an orderly sequence with one or more swap operators, meanwhile the order between swap operators is meaningful. Different swap sequence operate on the same solution may generate the same new solutions, the equivalent set of swap sequence is the set of swap sequence which has the same effect. Among all the equivalent sets of swap sequence, the swap sequence with least swap operators is called basic swap sequence. An array with N cities denotes the particle’s position X, All the possible arrays compose the state space of the problem.

Based on vectors, functions and operations defined above, the traditional updating equations will be changed in the following versions (Huang et al., 2003):

in which, α,β(α,β)∈[0,1] are random numbers. α(Pid−Xid) denotes that all the swap operators in the basic swap sequence (Pid−Xid) are reserved with a probability of α , similarly, β(Pgd−Xid) denotes that all the swap operators in the basic swap sequence (Pgd−Xid) are reserved with a probability of β . Thus the greater the value of α is, the more swap operators that (Pid−Xid) will be reserved, and the greater the impact of Pid is. Similarly, the greater the value of β is, the more swap operators that (Pgd−Xid) will reserve, and the greater the impact of Pgd is. The definition of operator “ ⊕ ” is the merger operator of two swap sequences. Operator “+” denotes the implementation of swap operation, operator “-” denotes to obtain the basic swap sequence of two sequences. For example：A = (1 2 3 4 5), B = (2 3 1 5 4), as can be seen that A(1) = B(3) = 1, so the first swap operator is SO(1,3), B1 = B + SO (1,3), so one gets that B1: (1 3 2 5 4), A(2) = B1(3) = 1, so the second swap operator is SO(2,3), B = B1 + SO(2,3), so one gets that B2: (1 2 3 5 4). Similarly, the third swap operator is SO(4,5), B3 = B2 + SO(4,5) = A, thus one gets a basic swap sequence: SS = A – B = (SO(1,3), SO(2,3), SO(4,5)).

The steps of the PSO algorithm for solving the TSP can be described as follows:

Step 1.Initialize generation and all the particles, set each particle a random initial solution and a random swap sequence.

If the maximum iteration or the minimum error criteria is met, turn to Step 5.

According to the particle’s current position Xidk , calculate the next position Xidk+1 , namely the new solution.

Calculate the difference between Pid and Xid , A=Pid−Xid , in which A is a basic swap sequence.

Calculate the difference between Pgd and Xid , B=Pgd−Xid , in which B is a basic swap sequence.

Calculate the velocity vidk+1 according to Eq. (3), and convert the swap sequence vidk+1 to a basic swap sequence.

Calculate the new solution xibk+1 according to Eq. (4).

Choose the particle with the best fitness value of all the particles as Pgd , turn to Step 2.

Show the result that obtained.

In the settlement of solving TSP, basic PSO generates new individual through Eq. (3) and Eq.(4), from which one can see that a basic swap sequence generated by Eq. (3) is in fact equivalent to the swap operator to a route, but the route between the two swap cities do not change, so it is easy to generate cross route which is illegal solution. To deal with the problem, inspired by the mutation operator in evolutionary algorithm, we add the mutation operator to the PSO. The specific approach is: after generating a new route by basic PSO approach during each iteration one does the mutation operator to the new route. More specific, change Step 4 as follows:

Step 4: Generate two mutate cities randomly, then reverse the order of all the cities between two mutate cities. If the length of new route is less than the original route, set the new route as xibk+1 . Otherwise, maintain the original route unchanged.

2.2. A hybrid algorithm of PSO and SA

The SA algorithm derived from the principle of solid annealing. Firstly, heat the solid to a sufficiently high temperature, and then cool it slowly. This process is based on an analogy from thermodynamics where a system is slowly cooled in order to achieve its lowest energy state. According to Metropolis criteria, the probability of particles balance at temperature T is e−ΔE/(kT) , where E is the internal energy at temperature T; ΔE is the increment of internal energy; k is the Boltzmann constant. Once one converts the internal energy E to the objective function value, and the temperature T to control parameter t, the SA algorithm of solving combinatorial optimization problems may be obtained with the initial solution i and initial control parameter t by repeating the iteration of “generate new solution → calculate the difference of objective function → accept or discard” to the current solution, and gradually reduce the value of control parameter t. The current solution is the approximate to the optimal solution when algorithm is terminated. It is a stochastic heuristic search process based on Monte Carlo iterative method. The process is controlled by Cooling Schedule, which includes initial control parameter t, attenuation factor Δt , iteration number for each t and the termination condition.

The SA algorithm used to solve the TSP can be described as follows:

Solution spaces：Solution spaces are all the routes of visiting each city once. The solution can be denoted as {ω1ω2...ωn} . ω1 , …, ωn in an array that is composed of 1 to n, which denotes one walks to start from the city ω1 , and visits along with the cities ω2 ,…, ωn orderly, then returns to the city ω1 .

Objective function: Objective function is the total distance length of the route pass through all the cities. The objective function value of the optimal route is the least one. Objective function is also called fitness function.

Criteria of new solution generation and acceptance: Here we use reverse operator to generate new solution, more specifically, choose two different number k and m between 1 to n randomly, moreover, k is smaller than m, then one swaps the cities between k to m, that is to convert {ω1ω2...ωk,ωk+1,...,ωm,...,ωn} into {ω1ω2...ωm,ωm+1,...,ωk+1,ωk,...,ωn} . Once the new route is generated, calculate the difference of distance length between new route and current route, if the length of new route is smaller, that is, Δf=f(xj)−f(xi)≤0 , set the new route as xibk+1 , if the length of new route is bigger, but exp(−Δf/t)random(0,1) , still set the new route as xibk+1 , otherwise maintain the current route unchanged.

The SA algorithm in hybrid algorithm of PSO and SA proposed calculates alternately by two steps:

1. Generate a new solution by stochastic perturbation and calculate the change of the objective function.

Decide whether the new solution is accepted or not. In the case at a high temperature, the solution that increases the objective function may be accepted by means of decreasing the temperature slowly, which may avoid to trap in local minima. In such a way the algorithm can converge to the global best solution.

The nature of basic PSO is the use of individual and global maximum to guide the position of the next iteration, which can converge fast at the early iteration. But, after several iterations current solution, personal best value and global best value tend to the same, which leads to the loss of population diversity, and makes the solution be trapped in local minima. In order to avoid this phenomenon, inspired by the SA algorithm, we redesign the algorithm’s framework: when basic PSO converges to a solution pg , use the solution pg as the initial solution of SA, accept the new solution in accordance with the Metropolis criteria. If there is such a solution y satisfies f(y)f(pg) , that is, the solution calculated by basic PSO is not the global optimal solution. Then a new solution y can be used to randomly replace a particle in the swarm, and the evolution of PSO continues, which can increase the population diversity as well as retain the previous operation procedure. If there is not such a solution y, then let f(y)f(pg) , that is, no better solution than pg has been found until the convergence of SA, which indicates that pg is the global optimal solution.

2.3. The C-TSP application and results analysis

The TSP is a well-known combinatorial optimization problem with typical NP-hard and is often used to verify the superiority of some intelligent heuristic algorithm. The mathematical description of the TSP is: Given a list of n cities in order of visiting as X={x1,x2,...,xn} and xn+1=x1 , the task is to find the shortest possible tour distance of min∑i=1,x∈Ωndxixi+1 that visits each city exactly once. The TSP can be modeled as a graph: the graph’s vertices correspond to cities and the graph’s edges correspond to connections between cities, the length of an edge is the corresponding connection’s distance. A TSP tour is now a Hamiltonian cycle in the graph, and an optimal TSP tour is the shortest Hamiltonian cycle.

Chinese Traveling Salesman Problem (C-TSP) is a typical symmetric TSP problem. Its simple description is: Given a list of 31 Chinese capital cities and their pairwise distances, the task is to find the shortest possible tour that visits each city exactly once. The C-TSP is a medium-scale TSP problem, which has (31−1)!/2=1.326*1032 possible routes.

The algorithms of solving TSP problem are divided into two categories: Exact algorithms and approximation algorithms. The exact algorithms used frequently includes branch and bound algorithms, linear programming and exhaustion method, etc.. The running time for this approach lies within a polynomial factor of O(n!), the factorial of the number of cities, so this solution becomes impractical even for only 20 cities. Approximation algorithm is divided into tour route construction algorithm, tour route optimization algorithm and heuristic algorithm, in which heuristic algorithm is the most spectacular including genetic algorithm, ant colony algorithm, particle swarm optimization, simulated annealing, differential evolution, etc.. Currently, to C-TSP problem, simulated annealing, improved genetic algorithm, differential evolution all achieve the optimal solution of 15404 km. The SA has the advantages of high, quality, robust, easy to achieve but with slow convergence.

The hybrid algorithm proposed is applied into the C-TSP. At the same time, basic PSO, SA and PSO with mutation are also applied to do the comparisons. The programming language is Matlab 7.0, and it runs on Win XP with Intel Core2 Duo 2.10 GHz CPU. We use the discrete PSO (DPSO) that proposed by Huang et al. in 2003 as the basic PSO, parameter settings are as follows: Particle number m = 200, maximum inertia weight wmax = 0.99, minimum inertia weight wmin = 0.09, acceleration coefficient c1 = 0.8, c2 = 0.4, iteration number k = 2000. The parameters of mutation in the PSO are the same as in DPSO. The parameter settings of SA are as follows: Initial temperature T0 = 5000, termination temperature Tf = 1, cycle constant L = 31000, attenuation factor α = 0.99. The parameter settings of the hybrid algorithm are as follows：Particle number m = 200, maximum inertia weight wmax = 0.99, minimum inertia weight wmin = 0.09, acceleration coefficients c1 = 0.8, c2 = 0.4, initial temperature T0 = 5000, termination temperature Tf = 1, cycle constant L = 31000, attenuation factor α = 0.95. Each algorithm has been run for 20 times, the numerical results of four algorithms are listed in Table 1, in which “Worst” denotes the worst solution in 20 runs, “Best” represents the best solution in 20 runs, “Average ” is the average fitness in 20 runs, “Optimal rate” denotes the rate that gets the times of the optimal solution (15404) over 20 times runs.

One can see from Table 1 that the result obtained by DPSO is not satisfactory. It is unable to find the optimal solution of 15404, and its average value of solutions is also away from the optimal solution, which is because current solution Xid , personal best value pib and global best value pgb tend to the same after several iterations. The DPSO is difficult to get new effective edge and new effective route, and it is easy to be trapped in local minima. On the other hand, the SA, PSO with mutation and the hybrid algorithm can obtain the optimal solution. The average and worst distances of hybrid algorithm proposed are 15453.4 and 15587, respectively, which are the smallest in all those of the four algorithms. The solutions obtained by the hybrid algorithm proposed are the best and its optimal solution rate is 20%, which is the highest. In order to further compare SA, PSO with mutation and the hybrid algorithm, we have studied the three algorithms in the fitness curve when finding the optimal solution. The results are shown in Figures 1-3, in which the x axes is the iteration number in unite of times, and the y axes is the global best route distance which is just the solution of C-TSP.

Figure 1.

The C-TSP optimization procedure of SA

From Figure 1 it is clear that SA converges to the optimal solution in about 500 iterations. Figure 2 indicates that PSO with mutation algorithm has oscillation at the early iteration, and it converges to the optimal solution in about 1000 iterations. From Figure 3 one can see that the hybrid algorithm converges to the optimal solution in about 100 iterations. The

Figure 2.

The C-TSP optimization procedure of PSO with mutation

Figure 3.

The C-TSP optimization procedure of hybrid algorithm proposed

Figure 4.

Optimal route of C-TSP

problem of local optimal solution is not only been solved, but also the convergence speed is greatly decreased. Figure 4 is the optimal route made by the hybrid algorithm for C-TSP problem. The optimal route is: Beijing - Harbin - Changchun - Shenyang - Tianjin - Jinan - Hefei - Nanjing - Shanghai - Hangzhou - Taipei - Fuzhou - Nanchang - Wuhan - Changsha - Guangzhou - Haikou - Nanning - Guiyang - Kunming - Chengdu - Lhasa - Urumchi - Xining - Lanzhou - Yinchuan - Xian - Zhengzhou - Shijiazhuang - Taiyuan - Hohhot – Beijing. The total distance length of the optimal route is 15404 km.

3.1. Ant colony optimization algorithms

Artificial ants of the ACO algorithms build solutions by performing random walk which use a certain stochastic rules on a completely connected graph GC=(C,L) whose nodes are the components C, and the set L fully connect the components C. Each connection of the map GC=(C,L) can be associated pheromone trail τij , and heuristic information ηij (The subscripts i and j are labeling of the nodes on the map).

It is important to note that artificial ants are in parallel movement independently. Although each ant is complex enough to find a (probably poor) solution to the problem under consideration, good-quality solutions can only emerge as the result of the collective interaction among the ants. This collaborative interaction is obtained via indirect communication mediated by the information ants read or write in the variables storing pheromone trail values. To some extent, this is a distributed learning process in which the single ant is not self-adaptive but, on the contrary, it can modify adaptively the way represented and perceived by other ants. Informally an ACO algorithm can be imagined as the interplay of three procedures (Dorigo & Stützle, 2004): Construction of Ants Solutions, Pheromone updating, and Daemon Actions.

1.Construction of Ants Solutions manages a colony of ants that concurrently and asynchronously visit adjacent states of the problem considered by moving through neighbor nodes of the problem’s construction graph G _C. They move by applying a stochastic local decision policy that makes use of pheromone trails and heuristic information. In such a way, ants incrementally build solutions of optimization problem. Once an ant has built a solution, the ant evaluates the solution that will be used by the Pheromone updating procedure to decide how much pheromone to deposit.

Pheromone updating is the procedure in which the pheromone trails are modified. The trails’ values move either increase, as ants deposit pheromone on the components or connections they use, or decrease due to pheromone evaporation. From a practical point of view, the deposit of new pheromone increases the probability whose components/ connections are used by either many ants or at least one ant. A very good solution produced will be used again in future ants. Differently, pheromone evaporation carries out a forgetting factor in order to avoid a too rapid convergence to a sub-optimal region, so pheromone evaporation has the advantage of generating new search areas.

Daemon Actions procedure is used to centralize the actions which can not be performed by single ants. Examples of daemon actions are the activation of a local optimization procedure, or the collection of global information used to decide whether it is useful or not to deposit additional pheromone to update the search process.

These three procedures should interact and take into account the characteristics of the problem considered. The TSP can be represented by a complete weighed graph GC=(C,L) with C being the set of nodes representing the cities, and L being the set of arcs. Each arc (i,j)∈L is assigned a value dij , which is the distance between cities i and j. In the symmetric TSP, dij=dji holds for all the arcs in L; but in the general case of the asymmetric TSP, the distance between a pair of nodes i, j is dependent on the direction of traversing the arc, that is, there is at least one arc (i, j) for which dij≠dji . More formally, TSP is described as: A finite set C{c1,c2,⋯,cN} of components is given, where N is the number of components. Set L={lij|ci,cj∈C} fully connects the components C. dij(i,j=1,2,⋯n) is the Euclid distance of arc lij :

In the TSP, G=(C,L) is a directed graph and the goal is to find a minimum length Hamiltonian circuit of the graph, where a Hamiltonian circuit is a closed path visiting each of the n nodes of G exactly once. In another way, an optimal solution to the TSP is a permutation R of the node indices {c1,c2,⋯,cn} such that the length F(R) is minimal, where F(R) is given by

The ACO can be applied to the TSP in a straightforward way. The construction graph GC=(C,L) , where the set L fully connects the components C, is identical to the problem graph; the set of states of the problem corresponds to the set of all possible partial tours; and the constraints Ω enforce that the ants construct only feasible tours that correspond to permutations of the city indexes. In all available ACO algorithms for the TSP, the pheromone trails are associated with arcs, so the τij refers to the desirability of visiting city j directly after city i. The heuristic information is chosen as ηij=1/dij thus the heuristic desirability of going from city i directly to city j is inversely proportional to the distance between the two cities. If there is the length of arc (i, j) equal to 0, then put the ηij set to a very small value. For implementation purposes, pheromone trails are usually collected into a pheromone matrix whose elements are the arcs’ pheromone trails. This can be done analogously for the heuristic information. Tours are constructed by applying the following simple construction procedure to each ant (Dorigo & Stützle, 2004): (1) choose a initial city according to some criterion. (2) make use of pheromone and heuristic values to probabilistically construct a tour by iteratively adding cities, to which the ant has not visited yet, until all cities have been visited; (3) go back to the initial city. After all ants have completed their tours, they may deposit pheromone on the tours they have followed. Sometimes, adding Daemon Actions such as the local search will improve the performance of algorithm.

3.2. Ant System (AS)

Ant System is created by Marco Dorigo and others in 1991 (Dorigo & Stützle, 2004, Dorido et al., 1991, Colorini et al., 1992, Dorigo 1992), and it is the first algorithm which is in line with the ACO algorithm framework. Initially three different versions of AS were proposed which are called ant-density, ant-quantity, and ant-cycle. These three versions are different on pheromone updating. Whereas in the ant-density and ant–quantity versions the ants update the pheromone directly after a move from one city to an adjacent city. In the ant-cycle version the pheromone deposited by each ant is set to be a function of the tour quality. The version of ant-cycle considers the quality of complete solution and uses pheromone updating method which has overall mechanism. Ant-cycle is better than the other two versions which just consider the single-step path information. Nowadays, when referring to AS, one actually refers to ant-cycle (AS described in this chapter also refers the ant-cycle version). The principle of AS is introduced as follows.

3.3. Elitist Ant System (EAS)

Elitist Ant System is the first improvement on the initial AS, which is called elitist strategy for Ant system (EAS) (Dorigo & Stützle, 2004, Dorido et al., 1991, Colorini et al., 1992, Dorigo, 1992). EAS gives the ant which has constructed so far the best path solution the elite logo, sets the so far the best path Tbs , and provides strong additional reinforcement to the arcs that is belong to the best tour Tbs found since the start of the algorithm.

In tour construction, the methods in EAS is the same as the methods in AS. In pheromone updating, the pheromone evaporation formula is the same as (3.4).The additional reinforcement of tour Tbs is achieved by adding a quantity e/Cbs to its arcs, where e is a parameter that defines the weight given to the best-so-far tour Tbs , and Cbs is its length. The following is the equation for the pheromone deposit:

where Δτijk is defined as in Eq. (10) in AS and Δτijbs is defined as follows:

The key is that the EAS has adopted a daemon action, which is the additional incentive of the elite ants. Although this operation belongs to the pheromone updating steps, it is a kind of additional guidance to the overall operation. One can image like this: After all the ants including the elite ant depositing pheromone on the arcs they have crossed over their tour, the elite ant give the Tbs additional release of pheromones.

Compared with EAS, the pheromone updating mechanism in the AS is weak indeed. Sometimes, the optimal path may be with a very small difference between the paths which are not so satisfactory, and the mechanism in the AS can not make a good distinction between them. This is because of the simple form of the pheromone depositing formula which allows all ants use the same weight for depositing pheromone. Usually, the level for the algorithm to explore the overall optimal solution is not enough. The EAS with the parameter e determining the weight gives the best-so-far tour, that is, the best-so-far solution has been improved in the course of the search status, and the algorithm attempts to search a better solution which around the best-so-far solution. From the other side of the coin, the EAS concentrates a smaller search space which is compressed from original search space. Such a smaller search space may have better solutions. The mechanism increases the probability for finding overall optimal solution, and at the same time it also speeds up the convergence. Experiments later in this chapter show that parameter e needs to be selected with a reasonable value: an appropriate value for parameter e allows EAS to both find better tours and have a lower number of iterations. If parameter e is too small, the elitist strategy will not have much effect because of the low discrimination for better ants, while if the parameter e is too large, the algorithm will be too dependent on the initial best-so-far tour, and have rapid convergence to a small number of local optimal solutions, which weaken algorithm’s ability of exploration.

3.4. Rank-Based Ant System (AS_rank)

Another improvement over AS is the rank-based version of AS (Dorigo & Stützle, 2004; Bullnheimer et al., 1997): AS_rank. In AS_rank, before updating the pheromone trails, the ants are sorted by increasing tour length and the quantity of pheromone deposited by an ant is weighted according to the rank of the ant. Usually in each iteration, only the (w-1) best ranked ants and the ant produced the best-so-far tour (this ant is not necessarily is belong to the set of ants of the current algorithm iteration) are allowed to deposit pheromone. The best-so-far tour gives the strongest feedback, with weight w; the r-th best ant in the current iteration contributes the pheromone updating with the value 1/Cr multiplied by a weight given by max {0, w-r}.

In tour construction, the methods in AS_rank are the same as the methods in AS. In pheromone updating, the pheromone evaporation formula is the same as in Eq. (8). The AS_rank pheromone update rule is:

where Δτijr=1/Cr and Δτijbs=1/Cbs ; Cr is the length of r-th solution and Cbs is the same as in Eq. (12)

Compared with AS the AS_rank selects w ants to deposit pheromone according to the rank of solutions’ quality, which is a new improved formula. It completely abolishes the national pheromone deposit mechanism, in other word, only the ant who has constructed a good enough solution can deposit pheromone and the amount of pheromone to deposit is decided by the rank. It can reduce the operation of the pheromone and get rid of bad solutions directly. The pheromone depositing mechanism in AS_rank is a group of elite ants award, and it is better than the mechanism in AS which just depends on the reciprocal of the tours. Totally, the performance of AS_rank is much better than AS.

AS_rank and EAS are different on the reward strategy of elitist ants. EAS just give the best-so-far solution an additional incentive, although it can find the good solutions, indirectly it is greatly weakened or even abandoned the second-best-so-far solutions whose neighborhood may have better solutions. In AS_rank, the algorithm gives a group of elitist ants award, but the award is according to the rank of solutions. On the one hand the mechanism takes into account the importance of the best-so-far solution, which ensures that the experience of the leading elitist ant is retained; on the other hand it considers the neighborhood of sub-optimal solutions, that is, it increases ability to explore the optimal solution.

3.5. MAX-MIN Ant System (MMAS)

Max-Min Ant System is a constructive amendment to AS (Dorigo & Stützle, 2004, , Stützle & Hoos 1996, Stützle & Hoos, 1997, Stützle & Hoos, 2000). It is one of the best ACO algorithms. There are four main modifications with respect to AS:

1.It strongly exploits the best tours found: only either the iteration-best ant, that is, the ant produced the best tour in the current iteration, or the best-so-far ant is allowed to deposit pheromone.

It limits the possible range of pheromone trail values in the interval [τmin,τmax] .

The pheromone trails are initialized to the upper pheromone trail limit together with a small pheromone evaporation rate.

Pheromone trails are reinitialized when the system approaches are stagnated or when no improved tour has been generated for a certain number of consecutive iterations.

The first point strongly exploits the best tours found, but may lead to a stagnation situation in which all the ants follow the same tour. To increase the effect, the second modification is introduced by MMAS. It makes that the pheromone in each arc does not accumulate too large or consume too small, so that it could ensure the sustainability of exploration. The third point makes MMAS have a stronger ability to explore at the initial stage. The fourth point introduces a new mechanism which can be applied to all the ACO algorithms, that is, through the resumption of initial pheromone and re-search thus it increases the possibility of finding the optimal solution.

Since only the optimal solution is allowed to deposit pheromones, the formula about pheromones depositing is very concise:

where Δτijbest=1/Cbest . Cbest is the length of optimal tour. It can be the best-so-far tour or iteration-best tour.

In general, in MMAS implementations both the iteration-best and the best-so-far updating rules are used, in an alternate way. Obviously, the choice of the relative frequency with which the two pheromone updating rules are applied has an influence on how greedy the search is: When pheromone updating is always performed by the best-so-far ant, the search focuses very quickly around the best-so-far solution, whereas when it is the iteration-best ant that updates pheromones, the number of arcs received pheromone is larger and the search is less directed.

In MMAS, there are two default daemon actions. One is the pheromone trails limits, the other one is the pheromone trails re-initialization. In MMAS, lower and upper limits τmin and τmax on the possible pheromone values of any arc are imposed in order to limit the probability pij of selecting a city j when an ant is in city i in the interval [Pmin,Pmax] , so that it could avoid searching stagnation and enhance the ability to explore. In the long run, the upper pheromone trail limit on any arc is bounded by 1/ρC* , where C* is the length of the optimal tour. Based on this result, MMAS uses an estimate of this value of 1/ρCbs , to define τmax : each time a new best-so-far tour is found, the value of τmax is updated. The lower pheromone trail limit is set to τmin=τmax/a , where a is a parameter which decides the proportion of upper limit and lower limit. In MMAS, in order to avoid stagnation, the lower pheromone trail limits play a more important role than upper limits. Pheromone trail re-initialization is typically triggered when the algorithm approaches the stagnation behavior or if no improved tour of a given number of algorithm iterations is found. It can increase the exploration of paths that have only a small probability of being chosen. By the way, the pheromone trail re-initialization also increases the probability of finding the global optimal solution (equivalent to cumulative probability).

3.6. Ant Colony System (ACS)

Ant Colony System (ACS) is an extension of AS (Dorigo & Stützle, 2004; Dorigo & Gambardella, 1997). The ACS exploits the search experience accumulated by the ants more strongly than AS. The rule is called pseudorandom proportional rule. When located at city i, ant k moves to a city j, the rule is given by

where q is a random variable uniformly distributed in [0,1] , J is a random variable selected according to the probability distribution given by Eq.(5) with α = 1. q0(0≤q0≤1) is a parameter with which the ant makes the best possible move as indicated by the learned pheromone trails from the heuristic information, while with probability (1−q0) it performs a biased exploration of the arcs.

This pseudorandom proportional has strong artificial operability because the parameters q₀ can be set to guide the algorithm’s preference. By tuning the parameter q_0, it is allowed to modify the degree of exploration and to select whether to concentrate the search of the system around the best-so-far solution or to explore other tours.

In ACS only the best-so-far ant is allowed to update pheromone after each iteration including the pheromone deposit and pheromone evaporation. In each time an ant uses an arc to move from city i to city j, which is called the local pheromone updating to remove some pheromone from the arc to increase the exploration of alternative paths. The global pheromone trail updating is described as the following equation:

where Δτijbs=1/Cbs . Cbs is the length of the best-so-far tour Tbs . It is worth to note that the deposited pheromone is discounted by a factor ρ, which results in the new pheromone trail becoming a weighted average between the old pheromone value and the amount of pheromone deposited.

The local pheromone trail updating is described as the following equation:

where ε, 0≤ε≤1 , ε is the local pheromone evaporation rate and the τ0 is the initial pheromone. The effect of the local updating rule is to reduce pheromone trail of an ant in some arc so that the arc becomes less desirable for the following ants. The role of local updating is to strengthen the capacity of artificial ant exploration.

3.7. Application of five ACO algorithms in C-TSP and results analysis

In order to demonstrate and verify the properties of five ACO algorithms mentioned above, in this section, we’ll apply them in the Chinese TSP with 31 capital cities, and compare their advantages and disadvantages in the medium-scale discrete optimization problem.

The important parameters of Ant System are as follows. N: the number of cities; m: the number of artificial ants; α: the parameter that controls the relative importance of pheromone trail; β: the parameter that controls the relative importance of heuristic information; ρ: pheromone evaporation rate; q: a constant represents the weight of the deposited pheromone; τ0 : initial pheromone trail; NC: preset number of iterations. In addition, in its improved algorithms the addition parameters are as follows. The e in EAS : a parameter that defines the weight given to the best-so-far tour. The w in AS_rank: the number of the artificial ants who are allowed to add pheromone. The τ_proportion, τ_max, τ_min and nowbest_p in MMAS: τ_proportion is a parameter which decides the proportion of upper limit and lower limit; τ_max is the upper limit of pheromone; τ_min is the lower limit of pheromone; nowbest_p is the frequency of selecting best-so-far tours rule of depositing pheromone. The pbest and local_pin ACS: pbest is a probability of choosing right path; local_p is the local pheromone evaporation rate.

A candidate list is first built to restrict the number of available choices considered at each construction step. In general, candidate lists contain a number of the best rated choices according to some heuristic criterion. First, configure for each city a nearest neighbor list which records the other cities sorted in ascending order by distance; Second, build the candidate lists for each city and set the parameter nn(0≤nn≤n) which decides the number of the nearest neighbors needed; Last, get the cities which are previous nn cities in the nearest neighbor list into the candidate list for each city. When an ant constructs solution, it gives priority to the candidate list of cities. In fact, the ant usually considers from the first city in the candidate lists that has not been visited and selects with random probability rule. When all the cities in the candidate list have been visited by an ant, the ant will consider other cities and select the city which has a maximum value of [τij]α[ηij]β , that is, the ant selects the city which is the best experience-oriented one. Set the parameter nn to a constant which is below the number of cities n, especially for a small value, the algorithm’s speed will be improved significantly. The mechanism is feasible, because in TSP a good path can not appear a city i connects another city j which has a long distances from city i. In other words, the ant in city i should choose the city j which nears the city i. It is worth to note that the parameter nn is important for candidate list. If the value of nn is too large, the effect of speeding up algorithm will be weakened. On the other hand, a too small nn will make the performance of algorithm very poor. However, it should be noted that the use of truncated nearest-neighbor lists can make it impossible to find the global optimal solution. The global optimal solution does not mean to be the combination of cities in the candidate lists. Perhaps in order to achieve the best, ants in some cities should choose far way cities to go and these cities are not in departure cities’ candidate lists.

The steps of Ant System (AS) algorithm is as follows (in the case no candidate lists):

Step 1.Enter an actual TSP, get the scale n of the problem, and transform the instance into a symmetric distance matrix distance (set the diagonal elements with small values to prevent the situation of NAN.)

Initialize all the parameters, including m, α, β, ρ, q, τ0 and NC. Set the iteration number to 0.

Initialize storage variable including best-so-far solution nowbest_opt = 2 * vicinity (vicinity is the solution comes form the nearest algorithm), best-so-far path nowbest_path = zeros (1, n), pheromone tails matrix τ = ones(n) * τ0 , and the matrix which describes the importance of heuristic information τβ = dist.^(- β).

Begin to circulate, and set the iteration number nc = nc +1.

The starting cities of ants are randomly distributed as begin_city = randperm (n); initialize tabu list tabu = ones (m, n) and taboo the starting city; initialize path matrix path = zeros (m, n) and add the starting city into the first column; build the matrix that describes the importance of pheromone τα = τ.^ α; build the comprehensive weight matrix τd=τα.*τβ .

Ant walking steps step=step+1 (initialize step = 0，and step≤n−1 ).

Artificial ants’ label k = k +1 (initialize k = 0, and k≤m )

Choose the next city in accordance with the probability formula, taboo the selected city in the ant taboo list tabu (k, :), and add this chosen city into the ant path sequence path (k, step +1).

If k < m, then turn to step 7; otherwise turn to step 10.

According to the pheromone updating formula of AS, update the pheromone trails τ, and update the optimal solution as best-so-far solution nowbest_opt and the optimal path as best-so-far tour nowbest_path which includes the cities of the best-so-far tour.

If the iteration number nc < NC, then turn to step 4; otherwise end the operation, export the data and image results.

Remark: The steps described above are just the AS algorithm. The improved algorithms mentioned above should have some changes, which are mainly in the steps of the parameter settings, constructing solutions and pheromone updating. In the process of building a solution, they use parallel mechanism. To add the candidate lists to the algorithm, first of all one needs to add a step in one of the first two steps in the algorithm: list the neighbor cities in ascending order by distance and configure nearest neighbor list for each city, set the parameter Neighbor_num of the candidate list, and then get the candidate list Neighbor_list from the nearest neighbor list; followed by modifying the step 8: first of all, begin to consider the first city in the candidate list that does not be visited, and select the target city of the next step with probability rules. When the cities of the candidate list have all been visited one compares the values τd=τα.*τβ of all the other remaining cities, and select the largest one as the next target city; taboo the selected city in the ant taboo list tabu (k,:), and add this city into the ant path sequence path (k, step +1).

Table 2 is the summary of the experimental results of the five ACO algorithms including adding the mechanism of the candidate list and pheromone re-initialization, etc.. Each system runs for 100 times. In Table 2, the algorithm with "_C" has the mechanism of candidate list, the algorithm with “_R” has the mechanism of pheromone re-initialization, Best or Worst mean the best or worst solution of the 100 solutions. “Opt. rate” is the rate of global optimal solution, Average is the average solution of all solutions. “Average converg.” is the average convergence iteration number, and Average time is the average time in 100 running solutions. Relative error is described as follows:

where the global optimal solution = 15404, and the average solution is the average solution in 100 running solutions.

One can obtain from Table 2 the following results:

1. In the test of all 12 kinds of algorithms, from the column "Best" solutions one can see except AS and AS_C which add the mechanism of candidate list to AS, the other 10 kinds of algorithms can detect the global optimal solution 15404.

Along all the algorithms, from the column "Opt. rate" one can see max-min ant system which has the mechanisms of the candidate list and pheromone initialization added in MMAS_C_R owns the highest rate of excellent, that is 73%, followed by MMAS_R; Except AS and AS_C, ACS and ACS_C which add the mechanism of candidate list to ACS have the worse performance, that is 40%.

Compare these five algorithms without any mechanism, from the "Opt. rate " one can also see that rank based version AS, AS_rank has the highest rate of excellent, that is 63%, followed by MMAS, EAS, ACS, and AS.

Compare these five algorithms without any mechanism to those five algorithms with mechanisms of the candidate list respectively, that is, compare XXX with XXX_N, one can see from the optimal rate and average running time except ACS and ACS_C are not obvious, the mechanism of the candidate list slightly improves the performance of the algorithms and greatly reduces the running time.

Compare the four cases of MMAS, that is, MMAS, MMAS_C, MMAS_R and MMAS_C_R, one can see from the optimal rate that every mechanism can improve the performance of algorithms, and max-min ant system only with the mechanisms of pheromone initialization the candidate list has a better performance than only with the mechanisms of the candidate list. Of course, max-min ant system with two mechanisms at the same time has the best performance.