c10-9 (Numerical Recipes in C), страница 4

(You must be sure thatthe best-ever point is not currently in the simplex when you do this!) We have foundproblems for which restarts — every time the temperature has decreased by a factorof 3, say — are highly beneficial; we have found other problems for which restartshave no positive, or a somewhat negative, effect.You should compare the following routine, amebsa, with its counterpart amoebain §10.4. Note that the argument iter is used in a somewhat different manner.10.9 Simulated Annealing Methods453psum=vector(1,ndim);tt = -temptr;GET_PSUMfor (;;) {ilo=1;Determine which point is the highest (worst),ihi=2;next-highest, and lowest (best).ynhi=ylo=y[1]+tt*log(ran1(&idum));Whenever we “look at” a vertex, it getsyhi=y[2]+tt*log(ran1(&idum));a random thermal fluctuation.if (ylo > yhi) {ihi=1;ilo=2;ynhi=yhi;yhi=ylo;ylo=ynhi;}for (i=3;i<=mpts;i++) {Loop over the points in the simplex.yt=y[i]+tt*log(ran1(&idum));More thermal fluctuations.if (yt <= ylo) {ilo=i;ylo=yt;}if (yt > yhi) {ynhi=yhi;ihi=i;yhi=yt;} else if (yt > ynhi) {ynhi=yt;}}rtol=2.0*fabs(yhi-ylo)/(fabs(yhi)+fabs(ylo));Compute the fractional range from highest to lowest and return if satisfactory.if (rtol < ftol || *iter < 0) {If returning, put best point and value inswap=y[1];slot 1.y[1]=y[ilo];y[ilo]=swap;for (n=1;n<=ndim;n++) {swap=p[1][n];p[1][n]=p[ilo][n];p[ilo][n]=swap;}break;}*iter -= 2;Begin a new iteration.

First extrapolate by a factor −1 through the face of the simplexacross from the high point, i.e., reflect the simplex from the high point.ytry=amotsa(p,y,psum,ndim,pb,yb,funk,ihi,&yhi,-1.0);if (ytry <= ylo) {Gives a result better than the best point, so try an additional extrapolation by afactor of 2.ytry=amotsa(p,y,psum,ndim,pb,yb,funk,ihi,&yhi,2.0);} else if (ytry >= ynhi) {The reflected point is worse than the second-highest, so look for an intermediateSample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5)Copyright (C) 1988-1992 by Cambridge University Press.Programs Copyright (C) 1988-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 copying of machinereadable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes books,diskettes, or CDROMsvisit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to trade@cup.cam.ac.uk (outside North America).crease temptr according to your annealing schedule, reset iter, and call the routine again(leaving other arguments unaltered between calls). If iter is returned with a positive value,then early convergence and return occurred. If you initialize yb to a very large value on the firstcall, then yb and pb[1..ndim] will subsequently return the best function value and point everencountered (even if it is no longer a point in the simplex).{float amotsa(float **p, float y[], float psum[], int ndim, float pb[],float *yb, float (*funk)(float []), int ihi, float *yhi, float fac);float ran1(long *idum);int i,ihi,ilo,j,m,n,mpts=ndim+1;float rtol,sum,swap,yhi,ylo,ynhi,ysave,yt,ytry,*psum;454Chapter 10.Minimization or Maximization of Functions}free_vector(psum,1,ndim);}#include <math.h>#include "nrutil.h"extern long idum;extern float tt;Defined and initialized in main.Defined in amebsa.float amotsa(float **p, float y[], float psum[], int ndim, float pb[],float *yb, float (*funk)(float []), int ihi, float *yhi, float fac)Extrapolates by a factor fac through the face of the simplex across from the high point, triesit, and replaces the high point if the new point is better.{float ran1(long *idum);int j;float fac1,fac2,yflu,ytry,*ptry;ptry=vector(1,ndim);fac1=(1.0-fac)/ndim;fac2=fac1-fac;for (j=1;j<=ndim;j++)ptry[j]=psum[j]*fac1-p[ihi][j]*fac2;ytry=(*funk)(ptry);if (ytry <= *yb) {Save the best-ever.for (j=1;j<=ndim;j++) pb[j]=ptry[j];*yb=ytry;}yflu=ytry-tt*log(ran1(&idum));We added a thermal fluctuation to all the currentif (yflu < *yhi) {vertices, but we subtract it here, so as to givey[ihi]=ytry;the simplex a thermal Brownian motion: It*yhi=yflu;likes to accept any suggested change.for (j=1;j<=ndim;j++) {psum[j] += ptry[j]-p[ihi][j];p[ihi][j]=ptry[j];}}free_vector(ptry,1,ndim);return yflu;}There is not yet enough practical experience with the method of simulatedannealing to say definitively what its future place among optimization methodsSample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5)Copyright (C) 1988-1992 by Cambridge University Press.Programs Copyright (C) 1988-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 copying of machinereadable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes books,diskettes, or CDROMsvisit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to trade@cup.cam.ac.uk (outside North America).lower point, i.e., do a one-dimensional contraction.ysave=yhi;ytry=amotsa(p,y,psum,ndim,pb,yb,funk,ihi,&yhi,0.5);if (ytry >= ysave) {Can’t seem to get rid of that high point.for (i=1;i<=mpts;i++) {Better contract around the lowestif (i != ilo) {(best) point.for (j=1;j<=ndim;j++) {psum[j]=0.5*(p[i][j]+p[ilo][j]);p[i][j]=psum[j];}y[i]=(*funk)(psum);}}*iter -= ndim;GET_PSUMRecompute psum.}} else ++(*iter);Correct the evaluation count.10.9 Simulated Annealing Methods455CITED 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. and Kirkpatrick, S. 1983, IEEE Transactions on Computer Aided Design, vol. CAD2, 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., Teller A., and Teller, E. 1953, Journal of ChemicalPhysics, vol.

