现代优化算法-模拟退火

现代优化算法-模拟退火

2020/1/19

简单理解模拟退火-求最大值

兔子们用酒将自己灌醉了，它们随机跳了很长的时间。在这期间，他们可能走向高处，也可能踏入平地。但是，随着时间的流逝，他们渐渐清醒了并朝着最高方向跳去。

---

模拟退火算法（SAA）来源于固体退火原理，将固体加温至充分高，再让其徐徐冷却，加温时，固体内部粒子随温升变为无序状，内能增大，而徐徐冷却时粒子渐趋有序，在每个温度都达到平衡态，最后在常温时达到基态，内能减为最小。模拟退火算法能够很好的解决NP问题。

模拟退火算法能够很好的解决NP问题。NP类问题是指一个复杂问题不能确定是否在多项式时间内找到答案，但是可以在多项式时间内验证答案是否正确。NP类问题数量很大，如完全子图问题、图着色问题、旅行商TSP问题等。

---

Metropolis算法

假设材料在状态$i$下的能量为$E(i)$，那么材料在温度$T$时从状态$i$进入状态$j$就遵循如下规律：

（1）如果$E(j)\le E(i)$，则接受该状态被转换。

（2）如果$E(j)>E(i)$，则状态转换被概率$e^{\frac{E（i）-E（j）}{KT}}$接受。

其中$K$是物理学中的玻尔兹曼常数。T是材料温度。

---

模拟退火算法

• 令$T=T_0$，随机生成一个初始解$x_0$。

• 令$T$等于冷却进度表中的下一个值$T_i$。（温度T的变化服从某种递减函数）。

• 根据当前解$x_i$进行扰动，产生一个新的解$x_j$。计算相应的目标函数值$f(x_j)$。

