Simulated Annealing for Traveling Salesman Problem - Artificial Intelligence C Program

/*C Implementation: (For speed it was necessary to reprogram the algorithm in C, using Mathematica
only to display the results. The following code managed to execute about 10,000 iterations per minute
on the 100-city problem using a 486-66 processor.) */
/* C code for Simulating Annealing on TSP -- David Bookstaber */
/* Artificial Intelligence */

#include <stdio.h>
#include <stdlib.h>
#include <conio.h>
#include <math.h>
#define N 100 // Number of cities
#define K 70 // Dimension of state graph (range of interchange)
#define ITS 100000L // Iterations
#define PMELT 0.7 // Fraction of melting point for starting temperature
#define TARGT 0.01 // Fraction of melting point for ending temperature
#define STAGFACTOR 0.05 // Fraction of runs to allow stagnant before reheating

typedef struct
{
int list[N]; float fit;} state;
float city[N][2];
float distance(int a, int b) {
return(pow(pow(city[a][0]-city[b][0],2)+pow(city[a][1]-city[b][1],2),0.5));
}
float fitness(state x) {
int i;
float sum = distance(x.list[0],x.list[N-1]);
for(i = 1;i < N;i++) sum += distance(x.list[i],x.list[i-1]);
return(sum);
}
state iterate(float temp,state x) {
int i, pt, sh;
state y = x;
pt = rand() % N;
sh = (pt+(rand() % K)+1) % N;
y.list[pt] = y.list[pt]^y.list[sh];
y.list[sh] = y.list[sh]^y.list[pt];
y.list[pt] = y.list[pt]^y.list[sh];
y.fit = fitness(y);
if(y.fit < x.fit) {
return(y);
}
else if((float)rand()/(1.0*RAND_MAX) < exp(-1.0*(y.fit-x.fit)/temp))
return(y);
else
return(x);
}
void main() {
int i, j, k, n;
long l, optgen;
double minf, maxf, range, Dtemp;
state x, optimum;
FILE *fp;
clrscr();
srand(1);
/* Initialization of city grid and state list */
for(i = 0,k = 0,n = sqrt(N);i < n;i++) {
for(j = 0;j < n;j++,k=n*i+j) {
city[k][0] = i; city[k][1] = j;
x.list[k] = k;
}
}
/* Randomization of state list--requires N Log[N] "shuffles" */
for(i = 0,k = rand()%(N-1)+1;i < N*log(N);i++,k = rand()%(N-1)+1) {
x.list[0] = x.list[0]^x.list[k];
x.list[k] = x.list[k]^x.list[0];
x.list[0] = x.list[0]^x.list[k];
}

/* Sample state space with 1% of runs to determine temperature schedule */
for(i = 0,maxf = 0,minf = pow(10,10),x.fit=fitness(x);i < max(0.01*N,2);i++) {
x = iterate(pow(10,10),x);
minf = (x.fit < minf) ? x.fit : minf;
maxf = (x.fit > maxf) ? x.fit : maxf;
}
range = (maxf - minf)*PMELT;
Dtemp = pow(TARGT,1.0/ITS);
/* Simulate Annealing */
for(optgen = l = 1,optimum.fit = x.fit;l < ITS;l++) {
x = iterate(range*pow(Dtemp,l),x);
if(x.fit < optimum.fit) {
optimum = x;
optgen = l;
}
/* Reheat if stagnant */
if(l-optgen == STAGFACTOR*ITS) Dtemp = pow(Dtemp,.05*l/ITS);
/* Graphics */
printf("\n\n\n\n\n\n\n");
gotoxy(1,1); printf("Iteration: %ld\t",l);
gotoxy(1,2); printf("Fitness: %f\t\tTemp: %f\t\t",x.fit,range*pow(Dtemp,l));
gotoxy(1,3); printf("Current Optimum %f found on %ld\t\t",optimum.fit,optgen);
gotoxy(1,4); printf("Global Optimum is %d",N);
gotoxy(1,5); printf("Sample Range: %f\tTemp decrement: %f",range,Dtemp);
/* End Graphics */
}
/* Output Solution */
fp = fopen("\\DOCS\\SA.OUT","wt");
fputc('{',fp);
for(i = 0;i < N;i++)
fprintf(fp,"%d,",optimum.list[i]);
fprintf(fp,"%f",optimum.fit);
fputc('}',fp);
fclose(fp);
}