Simulated Annealing(模拟退火算法)

/*
Simulated Annealing(模拟退火算法)
求解旅行商问题(TSP)
网上给的数据是31个省会的坐标，蚁群算法得到的结果是：15378
我算的结果中，最好的一次是：15495
*/
#include<iostream>
#include<cstdio>
#include<cstdlib>
//#include<ctime>
#include<cmath>
#include<time.h>
#define N 31       //城市个数
#define Tmax 8000  //初始温度
#define Tmin 1E-10 //终止温度
#define RATE 0.95  //温度衰减率
#define in_loop 13000 //内层循环次数
#define out_loop 2000 //外层循环次数
#define p_limit 10000 //概率选择次数

using namespace std;

//31个省会x和y的坐标
double x[N]={1304,3639,4177,3712,3488,3326,3238,4196,4312,4386,3007,2562,2788,2381,1332,3715,3918,4061,3780,3676,4029,4263,3429,3507,3394,3439,2935,3140,2545,2778,2370};
double y[N]={2312,1315,2244,1399,1535,1556,1229,1004,790,570,1970,1756,1491,1676,695,1678,2179,2370,2212,2578,2838,2931,1908,2367,2643,3201,3240,3550,2357,2826,2975};

//两个城市之间的距离
double d[N][N];

void init(){
    //初始化任意两个城市的距离
    for(int i=0;i<N;i++)for(int j=0;j<i;j++){
        d[i][j]=d[j][i]=sqrt((y[j]-y[i])*(y[j]-y[i])+(x[j]-x[i])*(x[j]-x[i]));
    }
    for(int i=0;i<N;i++)d[i][i]=0;

    srand((unsigned)time(0));
}

//路径
class path{
public:
    path(){}
    ~path(){}
    int city[N];    //依次经过的城市编号
    double dis;     //总的距离
    void totalLen(){
        dis=0;
        for(int i=0;i<N-1;i++)dis+=d[city[i]][city[i+1]];
        dis+=d[city[0]][city[N-1]];
    }
};

//最优解
path bestpath;

//产生新路径
path newpath(path prepath){
    path newp = prepath;
    int x,y,t;
    do{
        x=rand()%N;
        y=rand()%N;
    }while(x==y);
    t=newp.city[x];newp.city[x]=newp.city[y];newp.city[y]=t;
    newp.totalLen();
    return newp;
}

//Annealing
void Anne()
{
    init();
    //随机选一个初始解
    for(int i=0;i<N;i++)bestpath.city[i]=i;
    bestpath.totalLen();
    int out_t=0,p_t=0;
    double T=Tmax,delta,prob, rnd;
    path np,cp; //新路径，当前路径
    cp=bestpath;
    while(out_t<out_loop&&T>Tmin){
        for(int i=0;i<in_loop;i++){
            np=newpath(cp);
            if(np.dis<cp.dis){
                cp=np;
                out_t=0;
                p_t=0;
            }
            else{
                delta=np.dis-cp.dis;
                prob=exp(-delta/T);
                rnd=rand()%10000/10000.0;
                if(prob>rnd)cp=np;
                p_t++;
            }
            if(p_t>p_limit){
                out_t++;
                break;
            }
        }
        if(cp.dis<bestpath.dis)bestpath=cp;
        T*=RATE;
        printf("dis= %f\n",bestpath.dis);
    }
}

int main()
{
    Anne();
    return 0;
}