模拟退火解TSP问题

我根据这篇博客，实现了模拟退火用来解决TSP问题。

/*
// #SA,Simulated Annealing# 模拟退火算法
        while( T > T_min )
        {
            dE = J( Y(i+1) ) - J( Y(i) ) ;
            if ( dE >=0 ) // 表达移动后得到更优解，则总是接受移动
                Y(i+1) = Y(i) ; // 接受从Y(i)到Y(i+1)的移动
            else
            {
            // 函数exp( dE/T )的取值范围是(0,1) ，dE/T越大，则exp( dE/T )也越大
                if ( exp( dE/T ) > random( 0 , 1 ) )
                Y(i+1) = Y(i) ; //接受从Y(i)到Y(i+1)的移动
            }
            T = r * T ; //降温退火 ，0<r<1 。r越大，降温越慢；r越小，降温越快
            // 若r过大，则搜索到全局最优解的可能会较高，但搜索的过程也就较长。
            // 若r过小，则搜索的过程会很快，但最终可能会达到一个局部最优值
        }
        #产生新遍历路径的方法:
        #随机选择2个节点
        #将路径中这2个节点间的节点顺序逆转。
*/
#include <cstdio>
#include <cstdlib>
#include <ctime>
#include <iostream>
#include <string>
#include <cstring>
#include <fstream>
#include <cmath>
#include <iomanip>
using namespace std;

#define MAXN 64
clock_t start, finish;
double x[MAXN];  // x[i] : coordinate x of city i
double y[MAXN];  // y[i] : coordinate y of city i
double dist[MAXN][MAXN];  // dist[i][j] is the distance between city i&j
double calDist(double x1, double x2, double y1, double y2);  // distance between (x1, y1)&(x2, y2)
void printRunningTime(ofstream& foutput, const int& n);
int path[MAXN];  // record the current path
int tmpPath[MAXN];
double Result;  // the current distance of path
const string fName = "berlin52.tsp";  // input file name
double T, T_min = 1E-10;   // tempareture of system and minimal of it
double r;  // annealing rate
void printPathConsole(const int& n);
void printPathFile(ofstream& foutput, const int&n);

void copyPath(int*, int*, int);  // copy from s to t of length n.
void inversePath(int*, int, int, int);
double getPathDis(int*, int n);  // calculate the distance of current path

int main()
{
    ifstream finput;
    ofstream foutput("result.txt");
    finput.open(fName.c_str());
    int n;  // #cities
    finput >> n;
    r = 0.9999;
    for(int i = 0; i < n; i++)
        finput >> x[i] >> y[i];//,cout << x[i] << ' ' << y[i] << endl;
    for(int i = 0; i < n-1; i++){
        for(int j = i+1; j < n; j++){
            double tmpDist = calDist(x[i], x[j], y[i], y[j]);
            dist[i][j] = tmpDist;
            dist[j][i] = tmpDist;
        }
    }
    srand((unsigned)time(0)), Result = 0;
    for(int i = 0; i < n; i++) path[i] = i; path[n] = 0;
    for(int i = 0; i < n-1; i++) Result += dist[i][i+1]; Result += dist[n-1][0];
    T = 30000000;  // initial temperature

    start = clock();   // record start time
    while(T > T_min)
    {
        int p1, p2;
        while((p1 = rand()%n) == (p2 = rand()%n)); // random get 2 different positions
        if(p1 > p2) swap(p1, p2);
        // find the new path after swap points between 2 position
        copyPath(path, tmpPath, n);  // store current path in tmpPath
        inversePath(tmpPath, n, p1, p2);

        double tmpResult = getPathDis(tmpPath, n);
        // get the new Result
        double dE = (Result - tmpResult);
        if ( dE >= 0 ) // 表达移动后得到更优解，则总是接受移动
        {
            inversePath(path, n, p1, p2);
            Result = tmpResult;
        }
        else
        {
            // 函数exp( dE/T )的取值范围是(0,1) ，dE/T越大，则exp( dE/T )也越大
            double condition = (double)rand()/RAND_MAX;
            if ( exp(dE/T) <= 1 && exp(dE/T) > condition)
            {    //接受从Y(i)到Y(i+1)的移动
                inversePath(path, n, p1, p2);
                Result = tmpResult;
            }
        }
        T = r * T ; //降温退火 ，0<r<1 。r越大，降温越慢；r越小，降温越快
        /*
        若r过大，则搜索到全局最优解的可能会较高，但搜索的过程也就较长。
        若r过小，则搜索的过程会很快，但最终可能会达到一个局部最优值   */
    }
    finish = clock(); // record finish time
    cout << "Instance is: " << fName << endl;
    foutput << "Instance is: " << fName << endl;
    cout << "Result is: " << Result << endl;
    foutput << "Result is: " << Result << endl;
    printPathConsole(n);
    printPathFile(foutput, n);
    printRunningTime(foutput, n);
    finput.close();
    foutput.close();

   return 0;
}

double calDist(double x1, double x2, double y1, double y2){
    return sqrt( (x1-x2)*(x1-x2) + (y1-y2)*(y1-y2) );
}
void printRunningTime(ofstream& foutput, const int& n){
    double duration;
    duration = (double)(finish - start) / CLOCKS_PER_SEC;
    cout << fixed << setprecision(6) << duration << " seconds" << endl;
    foutput << fixed << setprecision(6) << duration << " seconds" << endl;
    return;
}
void printPathConsole(const int& n){
    cout << "The route is: ";
    for(int i = 0; i <= n; i++)
        cout << path[i] << ' ';
    cout << endl;
    return;
}
void printPathFile(ofstream& foutput, const int& n){
    foutput << "The route is: ";
    for(int i = 0; i <= n; i++)
        foutput << path[i] << ' ';
    foutput << endl;
    return;
}
void copyPath(int* s, int* d, int n){
    for(int i = 0; i <= n; i++)
        d[i] = s[i];
    return;
}
void inversePath(int* a, int n, int tmp1, int tmp2){
    while(tmp1 < tmp2)
        swap(a[tmp1++], a[tmp2--]);
    a[n] = a[0];
    return;
}
double getPathDis(int* a, int n){
    double dis = 0;
    for(int i = 0; i < n; i++)
        dis += dist[a[i]][a[i+1]];
    return dis;
}

这里是我用到的数据集：

52
565.0 575.0
25.0 185.0
345.0 750.0
945.0 685.0
845.0 655.0
880.0 660.0
25.0 230.0
525.0 1000.0
580.0 1175.0
650.0 1130.0
1605.0 620.0
1220.0 580.0
1465.0 200.0
1530.0 5.0
845.0 680.0
725.0 370.0
145.0 665.0
415.0 635.0
510.0 875.0
560.0 365.0
300.0 465.0
520.0 585.0
480.0 415.0
835.0 625.0
975.0 580.0
1215.0 245.0
1320.0 315.0
1250.0 400.0
660.0 180.0
410.0 250.0
420.0 555.0
575.0 665.0
1150.0 1160.0
700.0 580.0
685.0 595.0
685.0 610.0
770.0 610.0
795.0 645.0
720.0 635.0
760.0 650.0
475.0 960.0
95.0 260.0
875.0 920.0
700.0 500.0
555.0 815.0
830.0 485.0
1170.0 65.0
830.0 610.0
605.0 625.0
595.0 360.0
1340.0 725.0
1740.0 245.0