模拟退火算法

算法简介

模拟退火算法得益于材料统计力学的研究成果。统计力学表明材料中粒子的不同结构对应于粒子的不同能量水平。在高温条件下，例子的能量较高，可以自由运动和重新排列。在低温条件下，例子能量较低。如果从高温开始，非常缓慢地降温(这个过程被称为退火)，粒子就可以在每个温度下达到热平衡。当系统完全被冷却时，最终形成处于低能状态的晶体。

如果用粒子的能量定义材料的状态，Metropolis算法用一个简单的数学模型描述了退火过程。假设材料在状态\(i\)之下的能量为\(E(i)\),那么材料在温度\(T\)时从状态\(i\)进入状态\(j\)就遵循如下规律：

(1)如果\(E(j)\leq E(i)\),则接受该状态被转换。

(2)如果\(E(j) > E(i)\),则状态转换以如下概率被接受：

\[e^{\frac{E(i)-E(j)}{KT}}, \]

式中：K为物理学中的玻尔兹曼常数；T为材料温度。

在某一个特定温度下，进行了充分的转换之后，材料将达到热平衡。这时材料处于状态\(i\)的概率满足玻尔兹曼分布

\[P_T(x=i) = \frac{e^{- \frac{E(i)-E(j)}{KT}}}{\sum \limits_{j\in S}e^{\frac{-E(j)}{Kt}}} \]

式中: X为材料当前状态的随机变量；S为状态空间集合。显然

\[ \mathop{lim} \limits_{T\rightarrow \infty} \frac{e^{-\frac{E(i)}{KT}}}{\sum \limits _ {j\in S}e^{-\frac{E(j)}{KT}}} = \frac{1}{|S|} \]

式中：|S|为集合S中状态的数量。

这表明所有状态在高温下具有相同的概率。而当温度下降时，有

\[\begin{aligned} \mathop{lim} \limits_{T\rightarrow \infty} \frac{e^{-\frac{E(i)-E_{min}}{KT}}}{\sum \limits _ {j\in S}e^{-\frac{E(j)-E_{min}}{KT}}} &= \mathop{lim} \limits_{T\rightarrow \infty} \frac{e^{-\frac{E(i)-E_{min}}{KT}}}{\sum \limits _ {j\in S_{min}}e^{-\frac{E(j)-E_{min}}{KT}} + \sum \limits _ {j\notin S_{min}}e^{-\frac{E(j)-E_{min}}{KT}}}\\ & = \mathop{lim} \limits_{T\rightarrow \infty} \frac{e^{-\frac{E(i)-E_{min}}{KT}}}{\sum \limits _ {j\in S_{min}}e^{-\frac{E(j)-E_{min}}{KT}}}\\ &= \begin{cases} \frac{1}{|S_{min}|}, & i \in S_{min} \\ 0, & 其它 \end{cases} \end{aligned} \]

式中：\(E_{min} = \mathop{min}\limits _{j\in S} E(j)\) 且\(S_{min} = |i|E(i) = E_{min}|\)

上式表明当温度降至很低时，材料会以很大概率进入最小能量状态。

假定要解决的问题是一个寻找最小值的优化问题。将物理学中模拟退火的思想应用于优化问题就可以得到模拟退火寻优方法。

在每个\(T_{i}\)下，所得到的一个新状态\(x(k+1)\)完全依赖于前一个状态\(x(k)\),和前面的状态\(x(0),x(1),...,x(k-1)\)无关，因此这是一个马尔可夫过程。如果温度下降十分缓慢，而在每个温度都有足够多次的状态转移，使之在每一个温度下达到热平衡，则全局最优解将以概率1被找到。因此可以说模拟退火算法可以找到全局最优解。

应用举例

已知100个目标的经度、维度如表(省略)，我方有一个基地，经度和纬度为(70,40)。假设我方飞机的速度为\(1000km/h\)。我方派一架飞机从基地出发，侦察完所有目标，再返回原来的基地。在每一目标点的侦察时间不计，求该架飞机所花费的时间(假设我方飞机巡航时间可以充分长)。

这是一个旅行商问题。给我方基地编号为1，目标依次编号为\(2,3,...,101\),最后我方基地再重复编号为102（这样便于程序中计算）。距离矩阵\(D = (d_{ij})_{102\times102}\),其中\(d_{ij}\)表示\(i,j\)两点的距离，\(i,j = 1,2,\dots, 102\),这里\(D\)为实对称矩阵。则问题是求一个从点1出发，走遍所有中间点，到达点102的一个最短路径。

上面问题中给定的是地理坐标(经度和纬度)，必须求两点间的实际距离。设A，B两点的地理坐标分别为\((x_1,y_1),(x_2,y_2)\),过\(A,B\)两点的大圆的劣弧长即为两点的实际距离。以地心为坐标原点O，以赤道平面为XOY平面，以0度经线圈所在的平面为XOZ平面建立三维直角坐标系。则A,B两点的直角坐标分别为

