A new methodology of simulated annealing for the optimisation problems

Physica A 205 (1994) 367-374 North-Holland

~

~

S S D I 0378-4371(93)E0475-T

A new methodology of simulated annealing for the optimisation problems Simon C. Lin and James H.C. Hsueh Computing Centre, Academia Sinica, Taipei, Taiwan 11529, ROC

Complex optimisation problems with many degrees of freedom are often characterised by the enormously large configuration space, typically •(e N) or ~(N!). The idea of simulated annealing (SA) proposed by Kirkpatriek has been applied to the complex optimisation problems, which can be treated as annealing a statistical mechanical system from high temperature to low temperature; however, the SA is terribly slow for large problem sizes in typically (9'(N3 In N) time. We discover the hybrid algorithm (HA), which is based on a hybrid mechanism which combines conventional heuristics with low temperature simulated annealing (LTSA), which could be parallelised easily. The HA is a new approach of resolving optimisation problems with (7(N) complexity where information propagation can be inhibited by restraining the range of searches in the configuration space. We use the HA to resolve several famous combinatorial optimisation problems, including the travelling salesman problem (TSP) of large sizes up to 1 000 000 cities within 3 to 5 percent of the optimal value in linear time and other nonuniformly distributed TSPs as well. We shall also discuss the applicability of the HA to the optimisation problems in general.

1. Introduction Optimisation p r o b l e m s are ubiquitous. T h e y a p p e a r as P-type or N P - h a r d p r o b l e m s in combinatorics, as searching for the g r o u n d state in spin glasses, as learning m e c h a n i s m in neural networks, as partition, p l a c e m e n t , and wiring p r o b l e m s in the physical design of c o m p u t e r s , as action principles in physical p r o b l e m formulation, as linear p r o g r a m m i n g f o r m u l a t i o n , as p r o b l e m s in p a t t e r n recognition, even as natural language parser, etc. All the a b o v e p r o b l e m s are c o n c e r n e d with looking for the configuration with the m o s t optimal cost function, which is a function or functional defined o n the configuration space. A c o m m o n p r o p e r t y of the h a r d optimisation p r o b l e m s is the e n o r m o u s l y large size, typically ~7(eN) or ~ ( N I ) , of its configuration space. Since the optimisation p r o b l e m is essentially a searching p r o b l e m , for such an e n o r m o u s size of solution space, the exhaustive search b e c o m e s impossible and the heuristics are indispensable [1]. O n l y in recent years, the idea of simulated annealing (SA) p r o p o s e d by 0378-4371/94/$07.00 (~ 1994- Elsevier Science B.V. All rights reserved

368

S.C. Lin, J.H.C. Hsueh / Simulated annealing for the optimisation problems

Kirkpatrick [2] has been applied to complex optimisation problems, which can be treated as annealing a statistical mechanical system from high temperature to low temperature. SA is a stochastic searching strategy that prevents searches from becoming trapped in local minima such that the global optimal configuration of the optimisation problem could be found. In general, annealing schedules contain no specific knowledge about the nature of the particular optimisation problem. The cost function of the optimisation problem can be regarded as analogous to the internal energy of the statistical mechanical system, and the optimal configuration is analogous to the ground state (lowest energy state) configuration. The conventional fast SA (FSA) is at least an ~(N 3) algorithm when the quality of optimal configuration is not a major concern. If a good quality optimum is required, say, within 3 to 5 percent of optimum, then the typical complexity could be O(N 3 InN) and the upper bound could be around O(N 4 In N) [3]. The travelling salesman problem (TSP) is a generic problem belonging to the class of NP-complete problems in combinatorial optimisation, which are typically easy to state but hard to resolve [4]. In general, both the constructive and iterative heuristics for the approximate solution of the TSP are of the complexity G(N2) or 6(N3). These heuristics can achieve only around 30% above the optimum and only smaller size problems could ever be resolved [3]. Besides, parallelisation or vectorisation for the TSP are either too difficult or inefficient to implement [5]. Since the TSP is one of the most difficult problems in combinatorial optimisation, the concepts and techniques learnt from the treatment of the TSP may be applied to other popular optimisation problems and even to problems beyond the scope of traditional combinatorial optimisation. We shall use the TSP as a generic example to illustrate our hybrid algorithm [6] in the following section.

2. Scaling of cost function and control parameter

The probability distribution of energy at a certain temperature peaks around an equilibrium average value of energy in a Gaussian form. This Gaussian peak moves from high value of energy to lower energy as the system is annealed from high temperature. In other words, temperature provides a natural classification or partition of the configuration space that it is more probable to find high-energy configurations in the high temperature and low-energy configurations in lower temperatures. Therefore, temperature plays a similar role as the hash function in the search algorithms; temperature is also called the control parameter in the SA approach to the optimisation problem. Since

S.C. Lin, J.H.C. Hsueh / Simulated annealing for the optimisation problems

369

the optimisation problem is essentially a searching problem, the analogy of the control parameter to the hash function would imply that one should attempt to approach the zero control parameter directly. However, this is called quench in the language of statistical mechanics and the usual implication of quench is being trapped in the local minimum such that relaxation to a true global minimum could be extraordinarily long. This is the essence of the SA approach to the optimisation problem that one should anneal slowly, instead of quench, and stays long enough in each control parameter, so the equilibrium configuration at each temperature can be reached. One of the often neglected concerns of the complexity of algorithms is the consequence of the scaling of problem size on the quantities of interest, for example, the cost function of the TSP in a unit square scales as ~(VN) when the number of cities increases by N. In general, the scaling of cost function would imply the scaling of the control parameter as well; their relation can be qualitatively given by 1 dS -T = d E '

(1)

where T is temperature (control parameter), E is energy (cost function), S is entropy (information content), dS and dE are changes of the respective quantities. Note that entropy always scales as 6(N) since it is a logarithm of a quantity of the order of 6(eN). Taking the TSP in a unit square as an example, when the cost function scales as O(X/-N) then the control parameter would scale as 6(1/X/-N). In other words, T V ~ is a problem-size scale-invariant quantity and the corresponding value of T is smaller for larger problem sizes. This explains why the control parameter is a size dependent parameter in most treatments of the TSP by conventional FSA. Also note that the concept of scaling discussed here concerns the scaling of problem size, not the scaling of length scale as in the renormalisation group or critical phenomena [7]. We reckon that the scaling analysis, due to the increase of problem size, of cost function and related quantities in an optimisation problem is a crucial key to understand its algorithmic complexity issue.

3. Hybrid mechanism

We discover a general hybrid algorithm (HA) [6] that is based on a hybrid mechanism which combines conventional heuristics with low temperature simulated annealing (LTSA). In other words, the conventional SA can be

370

S.C. Lin, J.H.C. Hsueh / Simulated annealing for the optimisation problems

replaced by a two-phase H A , the front-end heuristic and the back-end L T S A with an adhesive hybrid mechanism. That hybrid mechanism is essentially a clever use of micro-canonical ensemble Monte Carlo (t~CEMC) simulation [8] which may serve as a thermometer to measure the temperature of configurations in SA. For the configurations generated by conventional heuristics, txCEMC simulation gives its temperature and the simulated annealing can be continued smoothly from that temperature. Since there exists very fast heuristics, but with relatively poor optima, the employment of our hybrid mechanism would save rather considerable computing time when compared with conventional SA for a near-optimal solution. For the TSP, many interesting front-end heuristics, for instance, spacefilling curve [9], nearestneighbour, and a general quench heuristic, have all been studied. A plot of the cost function against the control parameter for the TSP is shown in fig. 1, where the H A curve is overlaid upon the conventional SA curve. It seems that these two curves match almost perfectly from the plot. In the H A , we introduce a parameter 3, between 0 and 0.5 which confines the distances in Lin's opt-loop between two viable bonds [10]. The introduction of y gives rise to questions concerning the choice of 3, at different problem sizes. One finds from the simulation data that the decrease of y would cause the increase of acceptance ratio R, i.e.

175

P

/,J

15.0

i

Z

125

!0.0

<~> 7.5 Ii

5.0

/

/ 2.5

0.0 ~ 0.001

Y T

~

r 0.01

,'

/

/

~

, ,,,,iF ~ , , - i 1 10

0.1

, , , ,rTr, 100

T Fig. 1. Semilog plot of the average normalised optimal length versus temperature T for the N = 1000 uniform TSP. The data marked with crosses are the result of our hybrid algorithm.

S.C. Lin, J.H.C. Hsueh / Simulatedannealingfor the optimisationproblems

R(7, NM, T ~ ) = - ~

1

R(7, N, T V ~ ) ,

371

(2)

where R is taken at fixed problem size N. One also finds that decreases inversely with respect to N,

R(T/M, N, TX/'N) = MPR(7, N, TX/-N),

(3)

where R is fixed and the size changes from N to NM. Then, one can render R constant through 1

R('y/ M, NM, TX/-N-M) = - - ~ R( 7/ M, N, TX/N) = MP-~R(7, N, T V N ) , p -o'=0, In R(NM, T N ~ )

= -~r In M + In R(N, T V N ) ,

In R(7/M)

= p In M + In R(7 ) .

(4) (5) (6)

We find from the simulation data that indeed or = p - 1,

p = 0.988 +- 0.018,

or = 0.996 -+ 0.029.

Thus, the major problem of low acceptance ratio in the low temperature region occurring in the conventional SA has been cured by this scaling relation. Recall that our simulation is done in the very low temperature region, so most of the phenomenologically sophisticated annealing schedules are not optimised for this particular range of temperature. We argue that only simple annealing schedules are needed. The details of our hybrid algorithm appeared in other articles we have written [6]. We have pushed the H A up to 1 million cities; it is at least 100 times larger than anyone has ever tried [11]. For obvious reasons, it is not possible to show the 1 million cities configuration here. Instead, a typical 10 thousand cities configuration is shown in fig. 2. The total time spent on the N = 1 000 000 TSP is 50 hours on the HP750 workstation, the SFC heuristic plus (p~CEMC) simulation took 28 hours and the L T S A took only 22 hours to complete. The best possible normalised path length [12] is also extrapolated based on real simulation data to be 0.76. What we have achieved by our H A is to find a linear-time algorithm for the uniformly distributed NP-complete geometric TSP. Based on our result of extrapolation, it seems that our algorithm approaches asymptotically the optimal cost function in a controllable way. The hybrid algorithm is also capable of dealing with non-uniform distribution of city positions. We show the optimal configuration for exponen-

372

S.C. Lin, J.H.C. Hsueh l Simulated annealing for the optimisation problems N=I 0000 1.0o 2

0,0o 0. 00

I 00 =0. 748

Fig. 2. Plot of the N = 10000 uniform geometric TSP optimal configuration; it is produced by the

H A with SFC heuristic.

tial and normal distribution of city position in figs. 3 and 4, respectively. Applications of the hybrid mechanism to other heuristics and a parallelisable version of our scalable linear algorithm have also been done [6]. N=IO000

0.9

//

\,,

0 . 8

t

] r 1

\

\~'\

07L\ 0.6 0.5 0.4 0.3 0.2 0.1 0 0

0.2

0.4

0.6

0.8

1

L/sqrt(N) = 0.3234

Fig. 3. Plot of optimal configuration for N = 10000 exponentially distributed city positions.

S.C. Lin, J.H.C. Hsueh / Simulated annealing for the optimisation problems

373

N=I 0000 1.00

>

0 . 0 0

, ~

O. O0

r77

.....................

I~7

~ . . . . . . . . . . . ~. . . . . .

= O. 4 7 5 8

1. O0

Fig. 4. Plot of optimal configuration for N = 10000 normally distributed city positions.

4. Conclusion

The minimum complexity of typical scientific computing problems often requires at least ~7(N) when the problem size increases by N and at least another ¢(N) for the information propagation (IP). IP happens in an iterative approximation scheme and must be passed from one end to the other end of the system in order to ensure convergence or equilibrium of the computational processes. The secret to deal with the complexity of computation problems lies in the treatment of the IP. What we have presented in this article is a new approach of resolving optimisation problems with ~7(N) complexity where the IP can be inhibited by restraining the range of searching in the configuration space. The extraordinary speed-up of the H A has been proved to be applicable to geometric optimisation problems; however, the applicability of the H A to P2-type problems [13] remains an open question.

References [1] S.E Bradley, A.C. Hax and T.L. Magnanti, Applied Mathematical Programming ( A d d i s o n Wesley, MA, 1977); R.G. Parker and R.L. Rardin, Discrete Optimization (Academic Press, San Diego, 1988). [2] S. Kirkpatrick, C.D. Gelatt and M.P. Vecchi, Science 220 (1983) 671;

374

S.C. Lin, J.H.C. Hsueh / Simulated annealing for the optimisation problems

S. Kirkpatrick, J. Stat. Phys. 34 (1984) 975. [3] P.J.M. van Laarhoven and E.H.L. Aarts, Simulated Annealing: Theory and Applications (Kluwer, Dordrecht, 1987); E. Aarts and J. Korst, Simulated Annealing and Boltzmann Machines (Wiley, Chichester, 1989). [4] E.L. Lawler, ed., The Travelling Salesman Problem (Wiley, Chichester, 1985). [5] G. Fox et al., Solving Problems on Concurrent Processors, vol. 1 (Prentice Hall, N J, 1988). [6] S.C. Lin, H.C. Hsueh and R.S. Lin, An accelerated hybrid algorithm of the travelling salesman problem, preprint; S.C. Lin and H.C. Hsueh, Nearest-neighbour heuristic in accelerated algorithms of optimisation problems, preprint; A study on quench heuristic in the hybrid algorithms of the travelling salesman problem, preprint. S.C. Lin and H.H. Fu, Parallelisable algorithms for the travelling salesman problem, preprint. [7] S.K. Ma, Modern Theory of Critical Phenomena (Benjamin, Reading, MA, 1976). [8] M. Creutz, Phys. Rev. Lett. 50 (1983) 1411; S.C. Lin, C.H. Kao and S.T. Li, Monte Carlo simulation in microcanonical ensemble, presented at STATPHYS 17, Rio de Janeiro, unpublished. [9] L.K. Platzman and J.L. Bartholdi, J. ACM 36 (1989) 719; H.-O. Peitgen, H. Jurgens and D. Saupe, Fractals for the Classroom, part one (Springer, New York, 1992). [10] S. Lin, Bell Syst. Tech. J. 44 (1965) 2245. S. Lin and B. Kernighan, Oper. Res. 21 (1973) 498. [11] E. Bonomi and J-L Lutton, SIAM Rev. 26 (1984) 551. [12] J. Beardwood, J.H. Halton and J.M. Hammersley, Proc. Cambridge Philos. Soc. 55 (1959) 299; D. Stein, Scheduling Dial-a-Ride Transportation System: An Asymptotic Approach, Ph.D. thesis, Harvard Univ. (1977). [13] M. Mezard, G. Parisi and M.A. Virasoro, Spin Glass Theory and Beyond (World Scientific, Singapore, 1987); J.L. van Hemmen and I. Morgenstern, eds., Heidelberg Colloquium on Glassy Dynamics, Lecture Notes in Physics, vol. 275 (Springer, Berlin, 1987).