📄 acs.cpp
字号:
/*
ANT-CYCLE ALGORITHM FOR TSP
File: acs.c
Author: ehui928
Purpose: implementation of acs.h
Date: 2007-01-18
Reference: Dorigo M --"The Ant System: Optimization by a colony of cooperating agents" IEEE 1996
*/
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>
#include <memory.h>
#include "acs.h"
#include <assert.h>
Acs::Acs()
{
int i,j;
seed = 21345678; /* seed to generate random number */
alpha = 1; /* intensity of trail */
beta = 4.5; /* visibility of trail */
rho = 0.5; /* 1-rho is the evaporate rate */
tao0 = 1e-6; /* initial pheromone */
//tao=0会更好
Q = 100; /* const used to update pheromone */
q_0 = 0.25; /* used to decide which transition rule to be choosed */
memset(distances,10000,(CITY_NUM+1)*(CITY_NUM+1)*sizeof(double));
for(i=0;i<CITY_NUM+1;i++)
{
for(j=0;j<CITY_NUM+1;j++) distances[i][j]=10000.0;
DataPoints[i].x=0;
DataPoints[i].y=0;
}
for(i=0;i<NC_MAX;i++)
{
bestPerCircle[i]=0.0;
}
gLog.SetFileName("result.log");
if ( (f_out = fopen("result.txt", "w+")) == NULL)
{
fprintf(stderr, "can not create file result.txt\n");
exit(2);
}
}
Acs::~Acs()
{
gLog.Close();
fclose(f_out);
}
void Acs::read_data(char *filename)//1.读数据
{
/*
read city corordinates i,x, y refers to each record in
input file----(city_num, coordinate1,coordinate2)
*/
int i;
double x, y;
int n;
FILE *fp;
if ( (fp = fopen(filename, "r")) == NULL)
{
fprintf(stderr, "can not open file %s\n", filename);
exit(1);
}
n = 1;
while (n < CITY_NUM+1)
{
/*
since x, y all double numbers, so here the
format must be %lf , not %f, otherwise data reading will be error
*/
fscanf(fp, "%d %lf %lf", &i, &x, &y);
DataPoints[i].x = x;
DataPoints[i].y = y;
n++;
}
fclose(fp);
}
void Acs::initAll()//初始化
{
caculate_distances();
}
/*
FUNCTION: find the global best tour of TSP problem
INPUT: none
RETURN: the length of global best tour
SIDE EFFECTS: print the global best tour
NC --------------------iterator time
i,j--------------------city index
k----------------------ant index
s----------------------tabu index
ants-------------------ants which used to find tours
ants[0] is not used, so as the L[0],best_tour[0] etc.
start------------------starting city
best_tour--------------best tour ever find by ANT_NUM ants
best_ant_index---------index of ant which found the best tour
L[ANT_NUM+1]-----------------tour lengths ANT_NUM ant found
best_len---------------record best tour length all ants find one iteration
global_best_len--------global best length
converge---------------test convergence
*/
double Acs::find_best_tour()
{
int NC = 0;
int i, j;
int k=1;
int s=1;
Ant ants[ANT_NUM+1];
int start = 1;
int best_tour[CITY_NUM+1];
int best_ant_index = 1;
double L[ANT_NUM+1];
double best_len = 10000;
int converge = 0;
/* initialize pheromone on city edges */
initialize_pheromone();
memset(best_tour,0,(CITY_NUM+1)*sizeof(int));
memset(L,0,(CITY_NUM+1)*sizeof(double));
global_best_len = 10000;
while (NC < NC_MAX)
{
s = 1; /* tabu list index */
/* put each ant to starting city */
for (k = 1; k <= ANT_NUM; k++)//每次所有蚂蚁都从第一个城市出发
{
initial_ant(&ants[k]);
ants[k].cur = start;
ants[k].allow[start] = start;
ants[k].tabu[s] = start;
}
/*for (k = 1; k <= ANT_NUM; k++)//蚂蚁从不同城市出发
{
initial_ant(&ants[k]);
ants[k].cur = k;
ants[k].allow[k] = 1;
ants[k].tabu[s] = k;
}*/
/* CITY_NUM-1 cities to choose */
for (i = 1; i < CITY_NUM; i++)//每一只蚂蚁走CITY_NUM-1个城市为一个循环
{
s++;
for (k = 1; k <= ANT_NUM; k++)
{
int t; /* save next city to move */
//t = choose_next_city(&ants[k]);
t = choose_next_city_LunPanDu(&ants[k]);
ants[k].tabu[s] = t;
ants[k].allow[t] = 1;//1为已经访问过
ants[k].cur = t;
}
}
/*
caculate the length of every tour that ants found,
and update the pheromone.
*/
for (k = 1; k <= ANT_NUM; k++)
{
/*
print tour each ant found,just for debug
*/
//print_ant_tour(ants[k]);
L[k] = caculate_tour_length(ants[k]);
if (L[k] < best_len)
{
best_len = L[k];
best_ant_index = k;
}
}
fprintf(f_out, "%d--circle\t%f\n\n",NC+1,best_len);//显示循环周期及最小路径
/* compute delta_tao[i][j] on every edge of cities */
for (i = 1; i <= CITY_NUM; i++)
{
for (j = 1; j <= CITY_NUM; j++)
{
for (k = 1; k <= ANT_NUM; k++)
{
if (is_in_tour(i, j, ants[k].tabu, CITY_NUM))
delta_tao[i][j] += Q/L[k];
}
}
}
/* update pheromone on every edge on a tour */
for (i = 1; i <= CITY_NUM; i++)
{
for (j = 1; j <= CITY_NUM; j++)
{
if(i != j) tao[i][j] = rho * tao[i][j] + delta_tao[i][j];
}
}
/* check if converge to one best tour */
if (fabs(best_len - global_best_len) < EPSLON)
{
converge++;
if (converge >= COVERG_NUM)
break;
}
else
if (best_len < global_best_len) /* save best tour */
{
converge = 0;
global_best_len = best_len;
memcpy(best_tour,ants[best_ant_index].tabu,(CITY_NUM+1)*sizeof(int));
}
//if(NC<Real_CIrcle)
bestPerCircle[NC]=best_len;
NC++;
reset_delta_tao();/* set delta_tao[i][j] to 0 */
} /* end while */
/* free resources that ants are occupying */
for (k = 1; k <= ANT_NUM; k++)
destroy_ant(&ants[k]);
print_best_tour(best_tour, CITY_NUM,global_best_len);
return global_best_len;
}
int Acs::is_in_tour(int i, int j, int tour[], int n)
/*
FUNCTION: check whether edge(i,j) is in tour
INPUT: i, j---city index
tour[]----tour an ant build
n---------size of tour[]
RETURN: if edge(i,j) is in tour, then return 1;
otherwise return 0.
EFFECTS: print the global best tour
*/
{
int s;
for (s = 1; s < n; s++)
{
if (i == tour[s] && j == tour[s+1])
return 1;
}
if (i == tour[s] && j == 1)
return 1;
return 0;
}
/* test if distances were caculated correctly */
void Acs::print_distance(double distances[CITY_NUM+1][CITY_NUM+1])
{
int i, j;
i = 1;
for (j = 1; j < CITY_NUM+1; j++)
printf("d[%d][%d] = %f\n", i, j, distances[i][j]);
}
void Acs::print_best_tour(int tour[], int n,double global_best_len)
/*
FUNCTION: print global best tour
INPUT: tour[]----global best tour
n---------size of tour[]
RETURN: none
EFFECTS: print the global best tour
*/
{
int i,index=0;
char buf[1000];
printf("Find best tour:");
for (i = 1; i <= n; i++)
{
printf("%d-", tour[i]);
index+=sprintf(&buf[index], "%d-", tour[i]);
}
printf("%d\r%f\n", tour[1],global_best_len);
index+=sprintf(&buf[index], "%d\r\nFind best tour: %f", tour[1],global_best_len);
this->gLog.Log(buf,index+1);
}
double Acs::Euclid_distance(Point p1, Point p2)
{
double a = p1.x - p2.x;
double b = p1.y - p2.y;
return (double)sqrt(a * a + b * b);
}
void Acs::caculate_distances()
{
int i, j;
for (i = 1; i <= CITY_NUM; i++)
for (j = 1; j <= CITY_NUM; j++)
{
if (i != j)
distances[i][j] = Euclid_distance(DataPoints[i], DataPoints[j]);
}
}
/* function to test if data has been read correctly */
void Acs::print_coordinates(Point *p)
{
int i;
printf("Cities Coordinates:\n");
for (i = 1; i < CITY_NUM+1; i++)
printf("%d %f %f\n", i, p[i].x, p[i].y);
}
void Acs::initial_ant(Ant *pant)
{
pant->cur = 0;
pant->length = 10000.0;
pant->tabu = (int *) malloc(sizeof(int) * (CITY_NUM+1));
pant->allow = (int *) malloc(sizeof(int) * (CITY_NUM+1));
memset(pant->tabu, 0, sizeof(int) * (CITY_NUM+1));
memset(pant->allow, 0, sizeof(int) * (CITY_NUM+1));
/* if you want the seed changes every time when the
program runs, then uncomment the following line
*/
//seed = (long int)time(NULL);
}
/*
FUNCTION: choose the next city according to transition rule
INPUT: the poiter to Ant---pant
RETURN: next city been choosed
*/
int Acs::choose_next_city(Ant *pant)
{
int i;//city index
int j;//next city to choose
int k = 1;//index of f_tmp
int s = 1;//index of prob[]
int n;//actual size of candidate
int candidate[CITY_NUM+1];//candidate cities can be chosen
double tmp[CITY_NUM+1];//record the values of tao
double max = -10000.0;//max values of tmp
double sum = 0;//sum of tmp
double prob[CITY_NUM+1];//probabilities of each candidate city to be choosen
double q; /* random generate q to decide using which transition rule */
for(i=0;i<= CITY_NUM; i++)
{
candidate[i]=0;
tmp[i]=0.0;
prob[i]=0.0;
}
for (i = 1; i <= CITY_NUM; i++)
{
if (!pant->allow[i])
{
/* city i is not visited */
assert( pant->cur >0 && pant->cur <= CITY_NUM );
assert(distances[pant->cur][i] >=0);
candidate[s++] = i;
tmp[k] = pow(tao[pant->cur][i], alpha) * pow(1/distances[pant->cur][i], beta);
if (tmp[k] >= max)
{
max = tmp[k];
j = i;
}
sum += tmp[k];
k++;
}
}
n = k - 1;
assert(sum > 0);
for (i = 1; i <= n; i++)
{
assert(tmp[i] > 0);
prob[i] = tmp[i] / sum;
}
/* generate a random double number between [0,1] */
q = ran01( &seed );
if ( (q_0 > 0.0) && (q < q_0) )
/*
choose next city according to
pheromone maximum multiply heuristic values
*/
{
return j;
}
else
/* choose a random city according to probabilities */
{
double rnd;
double partial_sum;
rnd = ran01( &seed );
i = 1;
partial_sum = prob[i];
while ( partial_sum <= rnd )
{
//assert(i <= CITY_NUM);
i++;
partial_sum += prob[i];
}
return candidate[i];
}
}
int Acs::choose_next_city_LunPanDu(Ant* pant)//2 使用轮盘赌法找下一个城市
{
int i;
int k = 1;
int s = 1;
int n;
int candidate[CITY_NUM+1];
double tmp[CITY_NUM+1];
double sum = 0;
double prob[CITY_NUM+1] = {0};
double obj[CITY_NUM+1] = {0};
for (i = 1; i <= CITY_NUM; i++)
{
if (!pant->allow[i])
{
/* city i is not visited */
candidate[s++] = i;
tmp[k] = pow(tao[pant->cur][i], alpha) * pow(1/distances[pant->cur][i], beta);
sum += tmp[k];
k++;
}
}
n = k - 1;
/* 轮盘赌选择
tmp[]原始数组，选择概率如{10，30，20，40} 此时sum=100
prob[]规格数据到【0，1】内{0.1，0.3，0.2，0.4}
obj[]概率的长度表示{0.1，0.4，0.6，1}
*/
for (i = 1; i <= n; i++)
{
prob[i] = tmp[i] / sum;
obj[i]=prob[i];
}
for (i = 2; i <= n; i++)
{
obj[i] += obj[i-1];
}
double rnd=0;
rnd = ((double)(rand()%1000))/1000.0;
for(i=1;i<=CITY_NUM;i++)
{
if(rnd < obj[i]) break;
}
return candidate[i];
}
void Acs::destroy_ant(Ant *pant)
{
free(pant->tabu);
free(pant->allow);
}
double Acs::caculate_tour_length(Ant ant)
{
int i;
double len = 0.0;
for (i = 1; i < CITY_NUM; i++)
len += distances[ant.tabu[i]][ant.tabu[i+1]];
len += distances[ant.tabu[CITY_NUM]][ant.tabu[1]];
return len;
}
void Acs::print_ant_tour(Ant ant)
/*
FUNCTION: print a tour an ant found
INPUT: Ant ant
RETURN: none
*/
{
int i;
for (i = 1; i < CITY_NUM+1; i++)
fprintf(f_out, "%d-", ant.tabu[i]);
fprintf(f_out, "%d\n", ant.tabu[1]);
/*int i,index=0;
char buf[1000];
for (i = 1; i < CITY_NUM+1; i++)
{
index+=sprintf(&buf[index], "%d-", ant.tabu[i]);
fprintf(f_out, "%d-", ant.tabu[i]);
}
index+=sprintf(&buf[index], "%d", ant.tabu[1]);
//this->gLog.Log(buf,index+1);
fprintf(f_out, "%d\n", ant.tabu[1]);*/
}
void Acs::initialize_pheromone()
/*
FUNCTION: initial pheromone on all edge(i,j) to tao0
INPUT: none
RETURN: none
*/
{
int i, j;
for (i = 1; i < CITY_NUM+1; i++)
for (j = 1; j < CITY_NUM+1; j++)
{
tao[i][j] = tao0;
delta_tao[i][j]=0.0;
}
}
/*
constants for a random number generator,
for details see numerical recipes in C
*/
double Acs::ran01( long *idum )
/*
FUNCTION: generate a random number that is uniformly distributed in [0,1]
INPUT: pointer to variable with the current seed
OUTPUT: random number uniformly distributed in [0,1]
(SIDE)EFFECTS: random number seed is modified (important, this has to be done!)
ORIGIN: numerical recipes in C
*/
{
#define IA 16807
#define IM 2147483647
#define AM (1.0/IM)
#define IQ 127773
#define IR 2836
long k;
double ans;
k =(*idum)/IQ;
*idum = IA * (*idum - k * IQ) - IR * k;
if (*idum < 0 ) *idum += IM;
ans = AM * (*idum);
return ans;
}
//FUNCTION: clear the delta_tao matrix, set it all to 0
void Acs::reset_delta_tao()
{
int i, j;
for (i = 1; i < CITY_NUM+1; i++)
for (j = 1; j < CITY_NUM+1; j++)
delta_tao[i][j] = 0.0;
}⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -