f10-9

436 Chapter10. MinimizationorMaximizationofFunctionsSample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).to the usual method by any factor substantially larger than the “tender-loving-care factor” (which reflects the programming effort of the proponents). Problemswheretheobjectivefunctionand/oroneormoreoftheconstraintsare replacedbyexpressionsnonlinearinthevariablesarecalled nonlinearprogramming problems. Theliteratureonsuchproblemsisvast,butoutsideourscope. Thespecial case of quadraticexpressions is called quadraticprogramming . Optimization prob- lemswherethevariablestakeononlyintegervaluesarecalled integerprogramming problems, a special case of discrete optimization generally. The next section looks at a particular kind of discrete optimization problem. CITED REFERENCES AND FURTHER READING: Bland, R.G. 1981, Scientific American , vol. 244 (June), pp. 126–144. [1] Dantzig, G.B. 1963, Linear Programming and Extensions (Princeton, NJ: Princeton University Press). [2] Kolata, G. 1982, Science, vol. 217, p. 39. [3] Gill,P.E., Murray, W.,andWright, M.H. 1991, NumericalLinearAlgebraandOptimization , vol. 1 (Redwood City, CA: Addison-Wesley), Chapters 7–8. Cooper,L.,andSteinberg,D.1970, IntroductiontoMethodsofOptimization (Philadelphia:Saun- ders). Gass, S.T. 1969, Linear Programming , 3rd ed. (New York: McGraw-Hill). Murty, K.G. 1976, Linear and Combinatorial Programming (New York: Wiley). Land,A.H., andPowell,S.1973, FortranCodesforMathematicalProgramming (London:Wiley- Interscience). Kuenzi, H.P., Tzschach, H.G., and Zehnder, C.A. 1971, Numerical Methods of Mathematical Optimization (New York: Academic Press). [4] Stoer,J.,andBulirsch,R.1980, IntroductiontoNumericalAnalysis (NewYork:Springer-Verlag), §4.10. Wilkinson, J.H., and Reinsch, C. 1971, Linear Algebra , vol. II of Handbook for Automatic Com- putation(New York: Springer-Verlag). [5] 10.9 Simulated Annealing Methods Themethodofsimulatedannealing [1,2]isa techniquethathasattractedsignif- icant attention as suitable for optimization problems of large scale, especially ones wherea desiredglobalextremumis hiddenamongmany,poorer,localextrema. For practicalpurposes,simulatedannealinghaseffectively“solved”thefamous traveling salesman problem of finding the shortest cyclical itinerary for a traveling salesman who must visit each of Ncities in turn. (Other practical methods have also been found.) Themethodhasalsobeenusedsuccessfullyfordesigningcomplexintegrated circuits: The arrangement of several hundred thousand circuit elements on a tiny siliconsubstrateisoptimizedsoas tominimizeinterferenceamongtheirconnectingwires [3,4]. Surprisingly,the implementationof thealgorithmis relativelysimple. Notice that the two applications cited are both examples of combinatorial minimization . Thereisanobjectivefunctiontobeminimized,asusual;butthespace over which that function is defined is not simply the N-dimensional space of N 10.9SimulatedAnnealingMethods 437Sample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).continuouslyvariableparameters. Rather,itisadiscrete,butverylarge,configuration space, like the set of possible orders of cities, or the set of possible allocations ofsilicon “real estate” blocks to circuit elements. The number of elements in the configurationspaceis factoriallylarge,sothattheycannotbeexploredexhaustively. Furthermore, since the set is discrete, we are deprived of any notion of “continuingdownhill in a favorable direction.” The concept of “direction” may not have any meaning in the configuration space. Below,wewillalsodiscusshowtousesimulatedannealingmethodsforspaces with continuous control parameters, like those of §§10.4–10.7. This application is actuallymorecomplicatedthanthecombinatorialone,sincethefamiliarproblemof“long,narrowvalleys”againasserts itself. Simulatedannealing,aswewillsee,tries “random” steps; but in a long, narrow valley, almost all random steps are uphill! Some additional finesse is therefore required. Attheheartofthemethodofsimulatedannealingisananalogywiththermody- namics, specifically with the way that liquids freeze and crystallize, or metals cool andanneal. Athightemperatures,themoleculesofaliquidmovefreelywithrespecttooneanother. Iftheliquidis cooledslowly,thermalmobilityis lost. Theatomsare often able to line themselves up and form a pure crystal that is completely ordered overadistanceuptobillionsoftimesthesizeofanindividualatominalldirections. Thiscrystalisthestateofminimumenergyforthissystem. Theamazingfactisthat, forslowlycooledsystems,natureisabletofindthisminimumenergystate. Infact,ifa liquidmetal is cooledquicklyor “quenched,”it doesnot reachthis state but rather ends up in a polycrystallineor amorphousstate havingsomewhathigherenergy. So the essence of the process is slowcooling, allowing ample time for redistribution of the atoms as they lose mobility. This is the technical definition of annealing ,and it is essential forensuringthata lowenergystate will be achieved. Although the analogy is not perfect, there is a sense in which all of the minimization algorithms thus far in this chapter correspond to rapid cooling or quenching. Inall cases, we havegonegreedilyforthe quick,nearbysolution: Fromthe starting point, go immediately downhill as far as you can go. This, as often remarked above, leads to a local, but not necessarily a global, minimum. Nature’s own minimization algorithm is based on quite a different procedure. The so-calledBoltzmann probability distribution, Prob (E)∼exp(−E/kT )( 10.9.1 ) expresses the idea that a system in thermal equilibrium at temperature Thas its energy probabilistically distributed among all different energy states E.E v e n a t low temperature, there is a chance, albeit very small, of a system being in a high energystate. Therefore,there is a correspondingchanceforthe system to get outof a local energyminimumin favor of findinga better, more global,one. The quantity k(Boltzmann’s constant) is a constant of nature that relates temperature to energy. In otherwords, the system sometimes goes uphillas well as downhill;but the lower the temperature, the less likely is any significant uphill excursion. In 1953, Metropolis and coworkers [5]first incorporated these kinds of prin- ciples into numerical calculations. Offered a succession of options, a simulated thermodynamicsystem was assumed to change its configurationfrom energy E1to energy E2withprobability p=e x p [ −(E2−E1)/kT]. Noticethatif E2<E 1,this probability is greater than unity; in such cases the change is arbitrarily assigned a 438 Chapter10. MinimizationorMaximizationofFunctionsSample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).probability p=1,i.e.,thesystem alwaystooksuchanoption. Thisgeneralscheme, of always taking a downhill step while sometimes taking an uphill step, has come to be known as the Metropolis algorithm. TomakeuseoftheMetropolisalgorithmforotherthanthermodynamicsystems, one must provide the following elements: 1. A description of possible system configurations. 2. A generator of random changes in the configuration; these changes are the “options” presented to the system. 3. An objective function E(analog of energy) whose minimization is the goal of the procedure. 4. A control parameter T(analog of temperature) and an annealing schedule which tells how it is lowered from high to low values, e.g., after how many random changes in configuration is each downward step in Ttaken, and how large is that step. The meaning of “high” and “low” in this context, and the assignment of a schedule, may require physical insight and/or trial-and-errorexperiments. CombinatorialMinimization: The TravelingSalesman A concrete illustration is provided by the traveling salesman problem. The proverbialsellervisits Ncitieswithgivenpositions (xi,yi),returningfinallytohisor hercityoforigin. Eachcityis tobevisitedonlyonce,andtherouteis tobemadeas shortas possible. This problembelongsto aclass knownas NP-complete problems, whosecomputationtimeforan exactsolutionincreaseswith Nasexp(const.×N), becomingrapidlyprohibitiveincostas Nincreases. Thetravelingsalesmanproblem alsobelongsto aclass ofminimizationproblemsforwhichtheobjectivefunction E has many local minima. In practical cases, it is often enough to be able to choose fromtheseaminimumwhich,evenifnotabsolute,cannotbesignificantlyimprovedupon. Theannealingmethodmanagestoachievethis,whilelimitingits calculations to scale as a small power of N. Asaprobleminsimulatedannealing,thetravelingsalesmanproblemishandled as follows: 1.Configuration. Thecitiesarenumbered i=1...Nandeachhascoordinates (x i,yi). A configuration is a permutation of the number 1...N, interpreted as the order in which the cities are visited. 2.Rearrangements. An efficient set of moves has been suggested by Lin [6]. The movesconsist of two types: (a) A section of path is removedand thenreplaced withthesamecitiesrunningintheoppositeorder;or(b)asectionofpathisremoved andthenreplacedinbetweentwocitiesonanother,randomlychosen,partofthepath. 3.Objective Function. In the simplest form of the problem, Eis taken just as the total length of journey, E=L≡N/summationdisplay i=1/radicalbig (xi−xi+1)2+(yi−yi+1)2 (10.9.2 ) with the convention that point N+1is identified with point 1. To illustrate the flexibility of the method, however, we can add the following additional wrinkle: Supposethatthe salesmanhasan irrationalfearofflyingovertheMississippi River. 10.9SimulatedAnnealingMethods 439Sample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).Inthatcase,wewouldassigneachcityaparameter µi,equalto +1ifitiseastofthe Mississippi, −1if it is west, and take the objective function to be E=N/summationdisplay i=1/bracketleftBig/radicalbig (xi−xi+1)2+(yi−yi+1)2+λ(µi−µi+1)2/bracketrightBig (10.9.3 ) A penalty 4λis thereby assigned to any river crossing. The algorithm now finds the shortest path that avoids crossings. The relative importance that it assigns to lengthofpathversusrivercrossingsis determinedbyourchoiceof λ. Figure10.9.1 shows the results obtained. Clearly, this technique can be generalized to include many conflicting goals in the minimization. 4.Annealingschedule. This requires experimentation. We first generate some random rearrangements, and use them to determine the range of values of ∆Ethat willbeencounteredfrommovetomove. ChoosingastartingvaluefortheparameterTwhich is considerably larger than the largest ∆Enormally encountered, we proceed downward in multiplicative steps each amounting to a 10 percent decrease inT. We holdeach new value of Tconstant for, say, 100Nreconfigurations,or for 10Nsuccessful reconfigurations, whichever comes first. When efforts to reduce E further become sufficiently discouraging, we stop. The following traveling salesman program, using the Metropolis algorithm, illustrates the main aspects of the simulated annealing technique for combinatorial problems. SUBROUTINE anneal(x,y,iorder,ncity) INTEGER ncity,iorder(ncity)REAL x(ncity),y(ncity) C USES irbit1,metrop,ran3,revcst,revers,trncst,trnspt This algorithm finds the shortest round-trip path to ncitycities whose coordinates are in the arrays x(1:ncity),y(1:ncity) . The array iorder(1:ncity) specifies the order in which the cities are visited. On input, the elements of iorderm a yb es e tt oa n yp e r - mutation of the numbers 1toncity. This routine will return the best alternative path it can find. INTEGER i,i1,i2,idec,idum,iseed,j,k,nlimit,nn,nover,nsucc,n(6), * irbit1 REAL de,path,t,tfactr,ran3,alen,x1,x2,y1,y2 LOGICAL ansalen(x1,x2,y1,y2)=sqrt((x2-x1)**2+(y2-y1)**2)nover=100*ncity Maximum number of paths tried at any temperature. nlimit=10*ncity Maximumnumberofsuccessfulpathchangesbeforecontinuing. tfactr=0.9 Annealing schedule: tis reduced by this factor on each step. path=0.0 t=0.5 do 11i=1,ncity-1 Calculate initial path length. i1=iorder(i)i2=iorder(i+1) path=path+alen(x(i1),x(i2),y(i1),y(i2)) enddo 11 i1=iorder(ncity) Close the loop by tying path ends together. i2=iorder(1) path=path+alen(x(i1),x(i2),y(i1),y(i2))idum=-1iseed=111 do 13j=1,100 Try up to 100 temperature steps. nsucc=0do 12k=1,nover 1 n(1)=1+int(ncity*ran3(idum)) Choose beginning of segment .. 440 Chapter10. MinimizationorMaximizationofFunctionsSample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).0. 510.51 0. 510.51 0. 510.51(a) (b) (c) Figure 10.9.1. Traveling salesman problem solved by simulated annealing. The (nearly) shortest path among 100randomly positioned cities is shownin (a). Thedotted line isariver, butthere isnopenalty in crossing. In (b) the river-crossing penalty is made large, and the solution restricts itself to the minimum number of crossings, two. In(c) the penalty has been made negative: the salesman is actually a smugglerwho crosses the river on the flimsiest excuse! 10.9SimulatedAnnealingMethods 441Sample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).n(2)=1+int((ncity-1)*ran3(idum)) ..and end of segment. if (n(2).ge.n(1)) n(2)=n(2)+1 nn=1+mod((n(1)-n(2)+ncity-1),ncity) nn isthenumber ofcitiesnotonthe segment. if (nn.lt.3) goto 1 idec=irbit1(iseed) Decide whether to do a segment reversal or transport. if (idec.eq.0) then Do a transport. n(3)=n(2)+int(abs(nn-2)*ran3(idum))+1n(3)=1+mod(n(3)-1,ncity) Transporttoalocationnotonthepath. call trncst(x,y,iorder,ncity,n,de) Calculate cost. call metrop(de,t,ans) Consult the oracle. if (ans) then nsucc=nsucc+1path=path+de call trnspt(iorder,ncity,n) Carry out the transport. endif else Do a path reversal. call revcst(x,y,iorder,ncity,n,de) Calculate cost. call metrop(de,t,ans) Consult the oracle. if (ans) then nsucc=nsucc+1 path=path+de call revers(iorder,ncity,n) Carry out the reversal. endif endif if (nsucc.ge.nlimit) goto 2 Finish early if we have enough successful changes. enddo 12 2 write(*,*) write(*,*) ’T =’,t,’ Path Length =’,path write(*,*) ’Successful Moves: ’,nsucct=t*tfactr Annealing schedule. if (nsucc.eq.0) return If no success, we are done. enddo 13 returnEND SUBROUTINE revcst(x,y,iorder,ncity,n,de) INTEGER ncity,iorder(ncity),n(6)REAL de,x(ncity),y(ncity) This subroutine returns the value of the cost function for a proposed path reversal. ncity isthenumberofcities,andarrays x(1:ncity),y(1:ncity) givethecoordinatesofthese cities. iorder(1:ncity) holdsthepresentitinerary. Thefirsttwovalues n(1)andn(2) ofarray ngivethestartingandendingcitiesalongthepathsegmentwhichistobereversed. On output, deis the cost of making the reversal. The actual reversal is not performed by this routine. INTEGER ii,jREAL alen,xx(4),yy(4),x1,x2,y1,y2 alen(x1,x2,y1,y2)=sqrt((x2-x1)**2+(y2-y1)**2) n(3)=1+mod((n(1)+ncity-2),ncity) Find the city before n(1).. n(4)=1+mod(n(2),ncity) .. and the city after n(2). do 11j=1,4 ii=iorder(n(j)) Findcoordinatesforthefourcitiesinvolved. xx(j)=x(ii)yy(j)=y(ii) enddo 11 de=-alen(xx(1),xx(3),yy(1),yy(3)) Calculate cost of disconnecting the segment atbothendsandreconnecting intheop-posite order.* -alen(xx(2),xx(4),yy(2),yy(4)) * +alen(xx(1),xx(4),yy(1),yy(4)) * +alen(xx(2),xx(3),yy(2),yy(3)) returnEND 442 Chapter10. MinimizationorMaximizationofFunctionsSample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).SUBROUTINE revers(iorder,ncity,n) INTEGER ncity,iorder(ncity),n(6) Thisroutine performsapathsegment reversal. iorder(1:ncity) isaninputarraygiving the present itinerary. The vector nhas as its first four elements the first and last cities n(1),n(2) of the path segment to be reversed, and the two cities n(3)andn(4)that immediately precede and follow this segment. n(3)andn(4)are found by subroutine revcst. On output, iorder(1:ncity) contains the segment from n(1)ton(2)in reversed order. INTEGER itmp,j,k,l,nn nn=(1+mod(n(2)-n(1)+ncity,ncity))/2 This many cities must be swapped to effect the reversal. do11j=1,nn k=1+mod((n(1)+j-2),ncity) Start at the ends of the segment and swap pairs of cities, moving toward the cen- ter.l=1+mod((n(2)-j+ncity),ncity) itmp=iorder(k) iorder(k)=iorder(l)iorder(l)=itmp enddo 11 return END SUBROUTINE trncst(x,y,iorder, ncity,n,de) INTEGER ncity,iorder(ncity),n(6) REAL de,x(ncity),y(ncity) Thissubroutinereturnsthevalueofthecostfunctionforaproposedpathsegmenttransport. ncityis the number of cities, and arrays x(1:ncity) andy(1:ncity) give the city coordinates. iorderis an array giving the present itinerary. The first three elements of array ngive the starting and ending cities of the path to be transported, and the point among the remaining cities after which it is to be inserted. On output, deis the cost of the change. The actual transport is not performed by this routine. INTEGER ii,j REAL xx(6),yy(6),alen,x1,x2,y1,y2alen(x1,x2,y1,y2)=sqrt((x2-x1)**2+(y2-y1)**2) n(4)=1+mod(n(3),ncity) Find the city following n(3).. n(5)=1+mod((n(1)+ncity-2),ncity) ..and the one preceding n(1).. n(6)=1+mod(n(2),ncity) ..and the one following n(2). do 11j=1,6 ii=iorder(n(j)) Determine coordinates for the six cities in- volved. xx(j)=x(ii) yy(j)=y(ii) enddo 11 de=-alen(xx(2),xx(6),yy(2),yy(6)) Calculate the cost of disconnecting the path segment from n(1)ton(2), opening a space between n(3)andn(4), connect- ing the segment in the space, and con- necting n(5)ton(6).* -alen(xx(1),xx(5),yy(1),yy(5)) * -alen(xx(3),xx(4),yy(3),yy(4)) * +alen(xx(1),xx(3),yy(1),yy(3)) * +alen(xx(2),xx(4),yy(2),yy(4))* +alen(xx(5),xx(6),yy(5),yy(6)) return END SUBROUTINE trnspt(iorder,ncity,n) INTEGER ncity,iorder(ncity),n(6),MXCITYPARAMETER (MXCITY=1000) Maximum number of cities anticipated. This routine does the actual path transport, once metrophas approved. iorderis an input arrayoflength ncitygivingthe present itinerary. Thearray nhasasitssixelements the beginning n(1)and end n(2)of the path to be transported, the adjacent cities n(3) andn(4)between which the path is to be placed, and the cities n(5)andn(6)that precede and follow the path. n(4),n(5),an d n(6)arecalculated by subroutine trncst. On output, iorderis modified to reflect the movement of the path segment. INTEGER j,jj,m1,m2,m3,nn,jorder(MXCITY) m1=1+mod((n(2)-n(1)+ncity),ncity) Find number of cities from n(1)ton(2) 10.9SimulatedAnnealingMethods 443Sample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).m2=1+mod((n(5)-n(4)+ncity),ncity) ...and the number from n(4)ton(5) m3=1+mod((n(3)-n(6)+ncity),ncity) ...and the number from n(6)ton(3). nn=1 do11j=1,m1 jj=1+mod((j+n(1)-2),ncity) Copy the chosen segment. jorder(nn)=iorder(jj) nn=nn+1 enddo 11 do12j=1,m2 Thencopy thesegment from n(4)ton(5). jj=1+mod((j+n(4)-2),ncity) jorder(nn)=iorder(jj)nn=nn+1 enddo 12 do13j=1,m3 Finally, the segment from n(6)ton(3). jj=1+mod((j+n(6)-2),ncity)jorder(nn)=iorder(jj) nn=nn+1 enddo 13 do14j=1,ncity iorder(j)=jorder(j) Copy jorderback into iorder. enddo 14 return END SUBROUTINE metrop(de,t,ans) REAL de,tLOGICAL ans C USES ran3 Metropolis algorithm. ansisa logicalvariable that issues averdict on whether to accept a reconfigurationthatleadstoachange deintheobjectivefunction e.I fde<0,ans=.true. , whileif de>0,ansisonly .true.withprobability exp(-de/t) ,where tisatemperature determined by the annealing schedule. INTEGER jdumREAL ran3SAVE jdum DATA jdum /1/ ans=(de.lt.0.0).or.(ran3(jdum).lt.exp(-de/t))return END ContinuousMinimizationby Simulated Annealing The basic ideas of simulated annealing are also applicable to optimization problems with continuous N-dimensional control spaces, e.g., finding the (ideally, global) minimum of some function f(x), in the presence of many local minima, wherexis an N-dimensionalvector. The four elements requiredby the Metropolis procedure are now as follows: The value of fis the objective function. The system state is the point x. The control parameter Tis, as before, something like a temperature,withanannealingschedulebywhichitisgraduallyreduced. Andthere must be a generatorof randomchanges in the con figuration,that is, a procedurefor taking a random step from xtox+∆x. 444 Chapter10. MinimizationorMaximizationofFunctionsSample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).Thelastoftheseelementsisthemostproblematical. Theliteraturetodate [7-10] describes several different schemes for choosing ∆x, none of which, in our view, inspire complete con fidence. The problem is one of ef ficiency: A generator of random changes is inef ficient if,when local downhill moves exist , it nevertheless almost always proposes an uphill move. A good generator, we think, should notbecomeinef ficientinnarrowvalleys;norshoulditbecomemoreandmoreinef ficient as convergence to a minimum is approached. Except possibly for [7], all of the schemes that we have seen are inef ficient in one or both of these situations. Ourownwayofdoingsimulatedannealingminimizationoncontinuouscontrol spacesistouseamodi ficationofthedownhillsimplexmethod( §10.4). Thisamounts to replacing the single point xas a description of the system state by a simplex of N+1points. The “moves”are the same as described in §10.4, namely re flections, expansions,and contractionsof the simplex. The implementationof the Metropolisprocedure is slightly subtle: We adda positive, logarithmically distributed random variable, proportional to the temperature T, to the stored function value associated with every vertex of the simplex, and we subtracta similar random variable from the function value of every new point that is tried as a replacement point. Like the ordinaryMetropolisprocedure,this methodalways accepts atruedownhillstep,but sometimesaccepts anuphillone. In thelimit T→0,this algorithmreducesexactly to the downhill simplex method and converges to a local minimum. Atafinitevalueof T,thesimplexexpandstoascalethatapproximatesthesize of the regionthat can be reachedat this temperature,andthen executesa stochastic, tumblingBrownianmotionwithinthatregion,samplingnew,approximatelyrandom, points as it does so. The ef ficiency with which a region is explored is independent of its narrowness (for an ellipsoidal valley, the ratio of its principal axes) and orientation. If the temperature is reduced suf ficiently slowly, it becomes highly likely that the simplex will shrink into that region containing the lowest relative minimum encountered. As in all applications of simulated annealing, there can be quite a lot of problem-dependent subtlety in the phrase “sufficiently slowly ”; success or failure is quite often determined by the choice of annealing schedule. Here are some possibilities worth trying: •Reduce Tto(1−/epsilon1)Tafter every mmoves, where /epsilon1/mis determined by experiment. •Budget a total of Kmoves, and reduce Tafter every mmoves to a value T=T 0(1−k/K )α,where kis thecumulativenumberofmovesthusfar, andαis aconstant,say1,2,or4. Theoptimalvaluefor αdependsonthe statistical distributionof relative minimaof various depths. Largervalues ofαspend more iterations at lower temperature. •Afterevery mmoves,set Ttoβtimes f1−fb,where βisanexperimentally determinedconstantoforder1, f1is thesmallest functionvaluecurrently represented in the simplex, and fbis the best function ever encountered. However, never reduce Tby more than some fraction γat a time. Anotherstrategicquestioniswhethertodoanoccasional restart,whereavertex of the simplex is discardedin favor of the “best-ever”point. (You must be sure that thebest-everpointisnotcurrentlyinthesimplexwhenyoudothis!) Wehavefound problemsforwhich restarts —everytime the temperaturehas decreasedbya factor of 3, say —are highly bene ficial; we have found other problems for which restarts 10.9SimulatedAnnealingMethods 445Sample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).have no positive, or a somewhat negative, effect. Youshouldcomparethefollowingroutine, amebsa,withitscounterpart amoeba in§10.4. Note that the argument iteris used in a somewhat differentmanner. SUBROUTINE amebsa(p,y,mp,np,ndim,pb,yb,ftol,funk,iter,temptr) INTEGER iter,mp,ndim,np,NMAXREAL ftol,temptr,yb,p(mp,np),pb(np),y(mp),funkPARAMETER (NMAX=200) EXTERNAL funk C USES amotsa,funk,ran1 Multidimensional minimization of the function funk(x)where x(1:ndim) is a vector in ndimdimensions, by simulated annealing combined with the downhill simplex method of Nelder and Mead. The input matrix p(1..ndim+1,1..ndim) hasndim+1rows, each an ndim-dimensional vector which is a vertex of the starting simplex. Also input is the vector y(1:ndim+1) ,whosecomponentsmustbepre-initializedtothevaluesof funkevaluatedat thendim+1vertices(rows) of p;ftol,thefractionalconvergence tolerance tobeachieved in the function value for an early return; iter,a n d temptr. The routine makes iter function evaluations at an annealing temperature temptr, then returns. You should then decrease temptraccording to your annealing schedule, reset iter, and call the routine again(leavingotherargumentsunalteredbetween calls). If iterisreturnedwithapositive value, then early convergence and return occurred. If you initialize ybto avery large value on the firstcall, then ybandpb(1:ndim) will subsequently return the best function value and point ever encountered (even if it is no longer a point in the simplex). INTEGER i,idum,ihi,ilo,j,m,nREAL rtol,sum,swap,tt,yhi,ylo,ynhi,ysave,yt,ytry,psum(NMAX), * amotsa,ran1 COMMON /ambsa/ tt,idum tt=-temptr 1d o 12n=1,ndim Enter here when starting or after overall contraction. sum=0. Recompute psum. do11m=1,ndim+1 sum=sum+p(m,n) enddo 11 psum(n)=sum enddo 12 2 ilo=1 Enterhereafterchangingasinglepoint. Findwhichpoint isthehighest(worst),next-highest,andlowest(best). ihi=2 ylo=y(1)+tt*log(ran1(idum)) Wheneverwe“lookat”avertex,itgetsarandomthermal fluctuation. ynhi=ylo yhi=y(2)+tt*log(ran1(idum)) if (ylo.gt.yhi) then ihi=1ilo=2ynhi=yhi yhi=ylo ylo=ynhi endif do 13i=3,ndim+1 Loop over the points in the simplex. yt=y(i)+tt*log(ran1(idum)) More thermal fluctuations. if(yt.le.ylo) then ilo=i ylo=yt endifif(yt.gt.yhi) then ynhi=yhi ihi=iyhi=yt else if(yt.gt.ynhi) then ynhi=yt endif enddo 13 rtol=2.*abs(yhi-ylo)/(abs(yhi)+abs(ylo)) 446 Chapter10. MinimizationorMaximizationofFunctionsSample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).Compute the fractional range from highest to lowest and return if satisfactory. if (rtol.lt.ftol.or.iter.lt.0) then If returning, put best point and value in slot 1. swap=y(1) y(1)=y(ilo)y(ilo)=swap do 14n=1,ndim swap=p(1,n)p(1,n)=p(ilo,n)p(ilo,n)=swap enddo 14 return endifiter=iter-2 Beginanewiteration. Firstextrapolatebyafactor −1through thefaceofthesimplexacross from the high point, i.e., reflect the simplex from the high point. ytry=amotsa(p,y,psum,mp,np,ndim,pb,yb,funk,ihi,yhi,-1.0) if (ytry.le.ylo) then Gives a result better than the best point, so try an additional extrapolation by a factor 2. ytry=amotsa(p,y,psum,mp,np,ndim,pb,yb,funk,ihi,yhi,2.0) else if (ytry.ge.ynhi) then The reflected point isworse thanthe second-highest, solookforanintermediate lowerpoint, i.e., do a one-dimensional contraction. ysave=yhi ytry=amotsa(p,y,psum,mp,np,ndim,pb,yb,funk,ihi,yhi,0.5) if (ytry.ge.ysave) then Can’tseemtogetridofthathighpoint. Better contract around the lowest (best) point. do 16i=1,ndim+1 if(i.ne.ilo)then do15j=1,ndim psum(j)=0.5*(p(i,j)+p(ilo,j))p(i,j)=psum(j) enddo 15 y(i)=funk(psum) endif enddo 16 iter=iter-ndimgoto 1 endif else iter=iter+1 Correct the evaluation count. endifgoto 2 END FUNCTION amotsa(p,y,psum,mp,np,ndim,pb,yb,funk,ihi,yhi,fac) INTEGER ihi,mp,ndim,np,NMAX REAL amotsa,fac,yb,yhi,p(mp,np),pb(np),psum(np),y(mp),funkPARAMETER (NMAX=200)EXTERNAL funk C USES funk,ran1 Extrapolates by a factor facthrough the face of the simplex across from the high point, tries it, and replaces the high point if the new point is better. INTEGER idum,j REAL fac1,fac2,tt,yflu,ytry,ptry(NMAX),ran1COMMON /ambsa/ tt,idumfac1=(1.-fac)/ndim fac2=fac1-fac do 11j=1,ndim ptry(j)=psum(j)*fac1-p(ihi,j)*fac2 enddo 11 10.9SimulatedAnnealingMethods 447Sample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).ytry=funk(ptry) if (ytry.le.yb) then Save the best-ever. do12j=1,ndim pb(j)=ptry(j) enddo 12 yb=ytry endifyflu=ytry-tt*log(ran1(idum)) Weaddedathermalfluctuationtoallthecurrentvertices, but wesubtractit here, so as to give the simplex a thermal Brownian motion: It likesto accept any suggested change.if (yflu.lt.yhi) then y(ihi)=ytry yhi=yfludo 13j=1,ndim psum(j)=psum(j)-p(ihi,j)+ptry(j) p(ihi,j)=ptry(j) enddo 13 endif amotsa=yflu returnEND There is not yet enough practical experience with the method of simulated annealing to say de finitively what its future place among optimization methods will be. The method has several extremely attractive features, rather unique when compared with other optimization techniques. First, it is not “greedy,”in the sense that it is not easily fooled by the quick payoffachievedby fallinginto unfavorablelocal minima. Providedthat suf ficiently general recon figurations are given, it wanders freely among local minima of depth less than about T.A s Tis lowered, the number of such minima qualifying for frequent visits is gradually reduced. Second, con figuration decisions tend to proceed in a logical order. Changes thatcause thegreatestenergydifferencesaresifted overwhenthe controlparameter Tis large. These decisions become more permanent as Tis lowered, and attention thenshiftsmoretosmallerre finementsinthesolution. Forexample,inthetraveling salesman problem with the Mississippi River twist, if λis large, a decision to cross the Mississippi only twice is made at high T, while the speci fic routes on each side of the river are determined only at later stages. The analogies to thermodynamics may be pursued to a greater extent than we have done here. Quantities analogous to speci fic heat and entropy may be de fined, and these can be useful in monitoring the progress of the algorithm towards an acceptable solution. Information on this subject is found in [1]. CITED REFERENCES AND FURTHER READING: Kirkpatrick, S., Gelatt, C.D., and Vecchi, M.P. 1983, Science, vol. 220, pp. 671–680. [1] Kirkpatrick, S. 1984, Journal of Statistical Physics , vol. 34, pp. 975–986. [2] Vecchi, M.P. andKirkpatrick, S. 1983, IEEE Transactions onComputer AidedDesign ,vol. CAD- 2, pp. 215–222. [3] Otten, R.H.J.M., and van Ginneken, L.P.P.P. 1989, The Annealing Algorithm (Boston: Kluwer) [contains many references to the literature]. [4] Metropolis,N.,Rosenbluth,A.,Rosenbluth,M.,TellerA.,andTeller,E.1953, JournalofChemical Physics, vol. 21, pp. 1087–1092. [5] Lin, S. 1965, Bell System Technical Journal , vol. 44, pp. 2245–2269. [6] Vanderbilt,D.,andLouie,S.G.1984, JournalofComputationalPhysics ,vol.56,pp.259–271.[7] Bohachevsky,I.O.,Johnson,M.E.,andStein,M.L.1986, Technometrics ,vol.28,pp.209–217.[8] 448 Chapter10. MinimizationorMaximizationofFunctionsSample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).Corana, A., Marchesi, M., Martini, C., andRidella,S. 1987, ACM Transactions on Mathematical Software, vol. 13, pp. 262–280. [9] B´elisle, C.J.P., Romeijn, H.E., and Smith, R.L. 1990, Technical Report 90–25, Department of Industrial and Operations Engineering, University of Michigan, submitted to Mathematical Programming . [10] Christofides, N., Mingozzi, A., Toth, P., and Sandi, C. (eds.) 1979, Combinatorial Optimization (London andNew York: Wiley-Interscience) [not simulated annealing, butother topics andalgorithms].