main.cpp

《智能优化算法》课的作业

#include <iostream>
#include <iomanip>
#include <vector>
#include <iterator>
#include <algorithm>
#include <fstream>
#include <math.h>
#include <time.h>

const int NCITY = 51;
struct CityCoor
{
	int sn, x, y;
};
struct Sequence
{
	double distance;
	std::vector<int> vSeqnc;
};

std::vector<CityCoor> vCityCoors(NCITY);
std::vector<std::vector<double>> vDistanceMatrix(NCITY);

double computDistance( const std::vector<int>& v );

namespace ga
{
	const int POPSIZE = NCITY*5;
	const int GENERATION = 5000;				//在变化显著的情况下为5000代
	const int SIMILAR_TIMES = 200;
	const double DIFF = 1e-15;					//量化变化不显著
	const double P_CROSSOVER = 0.3;				//交叉概率
	const double P_MUTATION = 0.1;				//变异概率
	int nGen = 0;

	std::vector<int> vCitySequence(NCITY);
	std::vector<Sequence> vChromsome(POPSIZE+1);	//第一个染色体只记录上一代遗传的最优解而不参与遗传运算
	std::vector<double> vProbability(POPSIZE+1);

	namespace itor
	{
		inline void printSequence(int sn) {std::cout << sn << " ";}
		inline void printCityCoors(const CityCoor& cc) {std::cout << cc.sn << " " << cc.x << " " << cc.y << "\n";}
		inline void printChromsome(const Sequence &s) {std::cout << s.distance << " ";}

		inline void itorDistance( Sequence& se ) {se.distance = computDistance(se.vSeqnc);}
		inline void initChromsome(Sequence& v);
		
		inline void testPrintSeqLen(const Sequence &s) {std::cout << s.vSeqnc.size() << " ";}
	}//itor::
	inline bool less_second(const Sequence &a, const Sequence &b) {return a.distance < b.distance;}
}//ga::

//////////////////////////////////////////////////////////////////////////
//迭代器
typedef std::vector<CityCoor>::iterator vitorCityCoor;
typedef std::vector<int>::iterator vitorSn;

typedef std::vector<int>::size_type vintSize;

typedef std::vector<CityCoor>::const_iterator cvitorCityCoor;
typedef std::vector<int>::const_iterator cvitorSn;
typedef std::vector<Sequence>::const_iterator cvitorChromsome;

//////////////////////////////////////////////////////////////////////////
std::vector<Sequence> vTemp(ga::POPSIZE+1);

//////////////////////////////////////////////////////////////////////////
//全局函数，具体作用可见Definition
bool readTSP_Data();
inline double myrand(double, double);
void crossover();
void mutation();
void selection();
void calcDistanceMatrix();

bool isValid(std::vector<int> v)
{
	sort(v.begin(), v.end());

	for (std::vector<int>::size_type i=0; i<v.size()-1; ++i)
		if (v[i] == v[i+1])
		{
			std::cout << "序列出错!!!!!!!!!\n";
			std::cout << "v[" << i << "]" << "=v[" << i+1 << "]\n";
			for_each(v.begin(), v.end(), ga::itor::printSequence);
			system("pause");
			return false;
		}
	std::cout << "序列正确\n";
	return true;
}

int main()
{
	//////////////////////////////////////////////////////////////////////////
	std::cout<<"#读入数据";
	if (!readTSP_Data()) 
	{
		std::cerr << "---读入城市坐标出错!";
		exit(1);
	}
	std::cout<<" OK\n";
	clock_t startTime = clock();

	std::cout<<"#计算城市距离矩阵";
	calcDistanceMatrix();
	std::cout<<" OK\n";

	//求适应函数
	ga::vProbability[0] = 0.05; 
	double a = 0.05;
	for (int i=1; i<=ga::POPSIZE; ++i)
	{
		a *= 0.95;
		ga::vProbability[i] = ga::vProbability[i-1] + a;
	}
	//////////////////////////////////////////////////////////////////////////
	std::cout<<"#初始化" << ga::POPSIZE << "个染色体：";
	for_each(ga::vChromsome.begin()+1, ga::vChromsome.end(), ga::itor::initChromsome);
	std::cout<<" OK\n";

	ga::vChromsome[0] = ga::vChromsome[1];
	//先淘汰一次，以优化初值
	sort(ga::vChromsome.begin()+1, ga::vChromsome.end(), ga::less_second);

	//////////////////////////////////////////////////////////////////////////
	std::cout<<"\n#开始遗传算法，结束条件为连续" << ga::SIMILAR_TIMES <<"代变化小于" << ga::DIFF <<"";

	int generWidth = static_cast<int>(ceil(log10(static_cast<double>(ga::GENERATION))));
	for (int i=1; i<=ga::GENERATION; ++i)
	{
		selection();
		double prevDistance = ga::vChromsome[0].distance;

		srand(static_cast<unsigned>(time(NULL)));
		crossover();
		mutation();

		//对每个染色体求距离
		for_each(ga::vChromsome.begin()+1, ga::vChromsome.end(), ga::itor::itorDistance);
		//进化，好解(小)在前，上一代的最优解也带入
		sort(ga::vChromsome.begin(), ga::vChromsome.end(), ga::less_second);

		std::cout << "\n    第" ;
		std::cout << std::setfill(' ') << std::setw(generWidth) << i;
		std::cout << "代       最短距离="<< ga::vChromsome[0].distance;

		if (prevDistance - ga::vChromsome[0].distance < ga::DIFF) ++ga::nGen;
		else ga::nGen = 0;

		if (ga::SIMILAR_TIMES < ga::nGen)	break;
	}

	std::cout<<"\n\n#检验结果->";
	isValid(ga::vChromsome[0].vSeqnc);

	for_each(ga::vChromsome[0].vSeqnc.begin(), ga::vChromsome[0].vSeqnc.end(), ga::itor::printSequence);
	std::cout<<"\n";

	std::cout<<"\n共使用时间: " << static_cast<double>(clock()-startTime) / CLOCKS_PER_SEC<< "秒\n";
	system("pause");
	return 0;
}

