[人工智能]模拟退火算法及实践

模拟退火算法

• 算法原理
• 算法伪码
• 解决TSP问题
• 算法可视化演示

算法原理

摸拟退火算法是基于随机搜索的，即在解的空间中展开随机搜索的。当问题的空间很大，而可行解比较多，并且对解的精度要求不高时，随机搜索是很有效的解决办法，因为其他的做法在这个时候时空效率不能让人满意。而借助演化思想和群集智能思想改进过的随机算法更是对解的分布有规律的复杂问题有良好的效果。

所谓退火是冶金专家为了达到某些特种晶体结构重复将金属加热或冷却的过程，该过程的控制参数为温度T。模拟退火法的基本思想是：在系统朝着能量减小的趋势这样一个变化过程中，偶尔允许系统跳到能量较高的状态，以避开局部极小点，最终稳定达到全局最小点。可以看到模拟退火不是单纯的采用贪心策略，它每获得一个解，对于该解有两种做法：若该解为更优解，则100%采纳；若该解为劣解，以一定的概率采纳该解，也就是说可能丢弃，可能采纳。所以在模拟退火算法的随机搜索过程中，当前的采纳解是时好时坏，呈现出一种不断波动的情况，但在总体的过程中又朝着最优的方向收敛。

评估解的好坏取决于问题的语境，例如旅行商TSP问题，我们用全部的距离加起来作为解的优劣情况。

算法伪码

（1）由一个产生函数从当前解产生一个位于解空间的新解；为便于后续的计算和接受，减少算法耗时，通常选择由当前新解经过简单地变换即可产生新解的方法，如对构成新解的全部或部分元素进行置换、互换等。

（2）计算与新解所对应的目标函数差。因为目标函数差仅由变换部分产生，所以目标函数差的计算最好按增量计算。

（3）判断新解是否被接受,判断的依据是一个接受准则，最常用的接受准则是Metropolis准则: 若ΔT<0则接受S′作为新的当前解S，否则以概率exp(-ΔT/T)接受S′作为新的当前解S。

（4）当新解被确定接受时，用新解代替当前解，这只需将当前解中对应于产生新解时的变换部分予以实现，同时修正目标函数值即可。此时，当前解实现了一次迭代。可在此基础上开始下一轮试验。而当新解被判定为舍弃时，则在原当前解的基础上继续下一轮试验。

//参考百度百科
s:=s0;e:=E(s)//设定目前状态为s0，其能量E(s0)
k:=0//评估次数k
while k<kmax and e>emax
//若还有时间（评估次数k还不到kmax）且结果还不够好（能量e不够低）则：
    sn:=neighbour(s)//随机选取一临近状态sn
    en:=Esn)//sn的能量为E(sn)
    if random() < P(e,en,temp(k/kmax)) then//决定是否移至临近状态sn
        s:=sn;e:=en//移至临近状态sn
    k:=k+1//评估完成，次数k加一
returns//回转状态s

模拟退火是一个不断迭代的过程，我们通过设定一个迭代次数来模拟时间，注意这个迭代次数，不能多也不能少，具体是通过实践出来的。如果设置得太小，那么模拟过程还没有收敛就已经结束了；而设置得太大，那么收敛之后的迭代式浪费时间，因为收敛之后已经不会再变了。（收敛就是达到目前的最优解状态）所以经过多次尝试调出一个恰当的迭代次数是非常有必要。

解决TSP问题

了解模拟退火的基本原理之后，实践是最好的学习方式。

TSP问题背景

TSP问题，假设一个旅行商人要去n个城市，给出n个城市的坐标，他必须经过且只经过每个城市一次，要求最后回到出发的城市，并且要求他选择的路径是所有路径中的最小值。

拟定算法

(1).随机生成一个城市序列作为初始解，比如1、2、… 140，这样的一个序列；设置合适的初温；
(2).通过交换两个城市的位置得到序列的领域，作为新解，如果温度为0，则转(6)；
(3).将新解与最优解比较，如果新解小于最优解，则将新解作为最优解,否则则以Metropolis 准则决定是否接受差解为最优解；
(4).如果系统处于平衡状态，则转(5)，否则接着执行(2)；
(5).降温，迭代计数器加1，返回(2)；
(6).输出最优解。

其中，我们给每个城市编号1-n，每一个解对应一个路径序列，代表通过这个路径走遍全部城市。我们评估函数就是这个路径的大小，最终目的就是尽可能找到一个路径长度最小的解。关于达到系统平衡状态，这个我们设置一个内置的循环，整个算法是双重循环，外循环表示退火过程，内循环迭代致使达到一个平衡状态。

