人工智能作业homework3——模拟退火解决TSP

受固体退火过程的启发，Kirkpatrick等人意识到组合优化问题与固体退火过程的类似性，将组合优化问题类比为固体的退火过程，提出了求解组合优化问题的模拟退火算法。

表7.3给出了组合优化问题与固体退火过程的类比关系。

表7.3：组合优化问题与退火过程的类比

在求解组合优化问题时，首先给定一个比较大的t值，这相当于给定一个比较高的温度T。随机给定一个问题的解i，作为问题的初始解。在给定的t下，随机产生一个问题的解j，j∈N(i)，其中N(i)是i的邻域。从解i到新解j的转移概率，按照Metropolis准则确定，即：

如果新解j被接受，则以解j代替解i，否则继续保持解i。重复该过程，直到在该控制参数t下达到平衡。与退火过程中的温度T缓慢下降相对应，在进行足够多的状态转移之后，控制参数t要缓慢下降，并在每个参数t下，重复以上过程，直到控制参数t降低到足够小为止。最终我们得到的是该组合优化问题的一个最优解。由于这样一个过程模拟的是退火过程，所以被称为模拟退火算法。

该算法有内外两层循环。内循环模拟的是在给定温度下系统达到热平衡的过程。每次循环随机的产生一个新解，然后按照Metropolis准则，随机的接受该解。算法中的Random(0, 1)，是一个在[0, 1]间均匀分布的随机数发生器，与从解i到劣解j的转移概率相结合，模拟系统是否接受了劣解j。外循环模拟的是温度的下降过程，控制参数t_k起到与温度T相类似的作用，表示的是第k次循环时系统所处的温度。算法中的Drop(t_k)是一个温度下降函数，它按照一定的原则实施温度的缓慢下降。

模拟退火算法与局部搜索算法很相似，二者最大的不同是模拟退火算法按照Metropolis准则随机的接受一些劣解，即指标函数值大的解。当温度比较高时，接受劣解的概率比较大，在初始高温下，几乎以接近100％的概率接受劣解。随着温度的下降，接受劣解的概率逐渐减少，直到当温度趋于0时，接受劣解的概率也同时趋于0。这样，将有利于算法从局部最优解中跳出，求得问题的全局最优解。

上述模拟退火算法只是给出了一个算法的框架，其中重要的三个条件：初始温度的选取，内循环的结束条件和外循环的结束条件，算法中都没有提及，而这正是模拟退火算法的关键所在。因为正像前面叙述过的那样，对于固体退火过程来说，要最终使得物理系统以概率1处于能量最小的一个状态，在退火过程中必须满足以下三个条件：

（1）初始温度必须足够高；

（2）在每个温度下，状态的交换必须足够充分；

（3）温度T的下降必须足够缓慢。

这三个条件刚好是与算法中未提及的三个重要条件相对应的。与固体退火过程一样，为了使得模拟退火算法以概率1求解到问题的最优解，则至少也要满足这三个条件。然而“初始温度必须足够高，状态交换必须足够充分，温度的下降必须足够缓慢”这样的条件是与我们试图给出求解组合优化问题的低复杂度算法的初衷相违背的。如果模拟退火算法仍然是一个指数复杂度的算法，则对于求解复杂组合优化问题不会带来任何意义下的帮助。现在的问题是，如何弱化一些条件，使得我们能够在一个多项式时间复杂度内，求得一个组合优化问题的满意解。在下一节我们将讨论这些问题，给出一些如何确定初始温度以及内、外循环结束条件的基本方法。

参数的确定

从以上的分析我们知道，模拟退火算法以概率1找到全局最优解的基本条件，是初始温度必须足够高，在每个温度下状态的交换必须足够充分，温度t的下降必须足够缓慢。因此初始的温度t₀，在每个温度下状态的交换次数、温度t的下降方法，以及温度下降到什么程度算法结束等参数确定，成为模拟退火算法求解问题时必须要考虑的问题。因为从理论上来说，模拟退火算法逐渐达到最优解的能力是以搜索过程的无限次状态转移为前提的，作为一种最优化的求解算法，其算法的时间复杂性仍然是指数时间的，无法用于大规模的组合优化问题求解。但是对于很多实际问题，正如我们在第一节已经讨论的那样，求解问题最优解的意义并不大，一个满意解就足够了。而是否能在一个多项式的时间内得到问题的满意解，则是我们关心的主要问题。

并不是任何一组参数，都能够保证模拟退火算法收敛于某个近似解，大量的实验表明，解的质量与算法的运行时间是成正比关系的，很难做到两全其美。下面我们给出模拟退火算法中一些参数或者准则的确定方法，试图在求解时间与解的质量之间作一个折中的选择。这些参数或者准则包括：

（1）初始温度t₀；

（2）温度t的衰减函数，即温度的下降方法；

（3）算法的终止准则，用终止温度t_f或者终止条件给出；

（4）每个温度t下的马尔可夫链长度L_k。

仿照固体的升温过程，也可以通过逐步升温的方法，得到一个合适的初始温度。方法如下：

