t AMERICAN JOURNAL OF MATHEMATICAL AND MANAGEMENT SCIENCES CopyllgMQ ISM by American Sclonur Pmr, Inc. NEIGHBORHOOD SIZE IN THE SIMULATED ANNEALING ALGORITHM Larry Goldstein Michael Waterman University of Southern California Department of Mathematics Loa Angeles, CA 90089-1113 SYNOPTIC ABSTRACT Simulated annealing is a probabilistic algorithm that has shown some promise when applied to combinatorially NP-hard problems. One advantage of the simulated annealing algorithm is that it is based on an analogy with statistical mechanics which is not problem-specific. How- ever, any implementation of the algorithm for a given problem requires that several specific choices be made. The success or failure of the proce- dure may depend on these choices. In this study we explore the effect of choice of neighborhood size on the algorithm’s performance when applied to the travelling salesman problem. Key Words and Phrases: Simulated Annealing, travelling salesman, neighborhoods 1988, VOL. 8, NOS. 3 & 4, 409-423 0196-6324/88/030409-15 $20.00 409 " 4 10 GOLDSTEIN & WATERMAN 1. INTRODUCTION. Certain discrete versions of the simulated annealing algorithm are probabilistic approachs to combinatorially NP-hard problems, that is, to a class of problems for which no polynomial time algorithms are known. The key to the simulated annealing algorithm is an analogy with sta- tistical mechanics which is not problem specific. In a typical discrete optimization problem, one is given a finite set S, (typically large), and a cost funtion f, and seeks u E S such that f(u) is a minimum. One may regard the function f as the energy function of some physical system; if one could now simulate the cooling of this system, a state of minimum energy would be obtained. Although this parallel is universal, any imple- mentation of the algorithm requires choices to be made that are specific to the problem at hand. In order to simulate a physical system, the algorithm proceeds e quentially by moving from one state to another by a certain probabilistic mechanism. h m an y given state s, there are a set of states, say N,, where transitions from s are allowed. We call N. the set of neighbors of 5 It is with the choice of neighborhoods N, that this study is con- cerned. It seems to be often overlooked that the performance of the simulated annealing algorithm depends critically on the choice of neigh- borhood structure, and more importantly, that one is free to choose a system that allows the algorithm to perform well. If the choice of neigh- borhood is too small, then the resulting simulated process will not be able to move around the set S quickly enough to reach the minimum in a SIMULATED ANNEALING 411 reasonable time. On the other hand, if the neighborhoods are too large, then the process essentially performs a random search through S, with next possible state chosen practically uniformly over S. The question now arises: what choice of neighborhoods N, will allow the algorithm to converge quickly? Intuitively, it seems that a neighborhood system that strikes a compromise between these extremes would be best. Neighborhood structure is not the only aspect of the simulated an- nealing algorithm that is free to be chosen in a way that improves the performance of the algorithm; the form of the energy function may also affect the behavior of the algorithm. For example, if f is nonnegative one may contrast an implementation of the simulated annealing algo- rithm that minimizes f with one that minimizes this problem is not studied here. 2. THE SIMULATED ANNEALING ALGORITHM. We now describe the simulated annealing algorithm in a general setting. We begin with the underlying neighborhood system. Consider a finite set S. For each s E S, suppose there is given a subset N. C S that satisfies 1. Vs E S,s E N,. 2. Vs,t E S,s E Nt if and only if t E N.. 4. Vs,t E S, there exists an integer m and ul,up,. .. ,urn in S such 412 GOLDSTEIN & WATERMAN for i = 1,2,. . . ,m - 1. (That is, we require the graph on S con- structed by joining two elements of S with an & whenever they lie in the same set N, for some s E S, to be connected.) We call such an indexed system of subsets {Nd},Esa neighborhood system. Now assume given a cost (or energy) function f, where f : S + R, it is required to locate the element of S that minimizes f. For each neighborhood system, we can consider an associated “greedy” algorithm as follows: Algorithm G ({N,} ,s E S) : Begin at any point SI E S. At stage n, choose sn+l to satisfy f(sn+l) = min{f(t) : t E Nd,.}. It is clear that after a finite number of iterations of algorithm G(N,), the state will become trapped in a local minimum of f. The simulated annealing algorithm is a probabilistic modification of the greedy algorithm G(N,) that does not get trapped in a local minimum. This is accomplished by occasionally accepting a new state that increasea the energy function. The idea of a simulation of this type was first introduced by Metropolis, Roaenbluth, Rosenbluth, Teller, and Teller (1953). For any given T > 0 the Gibbs distribution over S, assigns to s E S probability ezd-f (s)/T) *T(S) = ZT SIMULATED ANNEALING 413 where ZT (the “partition function”) is chosen so that the above prob- abilities sum to 1, that is ZT = Cexp(-f(s)/T). 8ES Note that for T > 0 small, the Gibbs distribution concentrates its mass on favorable states, that is, states s with small values of f(s),a nd this effect is more pronounced the smaller the value of T. One may easily construct a Markov chain that has the above distribution as its stationary law. As the Markov chain converges in distribution to this law, one may run the simulation for a time and find a state of low energy with high probability. The greedy algorithm G (N,)is essentially this Markov chain run for the case of T fixed at 0, whereas in the limit of high T all states are essentially weighted with the same probability and one is moving from a state to its neighbors uniformly. Of course, by the above mentioned analogy with statistical mechan- ics, T here is seen to play the role of temperature, and one may now suspect that T may be lowered as the simulation proceeds in order to force the system to a state of minimum energy. This idea is due to Kirk- Patrick, Gelatt, and Vecchi (1983). As with a physical system, tempera- ture may be lowered too rapidly and the system may become trapped in a local energy minimum, that is, the algorithm will too closely resemble the greedy algorithm G(N,). A theorem of Geman and Geman (1984) shows that if Tn = c/ log n, for c sufficiently large, then the system will in fact not be trapped. With this choice of Tn, the algorithm proceeds as follows. Algorithm SA ({N,},s E S): Choose an initial point s1 E S, uni- 414 GOLDSTEIN & WATERMAN formly over S. At time n, assume un given. h m th e set N,,, choose a point uniformly, say t. Calculate A = f(t) - f (sn) NOW,s et sn+l= t with probability p = erp( -A+/Tn), and set Sn+l = Sn - with the complementary probability 1 p, where A+={ A ifA>O 0 otherwise. In order to implement the simulated annealing algorithm SA (N,) one is required to furnish a neighborhood system N, and a cost function f. It is exactly these elements of the algorithm that are problem specific. We now turn to a specific problem, and a description of a neighborhood system for that problem. 3. LIN'S k-NEIGBORHOOD SYSTEM FOR THE TRAVELLING SALESMAN PROBLEM The travelling salesman problem models the salesman who is re- quired to visit a number of cities and return home covering minimal . . , distance. Let the "cities" c1, c2,. CN be independent and uniformly distributed in the unit square [0, l]', and let di,j denote the distance be- tween city i and city j. The finite set S over which we seek a minimum .. is the set of all permutations of { 1,2,. ,N );a given permutation gives the order in which the tour of the cities is to be taken. The cost (energy), function f that is to be minimized is the total length of the tour taken in the order dictated by the permutation u E S and can be written N-1 + f(s) = d,(i),,(i+l) d.(N),,(l)- i- 1 SIMULATED ANNEALING 415 This problem belongs to the class of NP hard problems, hence no poly- nomial time algorithm for its solution is known (see for example Garey and Johnson (1979)). In this particular problem, while there is a “natural” choice of an energy function, as stated in the introduction there is really no reason to believe that f will be prefered to some other function on S that attains fl its global minimum at the same optimal tour, such as or f2 for example. Given the function f as above, one is now only required to choose a neighborhood structure for the set of permutations S. In a study of deterministic algorithms for the travelling desman problem, Lin (1965) introduced the notion of k-optimality, which gives rise to a neighborhood structure for each k. In the terminology used here, a tour is k-optimal if it has the smallest cost of all tours in its neighborhood. The larger the value of k, the more neighbors any given tour will have. For h = 1 a tour is a neighbor of itself only and hence every tour is l-opt; for k = N every tour is a neighbor of every other tour, and hence only optimal tours are N-opt. For fixed k we define a system of neighborhoods M follows. Imagine that there is a link between any two cities in a tour. We say that two tours are neighbors if one can break k or less links in the one tour and reassemble to obtain the other. nom this definition it becomes clear that two toum are neighbors for k = 2 if and only if one tour can be obtained from the other by reversing the order of the cities in a portion of one of the tours. For k small relative to N, the number of k neighbors of any given tour is approximately (f)w2Thke fac.to r (f) counts the 416 GOLDSTEIN & WATERMAN number of ways k links may be broken from the N possible, the number of ways the k components of the broken tour may be reassembled, and the factor 2' counts the number of orientations possible for the k components. The formula is not exact since it ignores the possibility of having components of size 1, and so a factor of 2 should not be entered for this component; indeed, for k = N all components are of size 1 and C$ no factors of 2 enter yielding the correct answer of for the total number of possible tours. We have the term N - 1 as we are considering the tour to be in a loop, and we may consider it to begin at city 1; the factor of 2 takes care of the fact that a given tour and the same tour taken in reverse order are to be considered equivalent. For the cases of interest below, k is small relative to N and the formula above gives a reasonable approximation to the order of growth of the neighborhood size in k. With the above ingredients, that is, with a cost function and a neigh- borhood structure now fixed by a choice of &, we can implement the sim- ulated annealing algorithm SA (N#)O.u r interest below is to determine which value of k allows the simulated annealing algorithm of locate the minimum quickly. 4. EFFECT OF NEIGHBORHOOD SIZE ON SPEED OF CONVERGENCE. Bonomi and Lutton (1984) implemented a version of the simulated annealing algorithm with Lin's 2-opt neighborhoods and reported posi- tive results. In this study, we are interested in how di#erent choices of neighborhood system, that is different values of &, affect the performance SIMULATED ANNEALING 417 of the algorithm. In Bonomi and Lutton (1984), N points are laid down uniformly in the unit square [0,1l2 as described. This area is then subdivided into many smaller subsquares; using a path that tends to a space filling curve in the limit, a short path is found that tours the subsquares and the algorithm is then run independently among a group of subsquares. We will call this procedure “modified simulated annealing” for the travelling salesman problem. This modified procedure will speed up convergence to the minimum. In our study, we consider the unmodified version of the simulated annealing algorithm SA (N,).T hat is, we study the simulated annealing algorithm’s performance as a function of k without the above modifica- tion that speeds convergence. We have three reasons for making such a study. -First, t he heuristic used in Bonomi and Lutton (1984) is highly problem specific as it relies on the fact that points close together in [0,1]* are likely to be close together in the optimal tour. In fact, one takes advantage of knowing the average intercity distance in the optimal tour (see Bearwood, Halton, and Hammersley (1959) and the discussion below) and therefore, a priori, need only consider moves that result in intercity distances on this order. In many problems, among them even problems such as those in Goldstein and Waterman (1987) that bear significant resemblance to the travelling salesman problem, one does not have such a priori information about the solution, and therefore cannot build a heuristic that uses this information to advantage. Therefore we i 418 GOLDSTEIN & WATERMAN . retain more generality by considering the unmodified algorithm Second, it is preferable not to complicate the outcome of the sim- ulation with the choice of some particular heuristic that may affect the results of the study in an unknown way. That is, with the modification, the choice of k is confounded with the choice of heuristic. In short, our second reason for making this study is that throughout, we are more interested in the simulated annealing algorithm in general than its per- . formance for this problem in particular Lastly, even as applied to the problem at hand, if one were to adopt the subdivision approach for the travelling salesman problem, one would always be solving the unmodified version of the problem on subsquares anyway and would still like to be using the best value of k on each sub- problem. (In any implementation designed to actually solve the travel- ling salesman problem it would certainly be advisable to adopt a heuristic such as the subdivision approach in order to speed convergence). In most minimization problems, one is not usually given in advance the value of the cost function at the minimum. In fortuitous cases where this value is known, the information can be used to devise a stopping rule for a procedure to halt when it gets sufficiently close to the minimum. In addition, this value can be used to gauge how well an algorithm performs against such a standard. The travelling salesman problem with a uniform city distribution is an example of a problem where the value of the cost function is known at the minimum (that is, the optimal tour length is known), in a probabilis- tic senac in the limit for many cities. (For an example where the value of