21, pp. 1087–1092. [5]Lin, S. 1965, Bell System Technical Journal, vol. 44, pp. 2245–2269. [6]Vanderbilt, D., and Louie, S.G. 1984, Journal of Computational Physics, vol. 56, pp. 259–271. [7]Bohachevsky, I.O., Johnson, M.E., and Stein, M.L. 1986, Technometrics, vol. 28, pp. 209–217. [8]Corana, A., Marchesi, M., Martini, C., and Ridella, S. 1987, ACM Transactions on MathematicalSoftware, vol. 13, pp. 262–280.

[9]Bélisle, C.J.P., Romeijn, H.E., and Smith, R.L. 1990, Technical Report 90–25, Department ofIndustrial and Operations Engineering, University of Michigan, submitted to MathematicalProgramming. [10]Christofides, N., Mingozzi, A., Toth, P., and Sandi, C. (eds.) 1979, Combinatorial Optimization(London and New York: Wiley-Interscience) [not simulated annealing, but other topics andalgorithms].Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5)Copyright (C) 1988-1992 by Cambridge University Press.Programs Copyright (C) 1988-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 copying of machinereadable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes books,diskettes, or CDROMsvisit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to trade@cup.cam.ac.uk (outside North America).will be.

The method has several extremely attractive features, rather unique whencompared with other optimization techniques.First, it is not “greedy,” in the sense that it is not easily fooled by the quickpayoff achieved by falling into unfavorable local minima. Provided that sufficientlygeneral reconfigurations are given, it wanders freely among local minima of depthless than about T . As T is lowered, the number of such minima qualifying forfrequent visits is gradually reduced.Second, configuration decisions tend to proceed in a logical order.

Changesthat cause the greatest energy differences are sifted over when the control parameterT is large. These decisions become more permanent as T is lowered, and attentionthen shifts more to smaller refinements in the solution. For example, in the travelingsalesman problem with the Mississippi River twist, if λ is large, a decision to crossthe Mississippi only twice is made at high T , while the specific routes on each sideof the river are determined only at later stages.The analogies to thermodynamics may be pursued to a greater extent than wehave done here.

Quantities analogous to specific heat and entropy may be defined,and these can be useful in monitoring the progress of the algorithm towards anacceptable solution. Information on this subject is found in [1]..