编程实现

变量表示:

const int nCities = 130; //城市数量  
const double SPEED = 0.98;//退火速度  
const int INITIAL_TEMP = 1000;//初始温度  

struct node{//表示一个城市
    int num;//编号
    double x;//坐标
    double y;
}nodes[nCities];

class GenbyOne {//仿函数，用于初始化一个序列1,2,3,4,5.....,n
public:
    GenbyOne(int _seed = -1) : seed(_seed) {}
    int operator() () { return seed += 1; }
private:
    int seed;
};

double length_table[nCities][nCities];//存储各点间距离

解的状态表示:

class answer //一个解  
{
public:
    answer() {
        length = 0;
        generate(path, path + nCities, GenbyOne(0));//自动生成形如0，1，2，3的编号序列
        random_shuffle(path, path + nCities);//将序列随机打乱作为初始路径
        Calculate_length();
        //cout << " length->" << length;
    }

    void Calculate_length() {
        length = 0;
        //遍历path
        for (int i = 1; i < nCities; i++) {
            length += length_table[path[i - 1] - 1][path[i] - 1];
            //cout << length << endl;
        }
        length += length_table[path[nCities - 1] - 1][path[0] - 1];
    }

    double length;//代价，总长度  
    int path[nCities];//路径  
    //用于比较解之间的优劣
    bool operator < (const  answer &other) const
    {
        return length < other.length;
    }
};

设置初温:根据接收准则倒推初温

double genInitialTemp(){
     srand(time(NULL));
     int numStatus = rand()%50+50;
     double minLength = INT_MAX;
     double maxLength = INT_MIN;
     //随机num个状态
     for(int x = 0;x < numStatus;x ++){
        answer temp;
        minLength = (minLength > temp.length)?temp.length:minLength;
        maxLength = (maxLength < temp.length)?temp.length:maxLength;
     }
     //获取最大最小差异
     double delta = maxLength - minLength;
     double p = 0.8;
     double result = 0.0;
     //倒推
     result = -delta/(log(1.0/p - 1.0));
     return result;
}

根据当前解获取它的一个邻域解:

void getNewSolution(answer &t) {
    int i = rand() % nCities;
    int j = rand() % nCities;
    if (i > j) {
        swap(i, j);
    }
    else if (i == j)return;
    //随机取路径中两点进行操作

    int choose = rand() % 3;
    switch (choose) {
    case 0://交换
        swap(t.path[i], t.path[j]); break;
    case 1://逆转
        reverse(t.path + i, t.path + j); break;
    default://旋转
        if (j == (nCities - 1)) return;
        rotate(t.path + i, t.path + j, t.path + j + 1);
    }
    t.Calculate_length();
}

求解模拟:

answer SA_TSP(bool is_SA) {
    srand(time(0));//之后用到很多随机数，设定种子
    int i = 0;
    double r = SPEED;
    double t = genInitialTemp();//获取初温
    const double t_min = 0.01;
    answer temp;
    answer best = temp;
    while (t > t_min) {
        //设定一个固定的循环次数
        int L = 100 * nCities;
        for (i = 0; i < L; i++) {
            answer next = temp;
            getNewSolution(next);//获取邻域解
            if (Accept(next, temp, t)) {
                temp = next;//接受新解为当前最优
                if(next < best){//记录全局最优
                    best = next;
                }
            }
        }
        t *= r;//降温
    }
    return best;
}

接受解的函数:

bool Accept(answer &temp, answer &best, double t) {
    if (temp.length < best.length) return true;//优解必定接受
    else if (rand() / (double)RAND_MAX < 1 / (1 + exp(-(best.length - temp.length) / t)))
        return true;
    return false;
}

算法可视化演示

上面已经列出了算法的核心代码，我借用Qt的做出了一个可视化演算。

曲线显示的是当前解的路径距离，右边显示的当前解下不同城市的连线状态，可以看到摸拟退火整个过程是上下波动的，但整体是朝着一个减小的过程进行的。最终达到收敛后会得到一个优美的无交叉的曲线，与目前那个网站上算到的最优距离差距不到2%：

模拟退火算法是可以达到全局最优的，但是爬山法就不一定的，爬山法就是采用贪心的策略进行局部搜索，它得到结果就没那么好了，可以看到明显的线的交叉：

核心源码和测试数据下载（无GUI)：https://github.com/ZeusYang/AILearning
参考资料：《人工智能基础教程（第二版》作者：朱福喜