\[A(Rcosx_1cosy_1,Rsinx_1cosy_1,Rsiny_1),\\ B(Rcosx_2cosy_2,Rsinx_2cosy_2,Rsiny_2), \]

式中：R = 6370km为地球半径。

A,B两点的实际距离

\[d = Rarcos(\frac{OA\cdot OB}{|OA| \cdot |OB|}), \]

化简，得

\[ d = Rarccos[cos(x_1-x_2)cosy_1cosy_2+siny_1siny_2] \]

求解的模拟退火算法描述如下：
(1)解空间。解空间S可表示为{1,2,...,102}的所有固定起点和终点的循环排列集合，即\(S = {(\pi_1,\pi_2,...,\pi_{102})}\),这里先使用蒙特卡洛方法求得一个较好的初始解。

（2）目标函数。目标函数（或称代价函数）为侦察所有目标的路径长度。要求

\[ min f(\pi_1, \pi_2, ... , \pi_{102}) = \sum_{i=1}^{101}d_{\pi_i \pi_{i+1}} \]

而一次迭代由下列三部构成。
（3）新解的产生。设上一步迭代的解为\(\pi_1...\pi_{u-1} \pi _u \pi _{u+1} ... \pi _{v-1} \pi _v \pi _{v+1} ... \pi _{w-1} \pi _w \pi _{w+1} ... \pi _{102}\)

①2变换法 任选序号\(u,v\)和\(w\),将\(u\)和\(v\)之间的路径插到\(w\)之后，对应的新路径为

\[ \pi_1...\pi_{u-1} \pi _v \pi _{v-1} ... \pi _{u+1} \pi _u \pi _{v+1} ... \pi _{102} \]

②3变换法。任选序号\(u,v\)和\(w\),将\(u\)和\(v\)之间的路径插到\(w\)之后，对应的新路径为

\[ \pi_1...\pi_{u-1}\pi _{v+1} ... \pi _w \pi _{u} ... \pi _{v} \pi _{w+1} ... \pi _{102} \]

（4） 代价函数差。对于2变换法，路径差可表示为

\[ \Delta f = (d_{\pi_{u-1}}{\pi _v} + d_{\pi_u \pi_{v+1}}) - (d_{\pi_{u-1}\pi_u +} d_{\pi_{v}\pi_{v+1}}) \]

(5)接受准则

\[P = \begin{cases} 1, & \Delta f< 0,\\ exp(-\Delta f/T), & \Delta f \geq 0 \end{cases} \]

如果\(\Delta f < 0\)，则接受 新的路径；否则，以 概率\(exp(-\Delta f/T)\)接受新的路径，即用计算机产生一个\([0,1]\)区间上均匀分布的随机数\(rand\),若\(rand \leq exp(-\Delta f/T)\)则接受。

(6)降温。利用选定的降温系数$ \alpha $ 进行降温，取新的温度\(T\)为\(\alpha T\)(这里T为上一步迭代的温度)，这里选定\(\alpha = 0.999\)。

(7)结束条件。用选定的中止温度\(e = 10^{-30}\),判断退火过程是否结束。若\(T < e\),则算法结束，输出当前状态。

代码：

sj0 = load('data12_1.txt');
x = sj0(:,[1:2:8]); x = x(:);
y = sj0(:,[2:2:8]); y = y(:);
sj = [x y]; d1 = [70,40];
xy = [d1;sj;d1];
xy = xy * pi /180;%角度化成弧度
d = zeros(102); %距离矩阵d初始化
for i = 1 : 101
    for j = i + 1:102
        d(i,j) = 6370 * acos(cos(sj(i,1)-sj(j,1)) * cos(sj(i,2)) * cos(sj(j,2)) + sin(sj(1,2)) * sin(sj(j,2)));
    end
end
d = d + d';
path = []; long = inf;%巡航路径及长度初始化
for j = 1 : 1000 %求较好的初始解
    path0 = [1 1+randperm(100),102]; temp = 0;
        for i = 1 : 101
            temp = temp + d(path0(i), path0(i+1));
        end
        if temp < long
            path = path0;
            long = temp;
        end
end
e = 0.1^30; L = 20000; at = 0.999; T = 1;
for k = 1 : L % 退火过程
    c = 2 + floor(100 * rand(1,2)); %产生新解
    c = sort(c); c1 = c(1); c2 = c(2);
    %计算代价函数值的增量
    df = d(path(c1-1),path(c2)) + d(path(c1),path(c2+1))-d(path(c1-1),path(c1))...
    -d(path(c1-1), path(c1)) - d(path(c2),path(c2+1));
    if(df < 0)          %接受准则
        path = [path(1:c1-1), path(c2:-1:c1), path(c2+1 : 102)];
        long = long + df;
    elseif exp(-df/T)>= rand
            path = [path(1:c1-1), path(c2 : -1 : c1), path(c2 + 1 : 102)];
            long = long + df;
    end
    T = T * at;
    if(T < e)
        break;
    end
end
path, long 
xx = xy(path,1);    yy = xy(path,2);
plot(xx, yy, '-*')