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

iorder[1..ncity]is an input array giving the present itinerary. The array n has as its six elements the beginningn[1] and end n[2] of the path to be transported, the adjacent cities n[3] and n[4] betweenwhich the path is to be placed, and the cities n[5] and n[6] that precede and follow the path.n[4], n[5], and n[6] are calculated by function trncst.

On output, iorder is modified toreflect the movement of the path segment.{int m1,m2,m3,nn,j,jj,*jorder;jorder=ivector(1,ncity);m1=1 + ((n[2]-n[1]+ncity) % ncity);m2=1 + ((n[5]-n[4]+ncity) % ncity);m3=1 + ((n[3]-n[6]+ncity) % ncity);nn=1;for (j=1;j<=m1;j++) {jj=1 + ((j+n[1]-2) % ncity);jorder[nn++]=iorder[jj];}for (j=1;j<=m2;j++) {jj=1+((j+n[4]-2) % ncity);jorder[nn++]=iorder[jj];}for (j=1;j<=m3;j++) {jj=1 + ((j+n[6]-2) % ncity);jorder[nn++]=iorder[jj];}for (j=1;j<=ncity;j++)iorder[j]=jorder[j];Find number of cities from n[1] to n[2]...and the number from n[4] to n[5]...and the number from n[6] to n[3].Copy the chosen segment.Then copy the segment from n[4] ton[5].Finally, the segment from n[6] to n[3].Copy jorder back into iorder.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).float trncst(float x[], float y[], int iorder[], int ncity, int n[])This routine returns the value of the cost function for a proposed path segment transport. ncityis the number of cities, and arrays x[1..ncity] and y[1..ncity] give the city coordinates.iorder[1..ncity] is an array giving the present itinerary.

The first three elements of arrayn give the starting and ending cities of the path to be transported, and the point among theremaining cities after which it is to be inserted. On output, de is the cost of the change. Theactual transport is not performed by this routine.{float xx[7],yy[7],de;int j,ii;10.9 Simulated Annealing Methods451free_ivector(jorder,1,ncity);}#include <math.h>return de < 0.0 || ran3(&gljdum) < exp(-de/t);}Continuous Minimization by Simulated AnnealingThe basic ideas of simulated annealing are also applicable to optimizationproblems 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,where x is an N -dimensional vector.

The four elements required by the Metropolisprocedure are now as follows: The value of f is the objective function. Thesystem state is the point x. The control parameter T is, as before, something like atemperature, with an annealing schedule by which it is gradually reduced. And theremust be a generator of random changes in the configuration, that is, a procedure fortaking a random step from x to x + ∆x.The last of these elements is the most problematical. The literature to date [7-10]describes several different schemes for choosing ∆x, none of which, in our view,inspire complete confidence.

The problem is one of efficiency: A generator ofrandom changes is inefficient if, when local downhill moves exist, it neverthelessalmost always proposes an uphill move. A good generator, we think, should notbecome inefficient in narrow valleys; nor should it become more and more inefficientas convergence to a minimum is approached. Except possibly for [7], all of theschemes that we have seen are inefficient in one or both of these situations.Our own way of doing simulated annealing minimization on continuous controlspaces is to use a modification of the downhill simplex method (§10.4).

This amountsto replacing the single point x as a description of the system state by a simplex ofN + 1 points. The “moves” are the same as described in §10.4, namely reflections,expansions, and contractions of the simplex. The implementation of the Metropolisprocedure is slightly subtle: We add a positive, logarithmically distributed randomvariable, proportional to the temperature T , to the stored function value associatedwith every vertex of the simplex, and we subtract a similar random variable fromthe function value of every new point that is tried as a replacement point.

Like theordinary Metropolis procedure, this method always accepts a true downhill step, butSample 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).int metrop(float de, float t)Metropolis algorithm.

metrop returns a boolean variable that issues a verdict on whetherto accept a reconfiguration that leads to a change de in the objective function e. If de<0,metrop = 1 (true), while if de>0, metrop is only true with probability exp(-de/t), wheret is a temperature determined by the annealing schedule.{float ran3(long *idum);static long gljdum=1;452Chapter 10.Minimization or Maximization of Functions#include <math.h>#include "nrutil.h"#define GET_PSUM \for (n=1;n<=ndim;n++) {\for (sum=0.0,m=1;m<=mpts;m++) sum += p[m][n];\psum[n]=sum;}extern long idum;Defined and initialized in main.float tt;Communicates with amotsa.void amebsa(float **p, float y[], int ndim, float pb[], float *yb, float ftol,float (*funk)(float []), int *iter, float temptr)Multidimensional minimization of the function funk(x) where x[1..ndim] is a vector inndim dimensions, by simulated annealing combined with the downhill simplex method of Nelderand Mead.

The input matrix p[1..ndim+1][1..ndim] has ndim+1 rows, each an ndimdimensional vector which is a vertex of the starting simplex. Also input are the following: thevector y[1..ndim+1], whose components must be pre-initialized to the values of funk evaluated at the ndim+1 vertices (rows) of p; ftol, the fractional convergence tolerance to beachieved in the function value for an early return; iter, and temptr. The routine makes iterfunction evaluations at an annealing temperature temptr, then returns. You should then de-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).sometimes accepts an uphill one. In the limit T → 0, this algorithm reduces exactlyto the downhill simplex method and converges to a local minimum.At a finite value of T , the simplex expands to a scale that approximates the sizeof the region that can be reached at this temperature, and then executes a stochastic,tumbling Brownian motion within that region, sampling new, approximately random,points as it does so. The efficiency with which a region is explored is independentof its narrowness (for an ellipsoidal valley, the ratio of its principal axes) andorientation.

If the temperature is reduced sufficiently slowly, it becomes highlylikely that the simplex will shrink into that region containing the lowest relativeminimum encountered.As in all applications of simulated annealing, there can be quite a lot ofproblem-dependent subtlety in the phrase “sufficiently slowly”; success or failureis quite often determined by the choice of annealing schedule. Here are somepossibilities worth trying:• Reduce T to (1 − )T after every m moves, where /m is determinedby experiment.• Budget a total of K moves, and reduce T after every m moves to a valueT = T0 (1 − k/K)α , where k is the cumulative number of moves thus far,and α is a constant, say 1, 2, or 4.

The optimal value for α depends on thestatistical distribution of relative minima of various depths. Larger valuesof α spend more iterations at lower temperature.• After every m moves, set T to β times f1 −fb , where β is an experimentallydetermined constant of order 1, f1 is the smallest function value currentlyrepresented in the simplex, and fb is the best function ever encountered.However, never reduce T by more than some fraction γ at a time.Another strategic question is whether to do an occasional restart, where a vertexof the simplex is discarded in favor of the “best-ever” point.