basicant_段海滨.cpp

基本蚁群算法的源代码 经过一些改正 改正了原来的一些错误 可以运行

//*****************************************//
//   Basic Ant Colony Algorithm for TSP    //
//                c++ Version              //
// 摘自段海滨著<<蚁群算法原理及其应用>>    //
// 北京:科学出版社,2005.12,p421-p426       //
// ISBN 7-03-016204-8                      //
// 已经更正了缺少的函数initparameter()     //
// 并更正了其它的错误,可以直接运行         //
// Visual C++ 6.0 调试通过                 //
// 数据文件为tsp.dat                       //
// 每一行格式: num  x  y                   //
// 注意: N的值必须等于tsp.dat中的城市数量  //
//*****************************************//

#include <iostream.h>
#include <fstream.h>
#include <math.h>
#include <time.h>
#include <conio.h>
#include <float.h>        // 需要DBL_MAX,可看成无穷大
#include <stdlib.h>
#include <iomanip.h>

#define N 30             // city size,城市数目修改时要同时修改数据文件
#define M 30             // ant number,蚂蚁数量就随你喜欢

double inittao = 1;
double tao[N][N];
double detatao[N][N];
double distance[N][N];
double yita[N][N];
int tabu[M][N];
int route[M][N];
double solution[M];
int BestRoute[N];
double BestSolution = DBL_MAX;
double alfa, beta, rou, Q;
int NcMax;

void initparameter();                  // initialize the parameters of basic ACA
double EvalueSolution( int *a);        // evaluate the solution of TSP,
                                       // and calculate the length of path
void InCityXY(double x[], double y[]); // input the nodes' coordinates of TSP

void main()
{
	int NC = 0;
	initparameter();
	double x[N], y[N];
	InCityXY(x,y);
	int i,j;
	for(i = 0; i < N; i++)
	{
		for(j = i + 1; j < N; j++)
		{
			distance[j][i] = sqrt((x[i]-x[j])*(x[i]-x[j])+(y[i]-y[j])*(y[i]-y[j]));
			distance[i][j] = distance[j][i];
		}
	}
	// calculate the heuristic parameters
	for(i = 0; i < N; i++)
	{
		for(j = 0; j < N; j++)
		{
			tao[i][j] = inittao;
			if(j != i)
				yita[i][j] = 100 / (distance[i][j] + 50 - int(distance[i][j]) % 50);
		}
	}
	for(int k = 0; k < M; k++)
	{
		for(i = 0; i < N; i++)
			route[k][i] = -1;
	}

	srand( time(NULL) );
	for(k = 0; k < M; k++)
	{
		route[k][0] = (k + N) % N;
		//route[k][0] = rand() % N;
		tabu[k][route[k][0]] = 1;
	}
	// each ant try to find the optimal path
	do
	{
		int s = 1;
		double partsum;
		double pper;
		double drand;
		// ant choose one whole path
		while (s < N)
		{
			for(int k = 0; k < M; k++)
			{
				drand = rand() % 3000 / 3000.0;
				partsum = 0.0; pper = 0.0;
				for(int j = 0; j < N; j++)
				{
					if(tabu[k][j] == 0)
						partsum += pow(tao[route[k][s-1]][j],alfa) * pow(yita[route[k][s-1]][j],beta);
				}
				for(j = 0; j < N; j++)
				{
					if(tabu[k][j] == 0)
						pper += pow(tao[route[k][s-1]][j],alfa) * pow(yita[route[k][s-1]][j],beta)/partsum;
					if(drand < pper) break;   // 找到了下一个城市j,跳出
				}
				tabu[k][j] = 1;
				route[k][s] = j;
			}
			s++;
		}
		// the pheromone is updated
		for(i = 0;i < N; i++)
		{
			for(int j = 0; j < N; j++)
			{
				detatao[i][j] = 0;
			}
		}
		for(int k = 0; k < M; k++)
			solution[k] = EvalueSolution(route[k]);
		for(k = 1; k < M; k++)
		{
			if(solution[k] < BestSolution)
			{
				BestSolution = solution[k];
				for(s = 0; s < N; s++)
					BestRoute[s] = route[k][s];
			}
		}
		for(k = 0; k < M; k++)
		{
			for(s = 0; s < N-1; s++)
				detatao[route[k][s]][route[k][s]] += Q / solution[k];
			detatao[route[k][N-1]][route[k][0]] += Q / solution[k];
		}
		for(i = 0; i < N; i++)
		{
			for(j = 0; j < N; j++)
			{
				tao[i][j] = rou * tao[i][j] + detatao[i][j];
				if(tao[i][j] < 0.00001)
					tao[i][j] = 0.00001;
				if(tao[i][j] > 20)
					tao[i][j] = 20;
			}
		}
		for(k = 0; k < M; k++)
		{
			for(j = 1; j < N; j++)
			{
				tabu[k][route[k][j]] = 0;
				route[k][j] = -1;
			}
		}
		NC++;
	}while(NC < NcMax);
	// output the calculating results
	fstream result;
	result.open("Optimal_Results.dat",ios::app);
	if(!result)
	{
		cout << " Can't open the 'Optimal_Results.dat' file\n";
		abort();
	}
	
	result << endl << endl;
	result << endl << " alfa = " << alfa << ", beta = " << beta << endl;
	result << " Q = " << Q << ", The maximum iteration number of ACA is: " << NcMax << endl;
	result << " The evaporation parameter of ACA is: " << rou << endl;
	
	// BestSolution = EvalueSolution(BestRoute);
	result << endl << endl << " 最佳路径长度: " << BestSolution << endl;
	result << endl << " 最佳路径如下所示: " << endl;
	for(i = 0; i < N; i++)
	{
		if(i % 20 == 0) result << endl << " ";
		result << BestRoute[i] + 1 << "  "; // 城市编号1,2,3,...,N
	}
	result << endl;
	result.close();
}