• $f(x_j)<f(x_i)$，则接受新解。

  $f(x_j)>f(x_i)$，若$rand<p$，则接受新解$x_j$。$rand$是随机生成数。

  由Metropolis准则：
  $$p=\{\begin{array}{r}1,f(x_j)<f(x_i)\\ e^{-\frac{f(x_i)-f(x_j)}{T}},f(x_j)>f(x_i)\end{array}$$

扰动的方式有许多种，可以自行确定。

• 在温度$T_i$下，重复$L_k$次的扰动和接受过程。也就是上述过程中的3、4.
• 判断T是否已经达到$T_f$，如果是，则终止算法。否则继续下一个$T$值。

---

符号说明

符号	意义	备注
$T_0$	冷却开始时的温度
$T(t)$	温度衰减函数	将连续的温度点转换为离散点
$T_f$	控制参数T的终值	停止准则
$L_k$	马尔科夫链长度	任何一个温度的迭代次数

---

温度衰减函数$T(t)$

常用的简单下降方式：
$$t_{k+1}=\alpha t_k，其中（0<\alpha <1）$$
经典模拟退火算法的降温方式：
$$T(t)=\frac{T_0}{1g(1+t)}$$
快速模拟退火算法的降温方式：
$$T(t)=\frac{T_0}{1+\alpha t}$$

内、外循环终止准则（最大迭代次数）

内循环

也称Metropolis抽样稳定准则。用于决定在各温度下产生候选解的数目。

（1）固定长度：在每一个温度迭代相同的步数，步数的选取同问题的大小有关。

（2）随着温度的下降，将同一温度的迭代步数增加。（由接受和拒绝的比率来控制迭代步数。）

外循环

零度法：设置终止温度的阀值。很小的正数$\xi$。如0.000001.

循环总数控制法：设置外循环迭代次数。$T_{max}$=100；

不改进规则的控制法：算法搜索到的最优值连续若干步保持不变。

---

总结

整个优化过程就是不断寻找新解和缓慢降低温度的交替过程。

最终的解是对该问题寻优的结果。

缺点：如果降温速度过缓，求解过程因为算法太慢而效率减小。相对于简单的搜索算法不具有明显的优势。如果降温速度过快，很可能最终得不到全局最优解。

• 要确定在每一温度下状态转换的结束准则。 也就是零度法或者不改进规则法。

• 初始温度的选择还有确定某个可行解的邻域的方法也要恰当。

---

题目解析

题目原文：https://blog.csdn.net/astandoffish/article/details/77189620

我方有一个基地，经度和纬度为（70,40）。假设我方飞机的速度为1000公里/小时。我方派一架飞机从基地出发，侦察完敌方所有目标，再返回原来的基地。在敌方每一目标点的侦察时间不计，求该架飞机所花费的时间（假设我方飞机巡航时间可以充分长）。

代码符号	含义	备注
$d_{ij}$	点$i$与点$j$的距离	对称矩阵
$path$	路径	各点排列组成的路径
$long$	总路程	同时作为判断条件，判断新解是否被接受
$dj$	代价函数值	代价函数计算的是变动路程$\Delta d$

$$S=[1,2,3,\cdots ,102]$$

是所有路径点标签的集合。

（共有100个敌方基地，飞机需要从我方基地出发与降落，因此初始点1和中止点102均为我方基地。）

而每一个路径点都可能会被第二次，第三次,…,第101次观察。记${\pi }_i=j$ 表示第i次经过j点，$(j\in 1,2,\cdots ,102)$

其中，初始点和终止点${\pi }_1,{\pi }_{102}$必等于1和102。

产生初始解

程序用蒙特卡洛法产生初始解，具体操作见代码部分。

产生新解

随机生成两个序号$u、v$

现有解：
$${\pi }_1\cdots {\pi }_{u-1}{\pi }_u{\pi }_{u+1}\cdots {\pi }_{v-1}{\pi }_v{\pi }_{v+1}\cdots {\pi }_{102}$$
使用两点变换法，即颠倒$u,v$之间的路径顺序：
$${\pi }_1\cdots {\pi }_{u-1}{\pi }_v{\pi }_{v-1}\cdots {\pi }_{u+1}{\pi }_u{\pi }_{v+1}\cdots {\pi }_{102}$$
这样便产生了新的解。注意，解的含义是路径。

退火过程

产生新解的过程中，实际距离变化的只有$u、v$改变的边界处变动。

因此，代价函数可以表示为：
$$\Delta d=(d_{{\pi }_{u-1}{\pi }_v}+d_{{\pi }_u{\pi }_{v+1}})-(d_{{\pi }_{u-1}{\pi }_u}+d_{{\pi }_v{\pi }_{v+1}})$$
通过这个代价函数，结合metropolis接受准则，判断是否降温。

最终温度达到设定的值（极小），算法结束。

这里给出的是python代码，参考了matlab代码修改而成。

import pandas as pd
import matplotlib.pyplot as plt
import math
import numpy as np
import random
import matplotlib.pyplot as plt

pp = open('路径数据文件')
data = pd.read_table(pp,header=None,sep=' ')

df = pd.DataFrame(data)

# ### 转换为单列，便于合并

x = df.iloc[:, 0:8:2]
x = pd.DataFrame(x.values.reshape(-1,1))
y = df.iloc[:, 1:8:2]
y = pd.DataFrame(y.values.reshape(-1,1))

sj = pd.concat([x,y],axis=1,ignore_index=True)
sj = sj.rename(columns=({0:'x',1:'y'}))
sj = sj.values

loc0 = [70,40]
sj = np.insert(sj, 0 ,values=loc0,axis=0)
sj = np.insert(sj, len(sj) ,values=loc0,axis=0)
#增加初始列


length = len(sj)
length


sj=sj*math.pi/180 #角度化弧度


#距离计算
d = np.zeros([length,length])
for i in range(0,length-1):
    for j in range(i+1,length):
        temp = math.cos(sj[i,0]-sj[j,0])*math.cos(sj[i,1])*math.cos(sj[j,1])+math.sin(sj[i,1])*math.sin(sj[j,1])#计算余弦值
        d[i,j]=6370*math.acos(temp)
d = d+d.T #做成对称阵


path = []
long = np.inf
for i in range(1000): #蒙特卡洛法设定初始解
    array_rand = list(range(2,length))
    np.random.shuffle(array)
    path_0 = [1]+array_rand+[length]
    #产生介于1与长度之间的随机整数。
    temp = 0
    for j in range(length-2):
        temp+=d[path_0[j],path_0[j+1]]
    if temp<long:
        path = path_0
        long = temp
e=0.1**30 #极小值（零度法）
L=200000 #（迭代次数）
at=0.999
T=1; 


# ### 退火过程


#产生新解
for i in range(L):
    c = 1+np.floor((length-2)*np.random.rand(1,2))
    c = c.astype(np.int)
    c.sort()
    c1 = c[0,0]
    c2 = c[0,1]
    #计算代价函数值
    dj=d[path[c1-1]-1,path[c2]-1]+d[path[c1]-1,path[c2+1]-1]-d[path[c1-1]-1,path[c1]-1]-d[path[c2]-1,path[c2+1]-1]
    #接受准则
    if dj<0:
        tt = path[c1:c2+1]
        tt = tt[::-1]
        path=path[0:c1]+tt+path[c2+1:length]#将中间段倒置得到新路径
        long = long + dj
        T = T*at #接受并降温
    elif np.exp(-dj/T)>np.random.rand(1):#路径没有变短，以一定概率接受
        tt = path[c1:c2+1]
        tt = tt[::-1]
        path=path[0:c1]+tt+path[c2+1:length]
        long = long + dj
        T = T*at #接受并降温
    if T<e:
        break


xx = []
yy = []
for i in path:
    xx.append(sj[i-1,0])
    yy.append(sj[i-1,1])


plt.figure(figsize=(8,8),dpi=200)
plt.style.use("ggplot")  
plt.title('Total distance: '+ str(long))
plt.plot(xx, yy,'o-')
plt.show()

输出结果如下：