bool readTSP_Data()
{
	try
	{
		std::ifstream fin("TSP.data");
		for (vitorCityCoor i=vCityCoors.begin(); i != vCityCoors.end(); ++i)
		{
			fin >> (*i).sn >> (*i).x >> (*i).y;
		}
		fin.close();
	}
	catch (...)
	{
		return false;
	}
	return true;
}

//贪婪算法Right
int GreedyRight( std::vector<int> &cityrouter, int nCity )
{
	bool bFindCity = false;
	vitorSn iter_city;
	for(iter_city=cityrouter.begin();iter_city!=cityrouter.end(); ++iter_city)
		if( *iter_city == nCity )
		{
			bFindCity = true;
			break;
		}
	
	if( bFindCity )
	{
		++iter_city;
		if( iter_city == cityrouter.end() )
			iter_city = cityrouter.begin();

		return *iter_city;
	}
	else
		return -1;
}

//贪婪算法Left
int GreedyLeft( std::vector<int> &cityrouter, int nCity )
{
	bool bFindCity = false;
	vitorSn iter_city;
	for( iter_city=cityrouter.begin();iter_city!=cityrouter.end(); ++iter_city )
		if( *iter_city == nCity )
		{
			bFindCity = true;
			break;
		}
	
	if( bFindCity )
	{
		if( iter_city == cityrouter.begin() )
			return cityrouter.back();
		else
		{
			--iter_city;
			return *iter_city;
		}
	}
	else
		return -1;
}

void GreedyErase( std::vector<int> &cityrouter, int ndelcity )
{
	bool bFindCity = false;
	vitorSn iter_city;
	for( iter_city=cityrouter.begin();iter_city!=cityrouter.end(); ++iter_city )
		if( *iter_city == ndelcity )
		{
			bFindCity = true;
			break;
		}
	
	if( bFindCity ) cityrouter.erase( iter_city );
}

//************************************
// Method:    doCrossover
// FullName:  doCrossover
// Access:    public 
// Returns:   void
// Parameter: int nFatherA
// Parameter: int nFatherB
// Qualifier: 贪心交叉方式（Greedy Crossover），
// 具体算法可参见 谢胜利,等.求解TSP问题的一种改进的遗传算法[J].计算机工程与应用,2002(8):58～245。见下载目录
//************************************
void doCrossover(int nFatherA, int nFatherB)
{
	int randomcity, nowopcity, nextopcity, rightA, rightB, leftA, leftB;
	std::vector<int>  SonA, SonB;

	randomcity = static_cast<int>(myrand(1, NCITY));

	nowopcity = randomcity;
	SonA.push_back( nowopcity );

	std::vector<int> FatherA = ga::vChromsome[nFatherA].vSeqnc;
	std::vector<int> FatherB = ga::vChromsome[nFatherB].vSeqnc;

	while( FatherA.size() > 1 && FatherB.size() > 1 )
	{
		rightA = GreedyRight( FatherA, nowopcity );
		rightB = GreedyRight( FatherB, nowopcity );

		if( vDistanceMatrix[nowopcity-1][rightA-1]  < vDistanceMatrix[nowopcity-1][rightB-1] )
		{
			SonA.push_back( rightA );
			nextopcity = rightA;
		}
		else
		{
			SonA.push_back( rightB );
			nextopcity = rightB;
		}

		GreedyErase( FatherA, nowopcity );
		GreedyErase( FatherB, nowopcity );
		nowopcity = nextopcity;
	}

	nowopcity = randomcity;
	SonB.push_back( nowopcity );
	FatherA = ga::vChromsome[nFatherA].vSeqnc;
	FatherB = ga::vChromsome[nFatherB].vSeqnc;
	
	while( FatherA.size() > 1 && FatherB.size() > 1 )
	{
		leftA = GreedyLeft( FatherA, nowopcity );
		leftB = GreedyLeft( FatherB, nowopcity );

		if( vDistanceMatrix[nowopcity-1][leftA-1] < vDistanceMatrix[nowopcity-1][leftB-1])
		{
			SonB.push_back( leftA );
			nextopcity = leftA;
		}
		else
		{
			SonB.push_back( leftB );
			nextopcity = leftB;
		}

		GreedyErase( FatherA, nowopcity );
		GreedyErase( FatherB, nowopcity );
		nowopcity = nextopcity;
	}

	swap(ga::vChromsome[nFatherA].vSeqnc, SonA);
	swap(ga::vChromsome[nFatherB].vSeqnc, SonB);
}

