元启发式算法：模拟退火算法原理&TSP搜索优化MATLAB实现

模拟退火算法

#单点搜索 #元启发式 #接受劣解 #全局最优

学习内容主要来自书籍《25个经典的元启发式算法-从设计到MATLAB实现》

金属退火过程&模拟

金属退火是一种热处理工艺：将金属缓慢加热到一定温度保持足够时间后以适宜速度(一般是缓慢)冷却，使得金属具备一些之前特性。

退火过程的粒子运动：如下图[1] 当固体的温度很高的时候，内能较大，固体的内部粒子处于快速无序运动，但随着温度的不断下降，粒子的热运动逐渐减弱且趋于有序，当温度降至结晶温度，粒子运动变为围绕晶体格点的微小振动，液体凝固成固体的晶态。冷却过程中，系统的熵值不断减小，系统能量也趋于最小值，由此可以认为，粒子运动本身就是一种寻优过程。

退火的两种过程

冷却过程：温度下降，粒子热运动减弱，冷却过程系统能量随温度降低趋于最小值。
等温过程：固体粒子在每个恒定温度下都可以达到一个平衡态，最终常温时达到固体的基态。

$Metropolis$准则

金属粒子在每一个恒定温度下都可以充分达到热平衡(金属内部各部分温度均匀),这个过程可以使用Monte Carlo方法模拟但是工作量巨大，因而$Metropolis$提出以一定概率接受新状态的方式模拟等温过程：某粒子初始状态为$i$,温度下降会导致粒子发生迁移，得到新状态$j$,粒子迁移前后能量风别为$E_i$和$E_j$,若$E_j<E_i$,即新状态能量相对较小，则无条件接受新状态成为当前状态，否则新状态的接受概率为$p$.
$$p=e^{-(E_j-E_i)/KT}$$
$$K为玻尔兹曼常数，T为温度，E_i-->E_j$$

上图为接收概率$p$的函数图像，其中$\Delta f=E_j-E_i$, 可见新状态在高温时容易被接受，随着T的下降,新状态的内接受的概率越来越小。

模拟退火算法(Simulated Annealing Algorithm,SA)

模拟“基于$Metropoils$的模拟退火过程”温度下降过程中粒子自发热运动以使系统状态朝着自由能减少并最终达到自由能最小的平衡态的方式/规律进行单点搜索寻找最优目标函数值的算法。
也就是说，模拟退火算法实际上模拟的是粒子的运动，寻优过程即粒子运动过程。

|       ~~~~~      退火过程| \qquad \qquad       ~~~~~      优化问题 \qquad |
| :—: | :----: | :----: |
| \qquad 自由能| \qquad \qquad 目标函数$f$(x)|
| \qquad 能量状态$E_i$| \qquad \qquad 可行解$x$|
| \qquad 最低自由能| \qquad \qquad 最优目标值$f^{\ast }$|
| \qquad 最低能量状态| \qquad \qquad 最优解$x^{\ast }$|
| \qquad 加温过程| \qquad \qquad 设定初始温度$T_0$|
| \qquad 等温过程| \qquad \qquad 当前$T$下展开邻域搜索并以$Metropolis$准则接受劣解|
| \qquad 冷却过程| \qquad \qquad $T$减小|

状态：粒子状态，能量/自由能：系统自由能

算法构成

• 设置初始温度$T_0$，即加温过程。

• 外循环通过温度$T$的下降控制迭代规模，对应冷却过程。

    > 利用温度下降控制迭代规模。常见降温函数：
     	- 比率法：T = T * tau, tau可取值0.95~0.99，tau越大温度下降越慢.
    	- 步长法：T = T - detaT.
  

• 内循环进行邻域搜索，对应等温过程，邻域搜索的目的是要达到该温度下的平衡态，理论上的平衡状态很难达到，所以内循环次数往往是一个常数。

    >  温度T固定下，对当前解进行邻域搜索获取一批邻域解，并按照Metropolis准则
    	接受劣解为新的当前解。
  

$Metropolis$准则接受劣解的意义

$$p=e^{-(E_j-E_i)/KT}$$
$$K为玻尔兹曼常数，T为温度，E_i-->E_j$$

• $变量T$

    >   如下图[2], 当T很大时,接受概率几乎为1，算法接受任何邻域解即使这个解比当前解
    要差，此时有利于广域搜索。
  

	> 如下图[2],当T很小时，接受概率趋于0，仅接受当前邻域中更好的解，
	此时有利于局域范围的精细搜索。

• 变量$\Delta f(x)$

    > 温度一定时， ∆𝒇越小，概率越大，即目标函数相差越小接受的可能性越大。
  

参数控制

• 能量状态表达/状态编码，即可行解的编码方式：编码方式直接决定领域的构造和大小，合理的状态编码可以减少计算复杂性。

    > 0-1编码适用于：0-1背包问题、指派问题.   
      自然数编码适用于：TSP问题.   
      实数编码适用于：连续函数优化问题.
  

