[人工智能]模拟退火的实现思路、搜索模型及其可视化

模拟退火的实现思路、搜索模型及其可视化

原理与特点

模拟退火(Simulated Annealing) 是一种启发式的优化算法，可以针对较复杂的解析式实现快速搜索其最值，优化目标函数。算法的实现过程与金属的退火过程类似。

退火 是指将金属加热到一定温度直到液化再冷却下来的过程，在退火过程中原子按照局部能量最小化进行排列。在冶金学中常利用金属的退火进行原子重排，以得到韧性和硬度都更好的金属。¹

模拟退火实质上是一种随机搜索算法，其核心在于温度这一变量。在退火初期温度较高时，目标函数的状态变量通过随机无序运动尽可能探索状态空间，寻找目标函数更优值，并一定概率接受更差状态以便于跳出局部最优解，而随着温度降低，探索幅度下降，状态变量趋于稳定。

模拟退火常常用于解当算法无法遍历整个状态空间去寻找目标函数的最优解时的问题，也就是组合爆炸(combination explosion) 的情况。这个时候可以采用贪心算法去快速寻找一个近似解，但同样也会有概率陷入局部最优解。模拟退火算法一方面能够通过随机搜索的方式节省算力，为NP复杂性问题提供有效的近似解。另一方面，也能实现通过温度这一变量实现尽可能广的探索，使得收敛结果能够跳出局部最优解。‘

算法实现思路

1.冷却过程

在模拟退火中，温度作为一个全局变量使用。首先需要明确温度下降函数。（温度T减小的模型也称为冷却进度表）

温度的下降一般有几何冷却和线性冷却两种方式：

线性冷却

$$T^{\prime }=T?\alpha \text{\hspace{1em}(1)}$$

几何冷却

$$T^{\prime }=T\alpha \text{\hspace{1em}(2)}$$

其中$\alpha$为0到1之间的数
一般在模拟退火算法中采用几何冷却的方式模拟温度的下降

退火代码总体框架

double T0=5000; //初始温度
double alpha=0.01; //每次退火1%  
while(T>1e-4)
{
            //退火过程 Process of annealing

T=T*（1-alpha）; //冷却过程 Process of cooling
}

2.退火过程

（1）选择起点

也就是初始化自变量$x_0$

（2）与温度正相关的随机搜索

随机搜索的幅度一般选择与温度线性相关，用公式表示如下

$$dx=Rand?T\text{\hspace{1em}(3)}$$
$$x_{new}=x_{old}+dx\text{\hspace{1em}(4)}$$

其中$dx$表示随机搜索的幅度，$Rand$为使得$x_{new}$生成在定义域范围内满足平均分布的随机数,可以设置为未生成在定义域内则重新生成，但算法往往会在生成规范的随机数上花一定的时间，降低了退火过程的效率，也可以通过设计搜索模型实现一次性生成满足要求的随机数。

随机搜索是模拟退火的关键，决定了最终解的精确度。其核心是: 满足在退火的前期随机搜索的程度较大使得尽可能搜索状态空间，优化目标函数并跳出局部最优解，而在退火的后期应该减小随机搜索程度使得自变量在最优解附近的空间搜索，使得解更加精确。

（3）接受更好的状态或以一定概率接受更坏的状态

退火算法中，跳出局部最优解的核心在于利用Metropolis抽样法确定对新的状态的选择准则。当系统能量状态将要变化时，系统能量由E(n)变化到E(n+1)的概率为
$$r=1,E(n+1)<E(n)\text{\hspace{1em}(5)}$$

$$r=e^{?\frac{E(n)?E(n+1)}{T}},E(n+1)>E(n)\text{\hspace{1em}(6)}$$

也就是系统状态向更有利的方向（能量更小）发展时，一定接受新的状态；如果向不利方向发展，则以一定概率($e^{?\frac{E(n)?E(n+1)}{T}}$)接受。
（E(n)可以看作就是目标函数f(x)，E(n)越小，代表f(x)更优，对应的x更优)

可以看出，温度影响了对新的状态接受概率，当温度较高时，可以接受较差情况的解，而当温度降低时，接收更差情况的解概率更小，以便于稳定当前解。

实例分析

例如利用MATLAB求解函数$f(x)=sin(x)+6sin(2x)$ 在（-3，3）之间的最大值。

