city.cpp

蚁群算法 控制台程序 基于类 不是基本的TSP模型的ant-cycle算法 而是蚁群系统 ant colony system算法

#include <fstream>
#include <iostream>
#include <vector>
#include <algorithm>
#include <cmath>
#include <ctime>
//没有实现带候选表的 和 带局部搜索的 策略
//仅仅是基本的ACS 08.8.5
//局部更新中的p*X  X是c 并且是个常数 待求..............................
//局部更新公式中 应该是p*tao0 也就是p*c 但是c的最优取值 是1.0/(n*Lnn)..........待定......
#include "ant.h"
#include "city.h"
using namespace std;

bool city::Init()
{
	ifstream infile("oliver30Tsp.dat",ios::in);
	if (!infile)
	{
		cerr<<"数据文件无法打开";
		return false;
	}
	double x,y,z;
	size_t i(0);
	while(infile>>x>>y>>z)
	{
		cityCoor[i][0]=y;
		cityCoor[i][1]=z;	
		++i;
	}
	for(i=0;i!=NUMOFCITY;++i)
		for(size_t j=i+1;j!=NUMOFCITY;++j)//i==j  时 无意义
		{
		//	tmpInfo[i][j]=0;//每次循环之后要复原
			posInfo[i][j]=c;//每次循环后保持旧值
			//preInfo保持不变为距离的倒数
			preInfo[i][j]=1.0/\
				sqrt(pow((cityCoor[i][0]-cityCoor[j][0]),2.0)+pow((cityCoor[i][1]-cityCoor[j][1]),2.0));
		//	tmpInfo[j][i]=tmpInfo[i][j];
			posInfo[j][i]=posInfo[i][j];
			preInfo[j][i]=preInfo[i][j];
		}
	for(i=0;i!=NUMOFANTS;++i)
	{
		antSet[i].GetTabu().reserve(NUMOFCITY);
		antSet[i].GetTabu().clear();
		antSet[i].GetTabu().push_back(i);//蚂蚁i放在城市i中
		shortestPath[i]=0;
	}

	//以下代码 求距离某城市的最近的距离 数组元素shortestSibling[i]就是距城市i最近的.
	double twoCityDis(0);
	for (i=1;i!=NUMOFCITY-1;++i)
	{
		twoCityDis=1.0/preInfo[i][i+1];
		for (size_t j(0);j!=NUMOFCITY;++j)
		{
			if (!(i==j))
			{
				if (1.0/preInfo[i][j]<twoCityDis)
				{
					twoCityDis=1.0/preInfo[i][j];
				}
			}
		}
		shortestSibling[i]=twoCityDis;
	}

	twoCityDis=1.0/preInfo[0][1];
	for (i=2;i!=NUMOFCITY;++i)
	{
		if (1.0/preInfo[0][i]<twoCityDis)
		{
			twoCityDis=1.0/preInfo[0][i];
		}
	}
	shortestSibling[0]=twoCityDis;

	twoCityDis=1.0/preInfo[0][NUMOFCITY-1];
	for (i=1;i!=NUMOFCITY-1;++i)
	{
		if (1.0/preInfo[NUMOFCITY-1][i]<twoCityDis)
		{
			twoCityDis=1.0/preInfo[NUMOFCITY-1][i];
		}
	}
	shortestSibling[NUMOFCITY-1]=twoCityDis;	

	return true;
}

void city::Run()
{//运行一次程序叫一次实验trial 一次程序有NCMAX次迭代 每次迭代NUMOFCITY步step
	srand(unsigned(time(0))); //产生随机种子  种子不一样 随机数才不一样
	//放在迭代的最外面是为了尽量使得每次迭代时的种子数不一样  从而加大随机数不同的概率
	//程序运行一次进行NCMAX次迭代 每次迭代有NUMODCITY-1步
	for(size_t i=0;i!=NCMAX;++i)
	{//进入一次迭代 
		RunCycle();
		ReSet();
	}
}

void city::RunCycle()
{//对于每一次迭代
	static size_t NC(0);
	cout<<++NC<<" Time"<<endl;
	
	for(size_t i=0;i!=NUMOFCITY-1;++i)//i表示本次迭代的步数 一共是NUMOFCITY-1个
	{//每一步
		//time(0)返回当前的时间
		for(size_t k(0);k!=NUMOFANTS;++k)
		{
			size_t currCity(*(antSet[k].GetTabu().end()-1));//蚂蚁k当前城市
			double q=rand()/(double)RAND_MAX; //产生随机数---------------------q
			int tempJ(-1);//蚂蚁k将要转移的城市
			if (q<=q0)
			{//在可选的里面去最大的作为转移目标
				double tempTrans(0.0),maxTrans(0.0);
				for(size_t j(0);j!=NUMOFCITY;++j)
				{
					if (antSet[k].GetTabu().end()==find(antSet[k].GetTabu().begin(),
					antSet[k].GetTabu().end(),j))
					{
						tempTrans=ProbTrans(posInfo[currCity][j],preInfo[currCity][j],beta);
						if(tempTrans>maxTrans)
						{
							maxTrans=tempTrans;
							tempJ=j;
						}
					}
				}

			}//A情况的tempJ求得 
			else
			{
				double sumAvalible=ProbTransSum(antSet[k].GetTabu(),currCity,beta);//p分母
				double pAccumlate=0.0;
				double pRand=rand()/(double)RAND_MAX; //产生随机数 
				for (size_t j=0;j!=NUMOFCITY;++j)
				{			
					if (antSet[k].GetTabu().end()==find(antSet[k].GetTabu().begin(),
						antSet[k].GetTabu().end(),j))
					{//轮盘赌法
						pAccumlate+=ProbTrans(posInfo[currCity][j],
							preInfo[currCity][j],beta)/sumAvalible;
						if (pRand<=pAccumlate)
						{
							tempJ=j; break;
							
						}
						
					}
				}
			}
		antSet[k].GetTabu().push_back(tempJ);
		//蚂蚁走完之后要进行局部信息素更新  当前步的当前蚂蚁
		//currCity就是公式中的r tempJ就是s
		posInfo[currCity][tempJ]=posInfo[currCity][tempJ]*(1-p)//可能有问题
			+p*1.0/(NUMOFCITY*shortestSibling[currCity]);//******************p*c？？？？
		posInfo[tempJ][currCity]=posInfo[currCity][tempJ];
		}
	}
}