• 初始温度$T_0$

    > 初始温度T_0的设定要适当地高，以满足广域搜索以及降温过程。
  

• 内循环次数

    > 内循环次数应足够多以达到平衡状态。
   	- 取常数：调参.
   	- 不取常数：当T较大时，内循环次数可以较少，当T较小时，内循环次数适当增大。
  

• 降温策略

    > 调参  
     降温过快，算法很快从广域搜索状态转换为局部搜索状态，可能导致过早陷入局部最优，  
     降温过慢，增加算法计算时间。
  

爬山vs模拟退火

爬山算法是一种十分简单的贪婪搜索算法，该算法每次从当前解的临近解空间中选择一个最优解作为当前解，直到达到一个局部最优解。

> 如上图[1]当前解是C点，局部搜索后可能找到G点，然后再从G点搜索，可能会找到D点，
  但是D点两侧的值都小于D点，爬山算法将认为D点为最大值，从而陷入了局部最优。
  
> 模拟退火在局部搜索即使找到一个比当前解差的解，也会以一定的概率跳到这个解上，
  则在一定程度上避免 陷入局部最优的情况.

> 爬山算法比喻：兔子朝着比现在高的地方跳去,它找到了不远处的最高山峰,但是这座山未必是珠穆朗
玛峰,这就是爬山算法，它不能保证局部最优值就是全局最优值.

> 模拟退火算法比喻：兔子喝醉了,它随机地跳了很长时间,这期间，它可能走向高处，也可能踏入平地,
   但是，它渐渐清醒了并朝最高方向跳去。这就是模拟退火.逐渐清醒的过程就是由广域搜索转向局部
   搜索的过程。

SA优化搜索TSP

城市数量17,自然数状态编码，初始温度1000，降温系数0.95，内循环固定20次，一次邻域搜索
获取邻域解数量同城市数量，邻域搜索调用一次$2-opt$算子(随机交换两种城市,O(n^2))
解空间：17! = 3.55687428×10^14
目标最优值：7291

代码（matlab，源代码见文末)

%SA2TSP_17.m
%外循环-降温过程，初始温度1000，降温系数0.95
while T0>=Ts
    %内循环-等温过程，固定迭代次数20次
    %path是当前路径，即当前走到了哪里、现在在山上哪个位置
    %进入内循环之前path就是上一个等温过程的平衡态
    %由于概率接受，path不一定等于历史最佳值bestPath
    for in = 1:MaxInnerLoop %20
        e0 = distance_calc(M,path);
        %对于每轮循环开始前，path就是经过上一轮走到的位置
        %每轮循环从开始，对path进行邻域搜索
        %使用2-opt(随机交换两个城市)，得到neiborNum(==cityNum,17)个邻域解，O(n^2)
        neiboorSols = Neiborhood(path,cityNum,neiborNum);
        e2 = distance_calc(M,neiboorSols);
        [better,index] = sort(e2);
        %取邻域搜索所得佳路径newpath及最佳值
        e1 = better(1,1);
        newpath = neiboorSols(index(1),:);
        %目标值改进，即比当前位置好，无条件接受为新的当前解path,并作为下一轮内循环的起点
        if e1 <e0
            currentbestValue = e1; %currentbestPath:一直代表当前位置
            currentbestPath = newpath;
            path = newpath;
            %如果新的当前位置比历史最佳位置好，更新bestPath
            if bestValue >currentbestValue
                bestValue = currentbestValue;
                bestPath = currentbestPath;
            end
        %否则按照Metropolis准则接受newpath为下一轮内循环起点path
        else
            pt = min(1,exp(-(e1-e0)/T0));
            if pt > rand
                path = newpath;  e0 = e1;
            end
        end
    end
    %等温过程结束，降温(缓慢下降)
    T0 = tau*T0;
end

> 外循环迭代次数：134(1000*0.95^134 = 1.03) 
> 调用O(n^2)的2-opt算子：2680次 
> 产生邻域解：45560个

• 运行100次:bestValue变化
  

  • 平均目标值：7319.5
  • 目标值标准差：58.0688
  • 平均时间：0.0483s

一次模拟退火过程当前目标值变化特点

• 记录每次迭代后(外循环),当前方案/位置/解对应的目标函数值变化：初始阶段搜索得到目标值波动大，但是伴随着循环次数的增加、温度的降低，呈现逐渐减小且波动变小的趋势，最后趋于稳定。
• 实际上80次之后就收敛了

Q: 前期广域搜索用大范围邻域结构、后期局部搜索使用小邻域结构会不会更好？？？

代码地址

SA2TSP_17

参考博文

[1]超详细–模拟退火算法（Simulated Annealing）有源码！！
[2]速通模拟退火(Simulated Annealing)