//************************************
// Method:    crossover
// FullName:  crossover
// Access:    public 
// Returns:   void
// Qualifier: 交配
//************************************
void crossover()			
{
	std::vector<int> vecCrossoverIndexs;
	double random;

	for( int i=1;i<=ga::POPSIZE;i++ )
	{
		random = static_cast<double>(myrand(0,1));
		if( random < ga::P_CROSSOVER )
			vecCrossoverIndexs.push_back( i );
	}

	size_t CrossoverNumber = vecCrossoverIndexs.size();
	if( CrossoverNumber%2 != 0 )
		vecCrossoverIndexs.pop_back();

	CrossoverNumber = vecCrossoverIndexs.size();
	for(size_t i=0; i<CrossoverNumber; i+=2)
	{
		int nFatherA = vecCrossoverIndexs[i];
		int nFatherB = vecCrossoverIndexs[i+1];

		doCrossover( nFatherA, nFatherB);
	}
}

void ga::itor::initChromsome(Sequence& v)
{
	v.vSeqnc.resize(NCITY);

	for (int i=1; i <= NCITY; ++i)
		v.vSeqnc[i-1] = i;

	std::random_shuffle(v.vSeqnc.begin(), v.vSeqnc.end());
	v.distance = computDistance(v.vSeqnc);
	std::cout << ".";
}

//************************************
// Method:    myrand
// FullName:  myrand
// Access:    public 
// Returns:   double
// Qualifier:
// Parameter: double a
// Parameter: double b
//  Uniform Distribution
//  return a random num in area [a,b]
//************************************
double myrand(double a, double b) 
{
	double y;
	if(a>b) {
		printf("\nThe first parameter should be less than the second!");
		exit(1);
	}
	//////////////////////////////////////////////////////////////////////////
	//rand() can reture a number in [0, RAND_MAX]
	y = static_cast<double>(rand())/RAND_MAX;

	return (a+(b-a)*y);
}

//************************************
// Method:    mutation
// FullName:  mutation
// Access:    public 
// Returns:   void
// Qualifier: 变异，对两个随机位置之间的城市随机重排
//************************************
void mutation()
{
	vitorSn it;
	for (int i=1; i <= ga::POPSIZE; ++i)
		if (ga::P_MUTATION > myrand(0,1))
		{
			int left = static_cast<int>(myrand(0, NCITY/2));
			int right = static_cast<int>(myrand(NCITY/2, NCITY));
			it = ga::vChromsome[i].vSeqnc.begin();

			std::random_shuffle(it+left, it+right);
		}
}

double computDistance( const std::vector<int>& v )
{
	double tmp = 0.0;
	
	for (int it=0; it<NCITY-1; ++it) 
		tmp += vDistanceMatrix[v[it]-1][v[it+1]-1];

	tmp += vDistanceMatrix[v[NCITY-1]-1][v[0]-1];
	return tmp;
}

//************************************
// Method:    selection
// FullName:  selection
// Access:    public 
// Returns:   void
// Qualifier: 选择，轮盘赌，让前面的解（好解）以较大概率被选中
//************************************
void selection()
{
	double r;
	int label;

	vTemp[0] = ga::vChromsome[0];	//第一个仅仅作记录用
	for (int i=1; i <= ga::POPSIZE; ++i)
	{
		r = myrand(0, ga::vProbability[ga::POPSIZE]);
		label = 0;
		for (int j=0; j <= ga::POPSIZE; ++j)
		{
			if (r <= ga::vProbability[j])
			{
				label = j;
				break;
			}
		}
		vTemp[i] = ga::vChromsome[label];
	}

	swap(ga::vChromsome, vTemp);
}

//************************************
// Method:    calcDistanceMatrix
// FullName:  calcDistanceMatrix
// Access:    public 
// Returns:   void
// Qualifier: 预先计算出城市间的距离矩阵供后面查询，典型的空间换时间
//************************************
void calcDistanceMatrix()
{
	double dltX, dltY;

	std::vector<std::vector<double>>::iterator it = vDistanceMatrix.begin();
	for(; it != vDistanceMatrix.end(); ++it)
		(*it).resize(NCITY);

	for (int i=0; i<NCITY; ++i)
		for (int j=i; j<NCITY; ++j)
			if (i == j) vDistanceMatrix[i][j] = 0.0;
			else
			{
				dltX = vCityCoors[i].x - vCityCoors[j].x;
				dltY = vCityCoors[i].y - vCityCoors[j].y;
				vDistanceMatrix[i][j] = sqrt(static_cast<double>(dltX*dltX + dltY*dltY));
				vDistanceMatrix[j][i] = vDistanceMatrix[i][j];
			}
}