当然在数学建模应用中，目标函数的实际解析式会比较复杂，难以通过遍历的方式求得最值时，可以用模拟退火得到解的近似值。

以下的代码均为MATLAB代码格式

（1）初始化设置

%目标函数： 在（-3，3）区间找到f(x)最大值对应的x
f=@(x)(sin(x)+3*sin(2*x)); 

%退火初始化设置
T0=12000;%初始温度
T=T0;
dT=0.01;%每次退火当前温度的1%
x=0;%退火起点
x_final=x;%最优解   精确解为x=0.8411

%随机数初始化设置
rng shuffle; %为每次计算设置不同的随机数种子

(2)与温度正相关的随机搜索

%退火过程
while(T>0.0001）
    %以温度的正相关概率随机移动（难点）
    dx=(((-3-x)/T0)+(6/T0)*rand())*(T0+0.8*(T-T0)*exp(-T));  
    %在((-3-x)/T0，(3-x)/T0)范围内的随机移动

这里有几种方法设计dx=f(T)的函数，此代码中所设计的函数为：
$$dx=[\frac{?3?x}{T_0}+\frac{6}{T_0}rand()]?[T_0+K(T?T_0)e^{?T}]\text{\hspace{1em}(7)}$$

T₀ 为初始温度，x 为当前的解，rand() 为(0,1)之间的随机数，T 为当前温度，K 为常数。

注：生成范围在（a,b）之间的随机数可以表示为$a+(b?a)?rand()$表示，例如$$\frac{?3?x}{T_0}+\frac{6}{T_0}rand()$$表示生成在（$\frac{?3?x}{T_0},\frac{3?x}{T_0}$）之间的随机数，通过此方法生成的$dx$有$x+dx$满足定义域（-3，3）范围

此搜索模型可以一次性生成在定义域范围内的随机数，提高退火效率，也能控制搜索的衰减过程。
参数K可以控制无序运动的衰减程度，设置合适的参数K可以更有效灵活地控制退火的搜索过程，使得在退火前期充分搜索退火后期处于稳定状态。注意此处K的取值为（0，1）

当然也可以利用重复生成与温度正相关的随机数，直到生成定义域范围的方法。代码如下：

  dx=inf;
  while(dx+x>b||dx+x<a) %当新生成的点在定义域外时
        dx=(2*rand()-1)*T; %重新生成随机数直到满足要求
  end

（3）接受更好的状态或以一定概率接受更坏的状态

如果随机搜索得到的新的自变量对应的目标函数$f(x_{new})$状态更优，则接受这个自变量作为当前解($x_{final}$)，如果新的状态$f(x_{new})$更差则以一定概率接受。

%接受更好的状态
    if(f(x)>=f(x_final))
        x_final=x;  
    %以一定概率接受更差的状态
    else
        df=f(x)-f(x_final); %计算df 即两者的状态差距大小
        if(rand()<exp(df/T))%（0，1）的随机数小于exp(df/T)的概率为exp(df/T)
            x_final=x;
        end
    end

（4）降温

    T=T*(1-dT);
    end

完整代码

%模拟退火
clear all

%目标函数： 在（-3，3）区间找到f(x)最大值对应的x
f=@(x)(sin(x)+3*sin(2*x)); 

%退火初始化设置
T0=12000;%初始温度
T=T0;
dT=0.01;%每次退火1%
x=0;%退火起点
x_final=x;%最优解   精确解为x=0.8411

%随机数初始化设置
rng shuffle; %为每次计算设置不同的随机数种子

%退火过程
while(T>0.0001)

    %以温度的正相关概率随机移动（难点）
    dx=(((-3-x)/T0)+(6/T0)*rand())*(T0+0.8*(T-T0)*exp(-T));  %在((-3-x)/T0，(3-x)/T0)范围内的随机移动
    x=x+dx;
    
    %接受更好的状态
    if(f(x)>=f(x_final))
        x_final=x;  
    %以一定概率接受更差的状态
    else
        df=f(x)-f(x_final); %计算df 即两者的状态差距大小
        if(rand()<exp(df/T))%（0，1）的随机数小于exp(df/T)的概率为exp(df/T)
            x_final=x;
        end
    end   
    
    %退火
    T=T*(1-dT);
end

求解结果

进行十次求解，得到的结果如下表所示。
注：精确解为x=0.8411

实现较准确的结果

（1）初始温度与终止温度

通过设置较高的初始温度$T_0$和尽可能接近0的终止温度提高搜索的精度。

（2）温度下降模型

一般采用几何冷却的模型，也就是温度T按照指数函数规律下降，使得在搜索后期能够在精确解附近进行搜索。

（3）搜索模型

关于随机搜索$dx=f(T)$有几种模型可供参考：

3.1衰减可控模型

$$dx=[\frac{a?x}{T_0}+\frac{b?a}{T_0}rand()]?[T_0+K(T?T_0)e^{?T}]\text{\hspace{1em}(8)}$$
$T_0$ 为初始温度，$x$为当前的解，$rand()$为(0,1)之间的随机数，模型满足$dx+x$（新的搜索点）的范围为[a,b]，$T$为当前温度，$K$ 为常数。

可以通过控制参数 $K$实现对搜索范围衰减的控制。以下是在初始温度$T_0=12000$、衰减率$dT=0.01$，终止温度为0.001时的测试结果。

1）当K=0.98时，随机运动的衰减程度较高，容易过早地收敛从而无法找到更精确的全局解。

2）当K=0.8时，能够实现退火前期充分搜索，后期在精确解附近搜索，从而准确收敛。

3）当K=0.5时，搜索过程衰减较小，能够实现对自变量状态空间较充分的搜索，不会陷入局部最优解，到了退火后期，当前解在精确解附近，但由于搜索幅度仍然较大，导致有一定几率无法搜索到精确解附近的状态空间。

实例求解中采用此模型

2.线性相关模型

$$dx=T?Rand\text{\hspace{1em}(9)}$$

采用重复生成的方法，直到生成满足定义域的随机数。大小与温度线性相关，由于温度是指数衰减，也就是实际上无序运动按照指数规律衰减。

  dx=inf;
  while(dx+x>b||dx+x<a) %当新生成的点在定义域外时
        dx=(2*rand()-1)*T; %重新生成随机数直到满足要求
  end

同样此模型会一定程度降低算法效率
但如果采用初始温度标幺模型：

$$dx=[\frac{a?x}{T_0}+\frac{b?a}{T_0}rand()]?T\text{\hspace{1em}(10)}$$

虽然能够使得生成的随机数一次性满足定义域$（a,b）$,但通过收敛表可以看到搜索过程随着温度衰减而过早地衰减，以至于还没有搜索到精确解附近的状态空间时就已经陷入了局部最优解。因此不推荐此类模型。

模拟退火可视化

实例：求解$f(x)=sin(x)+6sin(2x)$在（-3，3）之间的最大值

求解过程如图所示：

可视化之前，需要对每个搜索过程的变量保存到工作区，最后统一画图。在画图脚本里，可以使用matlab命令hold on和hold off，使得每一帧图像只包含函数图像和当前解，之前的解不再出现。然后通过储存每一帧图像合成动图gif格式。具体可参考命令：help(imwrite)。

可视化代码：
注：x_final_temp变量储存了每一次的当前解

%画动图
    figure(4)
    i=-3:0.01:3;
    
    %画动态图像
    for k=1:count-1
    
    %画整个函数图像
    plot(i,f(i));
    hold on        
    grid on 
    title(['温度Temperature=',num2str(T_temp(k)),"℃"])%题目
    
    %用.标记出当前解，大小为15
    plot(x_final_temp(k),f(x_final_temp(k)),'.','MarkerSize',15);
    legend('f(x)','当前解') %图注
    hold off
    
    %储存每一帧图像到文件
    frame = getframe(figure(4));
    im{k} = frame2im(frame);
    filename = 'SA.gif'; 
    [A,map] = rgb2ind(im{k},256);
    if k == 1
        imwrite(A,map,filename,'gif','LoopCount',Inf,'DelayTime',0);
    else
        imwrite(A,map,filename,'gif','WriteMode','append','DelayTime',0);
    end
    end 

当然，目标函数的自变量也可以拓展到二维，例如将$f(x)=sin(x)+6sin(2x)$绕Z轴旋转得到曲面，模拟退火求解过程类似，不再叙述。由于gif文件过大，这里不给出可视化。