TSP_旅行商问题——模拟退火算法

#include<stdio.h>
#include<stdlib.h>
#include<string.h>
#include<time.h>
#include<math.h>
#include<iostream>
using namespace std;
#define T0 50000.0  // 初始温度
#define T_end (1e-8)
#define q  0.98   // 退火系数
#define L 1000  // 每个温度时的迭代次数，即链长
#define N 31  // 省会、首府城市数量
#define EARTH_RADIUS 6378.137
#define PI acos(-1)
int city_list[N]; // 用于存放一个解

//中国31个省会、首府城市的经纬度
double city_pos[N][2] =
    {
    {103.73, 36.03},{101.74, 36.56},{104.06, 30.67},{114.48, 38.03},
    {102.73, 25.04},{106.71, 26.57},{114.31, 30.52},{113.65, 34.76},
    {117.00, 36.65},{118.78, 32.04},{117.27, 31.86},{120.19, 30.26},
    {115.89, 28.68},{119.30, 26.08},{113.23, 23.16},{113.00, 28.21},
    {110.35, 20.02},{123.38, 41.80},{125.35, 43.88},{126.63, 45.75},
    {112.53, 37.87},{108.95, 34.27},{87.68, 43.77}, {116.46, 39.92},
    {121.48, 31.22},{106.54, 29.59},{117.20, 39.13},{111.65, 40.82},
    {108.33, 22.84},{91.11, 29.97}, {106.27, 38.47}
};
//省会城市名称
string city_name[N]={
	"兰州","西宁","成都","石家庄","昆明","贵阳","武汉","郑州","济南","南京",
	"合肥","杭州","南昌","福州","广州","长沙","海口","沈阳","长春","哈尔滨",
	"太原","西安","乌鲁木齐","北京","上海","重庆","天津","呼和浩特","南宁","拉萨",
	"银川"
};
//函数声明
double getDistance(double lat1, double lng1, double lat2, double lng2); 通过给定经纬度计算城市之间的距离，参考Google maps
double pathLen(int * arry); //给定城市顺序，计算城市环路径的长度
void init();//初始化城市序列
void swapCity();//交换两个城市顺序产生新序列

double rad(double d)
{
    return d * PI / 180.0;
}
 
double getDistance(double lat1, double lng1, double lat2, double lng2)
{
    double radLat1 = rad(lat1);
    double radLat2 = rad(lat2);
    double a = radLat1 - radLat2;
    double b = rad(lng1) - rad(lng2);
    double s = 2 * asin(sqrt(pow(sin(a/2),2) +	cos(radLat1)*cos(radLat2)*pow(sin(b/2),2)));
    s = s * EARTH_RADIUS;
    s = round(s * 10000) / 10000;
    return s;
}

double pathLen(int * arry){
	double path = 0; // 初始化路径长度
    int first_city = arry[0]; // 定位到第一个数字(城市序号)
    for(int i=0;i<N-1;i++)
    {
        int this_city = arry[i];
        int next_city = arry[i+1];
        double dis = getDistance(city_pos[this_city][0],city_pos[this_city][1],
        						 city_pos[next_city][0],city_pos[next_city][1]);
        path += dis;
    }
    int last_city = arry[N-1]; // 最后一个城市序号
    double last_dis = getDistance(city_pos[first_city][0],city_pos[first_city][1],
        						  city_pos[last_city][0],city_pos[last_city][1]);
    path = path + last_dis;
    return path; // 返回总的路径长度
}

void init(){
	for(int i=0;i<N;i++)
		city_list[i]=i;
}

void swapCity()
{
    double r1 = ((double)rand())/(RAND_MAX+1.0);
    double r2 = ((double)rand())/(RAND_MAX+1.0);
    int pos1 = (int)(N*r1); //第一个交叉点的位置
    int pos2 = (int)(N*r2);
    while(pos1==pos2){
    	r1 = ((double)rand())/(RAND_MAX+1.0);
    	r2 = ((double)rand())/(RAND_MAX+1.0);
    	pos1 = (int)(N*r1);
    	pos2 = (int)(N*r2);
    }
    int temp = city_list[pos1];
    city_list[pos1] = city_list[pos2];
    city_list[pos2] = temp;   // 交换两个点
}

int main(){

	srand((unsigned)time(NULL));//初始化随机数种子
	time_t start,finish;
	start = clock();//程序开始时间
	double T=T0;
	int count=0,countgap=0;//记录降温次数,差距次数
	init();
	cout<<"初始路径长度为:"<<pathLen(city_list)<<"km"<<endl;
	int city_list_old[N];//保存更新前的城市顺序
	double oldLen,newLen,gap;//oldLen为更新前城市环的路径距离，newLen为更新后城市环的路径距离，gap是两者差值
	double r;//0-1之间的随机数，用来决定是否接受新解
	while(T>T_end){
		bool change=false;
		for(int i=0;i<L;i++){
			memcpy(city_list_old,city_list,N*sizeof(int));
			swapCity();
			oldLen=pathLen(city_list_old);
			newLen=pathLen(city_list);
			gap=newLen-oldLen;
			if(gap<0)
				change=true;
			//Metropolis准则
			else{
				r=((double)rand())/(RAND_MAX+1.0);
				if(exp(-gap/T) <= r){
					memcpy(city_list,city_list_old,N*sizeof(int));//保留原来的城市路径
				}else{
					change=true;
				}
			}
		}
		count++;
		T*=q;
		if(!change)
			++countgap;
		else
			countgap=0;
		if(countgap==2)
			break;
		
	}
	finish=clock();

	double duration=((double)(finish-start))/CLOCKS_PER_SEC;
	printf("模拟退火算法，初始温度T0=%.2f,降温系数q=%.2f,每个温度迭代%d次,共降温%d次，得到的TSP最优路径为:\n",T0,q,L,count);
    for(int i=0;i<N-1;i++)  // 输出最优路径
    {
        cout<<city_name[city_list[i]]<<"-->";
    }
    cout<<city_name[city_list[N-1]]<<endl;
    double len = pathLen(city_list); // 最优路径长度
    cout<<"最优路径长度为:"<<len<<"km"<<endl;
    cout<<"程序运行耗时:"<<duration<<"sec"<<endl;
    return 0;
}