（1）给定一个希望的初始接受概率P₀，给定一个较低的初始温度t₀，比如t₀＝1；

（2）随机的产生一个状态序列，并计算该序列的接收率：

如果接收率大于给定的初始接受概率P₀，则转（4）；

（3）提高温度，更新t₀，转（2）；

（4）结束。

其中更新t₀可以采用每次加倍的方法：

t₀=2×t₀

也可以采用每次加固定值的方法：

t₀=t₀＋T

这里的T为一个事先给定的常量。

2．温度的下降方法

退火过程要求温度下降足够缓慢，常用的温度下降方法有以下三种：
（1）等比例下降。
该方法通过设置一个衰减系数，使得温度每次以相同的比率下降：
 ，k=0,1,....。                                （7.36）
其中tk是当前温度，tk+1是下一个时刻的温度， 是一个常数。 越接近于1，温度下降的越慢，一般可以选取0.8～0.95左右的一个值。该方法简单实用，是一种常用的温度下降方法。
（2）等值下降
该方法每次温度的下降幅度是一个固定值：
                                         （7.37）
设K是希望的温度下降总次数，则：
                                               （7.38）
其中t0是初始温度。
该方法的好处是可以控制总的温度下降次数，但由于每次温度下降的是一个固定值，如果设置过小，在高温时温度下降太慢，如果设置的大，在低温下温度下降的又过快。

3．每一温度下的停止准则
在每一个温度下，模拟退火算法都要求产生足够的状态交换，如果用Lk表示在温度tk下的迭代次数的话，则Lk应使得在这一温度下的马尔可夫链基本达到平稳状态。
有以下几种常用的停止准则：
（1）固定长度方法
这是最简单的一种方法，在每一个温度下，都使用相同的Lk。Lk的选取与具体的问题相关，一般与邻域的大小直接关联，通常选择为问题规模n的一个多项式函数。如在n城市的旅行商问题中，如果采用交换两个城市的方法产生邻域的话，邻域的大小为 ，则Lk可以选取如Cn、Cn2等，其中C为常数。
（2）基于接受率的停止准则
由前面我们对退火过程的分析知道，在比较高的温度时，系统处于每一个状态的概率基本相同，而且每一个状态的接受概率也接近于1。因此在高温时，即便比较小的迭代数，也可以基本达到平稳状态。而随着温度的下降，被拒绝的状态数随之增加，因此在低温下迭代数应增加，以免由于迭代数太少，而过早的陷入局部最优状态。因此一个直观的想法就是随着温度的下降适当的增加迭代次数。
一种方法就是，规定一个接受次数R，在某一温度下，只有被接受的状态数达到R时，在该温度下的迭代才停止，转入下一个温度。由于随着温度的下降，状态被接受的概率随之下降，因此这样的一种准则是满足随着温度的下降适当的增加迭代次数的。但由于在温度比较低时，接受概率很低，为了防止出现过多的迭代次数，一般设置一个迭代次数的上限，当迭代次数达到上限时，即便不满足接受次数R，也停止这一温度的迭代过程。
与上一种方法相类似的，可以规定一个状态接受率R，R等于该温度下接受的状态数除以生成的总状态数。如果接受率达到了R，则停止该温度下的迭代，转入下一个温度。为了防止迭代次数过少或者过多，一般定义一个迭代次数的下限和上限，只有当迭代次数达到了下限并且满足所要求的接受率R时，或者达到了迭代次数的上限时，才停止这一温度的迭代。
还可以通过引入“代”的概念来定义停止准则。在迭代的过程中，若干相邻的状态称为“一代”，如果相邻两代的解的指标函数差值小于规定的值的话，则停止该温度下的迭代。
在某一温度下的迭代次数与温度的下降是紧密相关的。如果温度下降的幅度比较小，则相邻两个温度之间的平稳分布相差也应该比较小。一些研究表明，过长的迭代次数对提高解的质量关系不大，只会导致增加系统的运算时间。因此一般选取比较小的温度衰减值，只要迭代次数适当大就可以了。

4．算法的终止原则

模拟退火算法从初始温度t₀开始逐步下降温度，最终当温度下降到一定值时，算法结束。合理的结束条件，应使得算法收敛于问题的某一个近似解，同时应能保证解具有一定的质量，并且应在一个可以接受的有限时间内算法停止。

一般有下面几种确定算法终止的方法。

（1）零度法

从理论上讲，当温度趋近于0时，模拟退火算法才结束。因此，可以设定一个正常数，当时，算法结束。

（2）循环总控制法

给定一个指定的温度下降次数K，当温度的迭代次数达到K次时，则算法停止。这要求给定一个合适的K。如果K值选择不合适，对于小规模问题将导致增加算法无谓的运行时间，而对于大规模问题，则可能难于得到高质量的解。

（3）无变化控制法

随着温度的下降，虽然由于模拟退火算法会随机的接受一些不好的解，但从总体上来说，得到的解的质量应该逐步提高，在温度比较低时，更是如此。如果在相邻的n个温度中，得到的解的指标函数值无任何变化，则说明算法已经收敛。即便是收敛于局部最优解，由于在低温下跳出局部最优解的可能性很小，因此算法可以终止。