//****************************************//
// initialize the parameters of basic ACA //
//****************************************//
void initparameter()
{
	cout << " 信息启发式因子alfa = "; cin >> alfa;
	cout << " 期望启发式因子beta = "; cin >> beta;
	cout << " 信息素强度Q = "; cin >> Q;
	cout << " 最大循环次数NcMax = "; cin >> NcMax;
	cout << " 信息素挥发系数rou = "; cin >> rou;
}

//***********************************//
// Function: Evaluation the solution //
//***********************************//
double EvalueSolution(int *a)
{
	double dist = 0.0;
	for(int i = 0; i < N; i++)
		dist += distance[*(a + (i + N) % N)][*(a + (i + N + 1) % N)];
	return dist;
}

//************************************//
// Function: Set the city coordinates //
//************************************//
void InCityXY(double x[],double y[])
{
	fstream inxyfile;
	inxyfile.open("tsp.dat",ios::in);
	if(!inxyfile)
	{
		cout << " Can't open the 'tsp.dat' file\n";
		abort();
	}
	char ch1,ch2;
	while(! inxyfile.eof())
	{
		inxyfile.get(ch1);
		if(ch1 >= '0') break;
		int i = 0, j = 0;
		x[0] = y[0] = 0;
		while(! inxyfile.eof())
		{
			inxyfile.get(ch1);
			if(ch1 >= '0' && ch1 <= '9')
			{
				ch2 = ch1;
				while(ch2 >= '0' && ch2 <= '9')
				{
					switch(i)
					{
					case 0: break;
					case 1: 
						x[j] = x[j] * 10 + (double(ch2) - 48);
						break;
					case 2: 
						y[j] = y[j] * 10 + (double(ch2) - 48);
						break;
					}
					inxyfile.get(ch2);
				}
				i = (++i) % 3;
				if(i == 0 && j < N-1)
				{
					j ++;
					x[j] = 0;
					y[j] = 0;
				}
			}
		}
	}
}