Python手把手构建模拟退火算法（SA）实现最优化搜索

2019年11月6日 0条评论 12,056次阅读 21人点赞

《Python手把手构建模拟退火算法（SA）实现最优化搜索》

简介

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

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

如果E(j)<E(i)，接受该状态被转换；
如果E(j)>=E(i)，则状态会以一定的概率被转换，转换的概率为e^{(E(i)-E(j))/KT}。其中K是物理学中的玻尔兹曼常数，T是材料的温度。

可以注意到，这个过程所得到的一个新状态E(j)完全依赖于前一个状态E(i)，和更前面的状态无关，因此这是一个马尔科夫过程。

在模拟退火算法中应注意以下问题：

理论上，降温过程要足够缓慢，要使得在每一温度下达到热平衡。但在计算机实现中，如果降温速度过缓，所得到的解的性能会较为令人满意，但是算法会太慢，相对于简单的搜索算法不具有明显优势。如果降温速度过快，很可能最终得不到全局最优解。因此使用时要综合考虑解的性能和算法速度；
要确定在每一温度下状态转换的结束准则。实际操作可以考虑当连续m次的转换过程没有使状态发生变化时结束该温度下的状态转换，或者按照最终温度达到某一目标低温，或连续几个温度下转换过程没有使状态发生变化算法就结束。

案例

我们还是以一个简单的例子来讲解模拟退火的算法实现。

假设现在我们有一个函数，形式为y=3x²-60x+9，我们可以观察下它的图像：

def x_function(x):
    return 3*x**2 - 60*x + 9


x = [i for i in np.linspace(0, 100)]
y = map(x_function, x)
plt.plot(x, list(y))
plt.show()

def x_function(x):

return 3*x**2 - 60*x + 9

x = [i for i in np.linspace(0, 100)]

y = map(x_function, x)

plt.plot(x, list(y))

plt.show()

得到如下图形：

《Python手把手构建模拟退火算法（SA）实现最优化搜索》

通过求导我们也能发现，在x=10处该函数取到最小值，有全局最优解。

接下来我们一步步建立模拟退火算法实现整个寻优过程。

（1）解空间

解空间就是我们所要求的定义域的空间范围，在这里我们定为[0, 100]。

（2）目标函数

目标函数就是我们要求的函数，这里即为y=3x²-60x+9。

（3）新解的产生

通过当前的解产生一个新解的方法有很多，在这里我们直接通过加上一个微小的偏差bias来微调这个值。

x_new = x + np.random.uniform(-1, 1)

1	x_new = x + np.random.uniform(-1, 1)

增加一个[-1,1]之间的实数。

（4）代价函数差

代价函数差就是前后两次函数值的差值：E(j)-E(i)。

（5）接受准则

接受准则就是算法最核心的地方：

if E(j) <E(i)，概率P=1;
else, 概率P=e^{(E(i)-E(j))/KT}。

如果E(j) <E(i)，就接受新的值；否则，以一定的概率接受新的值，即概率大于0-1之间的随机数则接受。在实际使用当中，我们可以把K看作1处理，甚至温度T都可以是任意尺度的值。比如初始T=1。

（6）降温

利用选定的降温系数a进行降温处理，T=a*T，从而得到一个新的温度。比如a=0.999。

（7）结束条件

可以选定一个结束的温度，当温度T不断衰减到某个值时，算法结束，输出当前状态。比如std=0.0000001。

程序实现

我们编写如下的代码：

# 初始温度
T = 1
x = np.random.uniform(0, 100)
# 终止温度
std = 0.00000001
# 衰减率
a = 0.999
while T > std:
    y = x_function(x)
    # 新值通过扰动产生
    x_new = x + np.random.uniform(-1, 1)
    if 0 <= x_new <= 100:
        y_new = x_function(x_new)
        if y_new < y:
            x = x_new
        else:
            p = np.exp((y - y_new) / T)
            r = np.random.uniform(0, 1)
            if p > r:
                x = x_new
                # print(x)
    T = T * a
print(x, x_function(x))

# 初始温度

T = 1

x = np.random.uniform(0, 100)

# 终止温度

std = 0.00000001

# 衰减率

a = 0.999

while T > std:

y = x_function(x)

# 新值通过扰动产生

x_new = x + np.random.uniform(-1, 1)

if 0 <= x_new <= 100:

y_new = x_function(x_new)

if y_new < y:

x = x_new

else:

p = np.exp((y - y_new) / T)

r = np.random.uniform(0, 1)

if p > r:

x = x_new

# print(x)

T = T * a

print(x, x_function(x))

最后输出结果为：

9.999874115470941 -290.99999995245935

1	9.999874115470941 -290.99999995245935

可见，我们通过模拟退火的算法寻找到了全局最优的最小值x=10。

总结

怎么来理解这个过程？事实上就是在解空间中先随机的选择一个解x，计算它的函数值，然后以一定的方式（可以是增减一个扰动）得到一个新的解x1，比较这两个函数值，取最小值。即使后者的函数值比前者大，也可以以一定的概率保留，只是说这个概率定义成了玻尔兹曼分布中的概率形式。和其他的启发式算法（遗传算法）一样，它们都是以概率为导向的迭代算法，说白了这个世界就是这么奇妙，冥冥之中自有概率在安排。

简介

案例