以上内容来源

马少平，朱小燕，人工智能，清华大学出版社

作业相关：

代码是参考上述伪代码写的，仍建议按着命令行方式调试，如下：

SA.h

#ifndef SA_H
#define SA_H


#include <iostream>
#include <fstream>
#include <vector>
#include<stdlib.h>
#include <time.h>

using namespace std;


const double INIT_T=300; //初始温度
const double RATE=0.92; //温度衰减率
const int IN_LOOP=2000; //内循环次数
const double FINAL_T=0.001;  //外层循环温度终止条件
const int LIMIT=1000;  //概率选择上限


typedef struct
{
	char name;
	double x;
	double y;
}_city;



class SA{
private:
	int CityNum;  //读取城市个数
	_city city;
	vector<_city> path; //每个温度下平衡状态
	double pointDist[20][20];
	char bestPath[20];
public:
	void input(const string &);
	double calPointDist(char,char);
	double calLength(vector<_city> );
	vector<_city> getnext(vector<_city>);
	void solve(const string &); //输出路径
	SA(){
		srand((int)time(0));//初始化随机数
	}
};


#endif

#include "SA.h"
#include<string>



void SA::input(const string &fileName){
	ifstream fin;
	fin.open(fileName.c_str());
	if(!fin){  //判断文件是否正常打开
		cout<<"Unable to open the file!"<<endl;
		exit(1);
	}
	string line;
	int i=0,j=0;
	getline(fin,line);
	CityNum=atoi(line.c_str());     //读取城市个数
	while(getline(fin,line)){       //逐行读取输入流中的文件，直至该输入流结束
		city.name=line[0];
		string::size_type pos1, pos2;
		pos1 = line.find_first_of('\t');//找到line第一个'\t'
		pos2 = line.find_last_of('\t'); //找到line最后一个'\t'
		city.x = atof(line.substr(pos1+1,6).c_str());
		city.y = atof(line.substr(pos2+1,6).c_str());
		path.push_back(city);
	}	
	
	double l;
	for(int i=0;i<CityNum;i++)  //用二维数组存储两城市间距离
		for(int j=i+1;j<CityNum;j++){
			l=(path[j].x-path[i].x)*(path[j].x-path[i].x);
			l+=(path[j].y-path[i].y)*(path[j].y-path[i].y);
			pointDist[i][j]=sqrt(l);
			pointDist[j][i]=sqrt(l);
		}


}

double SA::calPointDist(char c1,char  c2){  //返回两城市间距离
	return pointDist[c1-65][c2-65];
}

vector<_city> SA::getnext(vector<_city> curPath){
	 int i,j;
	 do 
	 {
		 i=(int)(CityNum*rand()/(RAND_MAX+1.0));  //生成[0,CityNum-1];
		 j=(int)(CityNum*rand()/(RAND_MAX+1.0));
	 } while (!(i<j));       //生成两个不同随机数

	 while (i<j)
	 {
		 swap(curPath[i++], curPath[j--]);  //逆序法生成新解
	 }
	 return curPath;
}

double SA::calLength(vector<_city> curPath){  //计算当前路径距离
	double totalLength=0;
	for(int i=0;i<CityNum-1;i++)
		totalLength+=calPointDist(curPath[i].name,curPath[i+1].name);
	totalLength+=calPointDist(curPath[CityNum-1].name,curPath[0].name);
	return totalLength;
}


void SA::solve(const string &fileName){
	vector<_city> curPath;
	vector<_city> newPath;
	double curLen,newLen;
	curPath=path;curLen=calLength(curPath);
	double delta;
	double p;
	double T=INIT_T;
	while(true){
		for(int i=0;i<IN_LOOP;i++){  //固定步数结束,即达到温度平衡条件	
			newPath=getnext(curPath);newLen=calLength(newPath);
			delta=newLen-curLen;  //判断指标函数增减
			if(delta<0){  //更新长度
				curPath=newPath;curLen=calLength(curPath);


			}else{
				p = (double)(1.0*rand()/(RAND_MAX+1.0));  //以一定概率接受差解
				double ee=exp(delta/T);
				if(exp(-delta/T) > p)  
				{  
					curPath = newPath;
					curLen=calLength(curPath);  
				} 

			}

		}

	
		path=curPath;

		for(int i=0;i<CityNum;i++){
			cout<<curPath[i].name<<" ";
		}
		cout<<curLen<<endl;

		static ofstream fout(fileName.c_str());
		if(!fout){
			cout<<"error!"<<endl;
			exit(1);
		}
		for(int i=0;i<CityNum;i++)
			fout<<path[i].name<<" ";  //逐个写入输出文件
		fout<<calLength(path)<<endl;

		if(T<FINAL_T)
			break;
		T=T*RATE;
	}


}

#include "SA.h"

void main(int argc,char *argv[]){
	SA tsp;
	//tsp.input("G:\\homework3\\TSP20.txt");
	tsp.input(argv[1]);
	//tsp.solve("G:\\homework3\\O.txt");
	tsp.solve(argv[2]);
    
	//getchar();

}