void city::ReSet()
{	
	for(size_t i=0;i!=NUMOFANTS;++i)
	{//对于蚂蚁
		for(size_t j=0;j!=NUMOFCITY-1;++j)
		{//对于蚂蚁的禁忌表
			antSet[i].GetLength()+=
				1.0/preInfo[antSet[i].GetTabu().at(j)][antSet[i].GetTabu().at(j+1)];
			//cout<<antSet[i].GetTabu().at(j);
		}
		antSet[i].GetLength()+=
			1.0/preInfo[antSet[i].GetTabu().back()][antSet[i].GetTabu().front()];
	}

	//求出该次迭代的n个蚂蚁的最短
	double shortestPathEach=antSet[0].GetLength();
	size_t shortestAnt(0);
	for(i=1;i!=NUMOFANTS;++i)
	{
		if(antSet[i].GetLength()<shortestPathEach)
		{
			shortestPathEach=antSet[i].GetLength();
			shortestAnt=i;
		}
		//	cout<<antSet[i].GetLength()<<endl;
	}	
	cout<<"this time "<<shortestPathEach<<endl;
	
	//求出历史最短路径和路径图
	if (shortestPathEach<shortestPathFinal)
	{
		shortestPathFinal=shortestPathEach;
		//更新历史最短路径图
		for (size_t i(0);i!=NUMOFCITY;++i)
		{
			shortestPath[i]=antSet[shortestAnt].GetTabu().at(i);
		}
	}
	cout<<"SO far "<<shortestPathFinal<<endl;
	double tempInfo(0.0);
	//全局更新规则
	for(i=0;i!=NUMOFCITY;++i)
		for(size_t j=i+1;j!=NUMOFCITY;++j)
		{//对于每一条路径		
			if(AntInRoute(i,j,shortestPath,shortestPath+NUMOFCITY))
			{//判断时候ij边或ji边出现在历史最优路径中
				tempInfo=1.0/shortestPathFinal;
			}
			else
			{
				tempInfo=0.0;
			}
			//开始更新该ij边的信息素
			posInfo[i][j]=(1-alpha)*posInfo[i][j]+alpha*tempInfo;
			//		cout<<"更新后ij边"<<posInfo[make_pair(i,j)]<<" ";
			posInfo[j][i]=posInfo[i][j];
		}
		
		//清空禁忌表---放到循环开始 以便于在while判断语句中可以提取该信息
		for(i=0;i!=NUMOFANTS;++i)
		{//相关状态复原 
			antSet[i].GetTabu().clear();
			antSet[i].GetTabu().push_back(i);//蚂蚁归位
			antSet[i].GetLength()=0.0;//此时要进行下一次循环 所以路径必须清空
		}

}

void city::ShowShortestPath()
{
	cout<<"shortest path length: "<<shortestPathFinal<<endl;
	cout<<"path like: "<<ends;
	for (size_t i(0);i!=NUMOFCITY;++i)
	{
		cout<<shortestPath[i]<<'-';
	}
	cout<<endl;
}

inline double city::ProbTrans(const double & pos,
				 const double & pre,const double & beta)
{
	return (pos*pow(pre,beta));
}

bool city::AntInRoute(const size_t & i,const  size_t & j,
				const vector<size_t>::const_iterator & first,
				const vector<size_t>::const_iterator & last)
{//判断ij或ji是否连续的出现在由first last指向的序列中
	vector<size_t>::const_iterator iter(first);
	for(;iter!=last-1;++iter)
	{
		if( (*iter==i&&*(iter+1)==j)
			|| (*iter==j&&*(iter+1)==i) )
			break;
	}
	if(iter!=last-1)
		return true;
	return false;
}

double city::ProbTransSum(const vector <size_t > & tabu,const size_t & curCity,
					const  double & beta )
{
	double sum(0.0);
	for (size_t i(0);i!=NUMOFCITY;++i)
	{
		if (tabu.end()==find(tabu.begin(),tabu.end(),i))
		{
			sum+=ProbTrans(posInfo[curCity][i],
				preInfo[curCity][i],beta);
		}
	}
	return sum;
}