机器学习之算法优化(二)


一、模拟退火算法概念

模拟退火算法(SA)来源于固体退火原理,是一种基于概率的算法。将固体加温至充分高的温度,再让其徐徐冷却,加温时,固体内部粒子随温升变为无序状,内能增大,分子和原子越不稳定,从而可以自由运动和重新排序。而徐徐冷却时粒子渐趋有序,能量减少,原子越稳定。从高温开始缓慢降温(这个过程称为退火),粒子就可以在每个温度下达到热平衡,最后在常温时达到基态,此时内能最小。

模拟退火算法从某一较高初温出发,伴随温度参数的不断下降,结合概率突跳特性在解空间中随机寻找目标函数的全局最优解,即在局部最优解能概率性地跳出并最终趋于全局最优。模拟退火算法是通过赋予搜索过程一种时变且最终趋于零的概率突跳性,从而可有效避免陷入局部极小并最终趋于全局最优的串行结构的优化算法。

二、模拟退火算法原理

模拟退火算法包含两个部分即Metropolis算法和退火过程,,分别对应内循环和外循环

  • 外循环就是退火过程,将固体达到较高的温度(初始温度T(0)),然后按照降温系数alpha使温度按照一定的比例下降,当达到终止温度Tf时,冷却结束,即退火过程结束。
  • Metropolis算法是内循环,即在每次温度下,迭代L次,寻找在该温度下能量的最小值(即最优解)。

Metropolis算法数学模型:假设某材料在状态 i 之下的能量为E(i),那么材料在温度T时从状态 i 进入状态 j ,遵循以下规律:

  • (1)如果E(j)≤E(i),则接受该状态被转换,接受概率为1。

  • (2)如果E(j)>E(i),则状态转换以如下概率被接受:在这里插入图片描述

  • T为材料温度

Metropolis算法中能量和迭代次数简单关系图:Metropolis算法就是如何在局部最优解的情况下让其跳出来(如下图中B、C、E为局部最优),是退火的基础。1953年Metropolis提出重要性采样方法,即以概率来接受新状态,而不是使用完全确定的规则,称为Metropolis准则,计算量较低。
在这里插入图片描述
图中所示即为在一次温度下,跌代L次,固体能量发生的变化。在该温度下,整个迭代过程中温度不发生变化,能量发生变化,当前一个状态E(i)的能量大于后一个状态E(j)的能量时,状态E(i)的解没有状态E(j)的解好,所以接受状态为E(j)。但是如果下一状态的能量比前一个状态的能量高时,该不该接受下一状态呢?在这里设置一个接受概率P,即如果下一状态的能量比前一个状态的能量高,则接受下一状态的概率为P,下面具体讲一下如何接受下一个状态。

  1. 假设开始状态在A,多次迭代之后更新到B的局部最优解,这时发现更新到B时,能力比A要低,则说明接近最优解了,因此百分百转移,状态到达B后,发现下一步能量上升了,如果是梯度下降则是不允许继续向前的,而这里会以一定的概率跳出这个坑,这各概率和当前的状态、能量等都有关系。

  2. Metropolis算法数学模型我们可以看到,如果能量减小了,那么这种转移就被接受(概率为1),如果能量增大了,就说明系统偏离全局最优值位置更远了,此时算法不会立刻将其抛弃,而是进行概率操作:首先在区间[0,1]产生一个均匀分布的随机数ε,如果ε<P,则此种转移接受,否则拒绝转移,进入下一步,往复循环。其中P以能量的变化量和T进行决定概率P的大小,所以这个值是动态的。

  3. 用固体退火模拟组合优化问题,将内能E模拟为目标函数值f,温度T演化成控制参数 t,即得到解组合优化问题的模拟退火算法:由初始解 i 和控制参数初值 t 开始,对当前解重复“产生新解→计算目标函数差→接受或舍弃”的迭代,并逐步衰减 t 值,算法终止时的当前解即为所得近似最优解,退火过程由冷却进度表(Cooling Schedule)控制,包括控制参数的初值t及其衰减因子Δt、每个 t 值时的迭代次数L和停止条件 Tf。而温度的作用就是来计算转移概率P的。当温度每次下降后,转移概率也发生变化,因此在所有温度下迭代L次的结果也都是不相同的。在每个温度下迭代L次来寻找当前温度下的最优解,然后降低温度继续寻找,直到到达终止温度,即转移概率P接近于0。

Metropolis算法中接受状态的三条原则:

  • (1)在固定温度下,接受使目标函数下降的候选解的概率要大于使目标函数上升的候选解概率。
  • (2)随着温度的下降,接受使目标函数上升的解的概率要逐渐减小;
  • (3)当温度趋于零时,只能接受目标函数下降的解。

三、模拟退火算法参数控制

Metropolis算法是模拟退火算法的基础,但是直接使用Metropolis算法 可能会导致寻优速度太慢,以至于无法实际使用,为了确保在有限的时间收敛,必须设定控制算法收敛的参数,在上面的公式中,可以调节的参数就是T,T如果过大,就会导致退火太快,达到局部最优值就会结束迭代,如果取值较小,则计算时间会增加,实际应用中采用退火温度表,在退火初期采用较大的T值,随着退火的进行,逐步降低,具体如下:

  • (1)初始的温度T(0)应选的足够高,使的所有转移状态都被接受。初始温度越高,获得高质量的解的概率越大,耗费的时间也就越长。
  • (2) 退火速率。 最简单的下降方式是指数式下降:T(n)=λT(n),n=1,2,3……。 其中λ是小于1的正数,一般取值为0.8到0.99之间。
    其它下降方式有:
    在这里插入图片描述
  • ( 3)终止温度:如果温度下降到终止温度或者达到用户设定的阈值,则退火完成。

四、模拟退火算法步骤

1.基本思想

  • (1) 初始化:初始温度T(足够大),初始解状态S(算法迭代的起点),每个T值的迭代次数L

  • (2) 对k=1, …, L做第(3)至第6步:

  • (3) 产生新解S′

  • (4) 计算增量ΔT=C(S′)-C(S),其中C(S)为目标函数,C(S)相当于能量

  • (5) 若ΔT<0则接受S′ 作为新的当前解,否则以概率exp(-ΔT/T)接受S′作为新的当前解.

  • (6) 如果满足终止条件则输出当前解作为最优解,结束程序。

  • (7) T逐渐减少,且T->0,然后转第2步。

2.新解的产生和是否接受

模拟退火算法新解的产生和接受可分为如下四个步骤:

  • 第一步:由一个产生函数从当前解产生一个位于解空间的新解;为便于后续的计算和接受,减少算法耗时,通常选择由当前新解经过简单地变换即可产生新解的方法,如对构成新解的全部或部分元素进行置换、互换等,注意到产生新解的变换方法决定了当前新解的邻域结构,因而对冷却进度表的选取有一定的影响。

  • 第二步:计算与新解所对应的目标函数差。因为目标函数差仅由变换部分产生,所以目标函数差的计算最好按增量计算。事实表明,对大多数应用而言,这是计算目标函数差的最快方法。

  • 第三步:判断新解是否被接受,判断的依据是一个接受准则,最常用的接受准则是Metropolis准则:
    若ΔT<0则接受S′作为新的当前解S,否则以概率exp(-ΔT/T)接受S′作为新的当前解S。

  • 第四步:当新解被确定接受时,用新解代替当前解,这只需将当前解中对应于产生新解时的变换部分予以实现,同时修正目标函数值即可。此时,当前解实现了一次迭代。可在此基础上开始下一轮试验。而当新解被判定为舍弃时,则在原当前解的基础上继续下一轮试验。

3.退火算法程序流程图

在这里插入图片描述

五、模拟退火算法案例

举一个简单问题简述模拟退火算法过程:
在这里插入图片描述
式中,f ( x )为目标函数;g ( x ) 为约束方程;D为定义域,简单的模拟退火算法如下:

在这里插入图片描述

模拟退火的直观理解是:在一个给定的温度,搜索从一个状态随机地变化到另一个状态。每一个状态达到的次数服从一个概率分布。当温度很低时,以概率1停留在最优解。

案例一

在这里插入图片描述

import numpy as np

# 位置矩阵
coordinates = np.array([[565.0, 575.0], [25.0, 185.0], [345.0, 750.0], [945.0, 685.0], [845.0, 655.0],
                        [880.0, 660.0], [25.0, 230.0], [525.0, 1000.0], [580.0, 1175.0], [650.0, 1130.0],
                        [1605.0, 620.0], [1220.0, 580.0], [1465.0, 200.0], [1530.0, 5.0], [845.0, 680.0],
                        [725.0, 370.0], [145.0, 665.0], [415.0, 635.0], [510.0, 875.0], [560.0, 365.0],
                        [300.0, 465.0], [520.0, 585.0], [480.0, 415.0], [835.0, 625.0], [975.0, 580.0],
                        [1215.0, 245.0], [1320.0, 315.0], [1250.0, 400.0], [660.0, 180.0], [410.0, 250.0],
                        [420.0, 555.0], [575.0, 665.0], [1150.0, 1160.0], [700.0, 580.0], [685.0, 595.0],
                        [685.0, 610.0], [770.0, 610.0], [795.0, 645.0], [720.0, 635.0], [760.0, 650.0],
                        [475.0, 960.0], [95.0, 260.0], [875.0, 920.0], [700.0, 500.0], [555.0, 815.0],
                        [830.0, 485.0], [1170.0, 65.0], [830.0, 610.0], [605.0, 625.0], [595.0, 360.0],
                        [1340.0, 725.0], [1740.0, 245.0]])
# 点的总数量
num = coordinates.shape[0]
# 两点之间的距离矩阵
distance = np.zeros((num, num))
for i in range(num):
    for j in range(i, num):
        distance[i][j] = distance[j][i] = np.linalg.norm(coordinates[i] - coordinates[j])
# 三大超参数
delta = 0.99 # 温度差系数
t = 100      # 迭代次数
tk = 1       

solution_new = np.arange(num)
#valuenew = np.max(num)
solution_current = solution_new.copy()
value_current = 99000  # np.max这样的源代码可能同样是因为版本问题被当做函数不能正确使用,应取一个较大值作为初始值
#print(valuecurrent)
solution_best = solution_new.copy()
value_best = 99000  # np.max
result = []  # 记录迭代过程中的最优解

# 模拟退火
while t > tk:
    for i in np.arange(1000):
        # 下面的两交换和三角换是两种扰动方式,用于产生新解
        if np.random.rand() > 0.5:  # 交换路径中的这2个节点的顺序
            while True:  # 产生两个不同的随机数
                loc1 = np.int(np.ceil(np.random.rand() * (num - 1)))
                loc2 = np.int(np.ceil(np.random.rand() * (num - 1)))
                if loc1 != loc2:
                    break
            solution_new[loc1], solution_new[loc2] = solution_new[loc2], solution_new[loc1]
        else:  # 三交换
            while True:
                loc1 = np.int(np.ceil(np.random.rand() * (num - 1)))
                loc2 = np.int(np.ceil(np.random.rand() * (num - 1)))
                loc3 = np.int(np.ceil(np.random.rand() * (num - 1)))
                if (loc1 != loc2) & (loc2 != loc3) & (loc1 != loc3):
                    break
            # 下面的三个判断语句使得loc1<loc2<loc3
            if loc1 > loc2:
                loc1, loc2 = loc2, loc1
            if loc2 > loc3:
                loc2, loc3 = loc3, loc2
            if loc1 > loc2:
                loc1, loc2 = loc2, loc1
            # 下面的三行代码将[loc1, loc2)区间的数据插入到loc3之后
            tmplist = solution_new[loc1:loc2].copy()
            solution_new[loc1:loc3 - loc2 + 1 + loc1] = solution_new[loc2:loc3 + 1].copy()
            solution_new[loc3 - loc2 + 1 + loc1:loc3 + 1] = tmplist.copy()
        value_new = 0
        for j in range(num - 1):
            value_new += distance[solution_new[j]][solution_new[j + 1]]
        value_new += distance[solution_new[0]][solution_new[51]]
        if value_new < value_current:  # 接受该解
            value_current = value_new
            solution_current = solution_new.copy()

            if value_new < value_best:
                value_best = value_new
                solution_best = solution_new.copy()
        else:  # 按一定的概率接受该解
            if np.random.rand() < np.exp(-(value_new - value_current) / t):
                value_current = value_new
                solution_current = solution_new.copy()
            else:
                solution_new = solution_current.copy()
    t = delta * t # 缓慢降低温度
    result.append(value_best)
    # print (t) #程序运行时间较长,打印t来监视程序进展速度
    
print('最佳路线', solution_best)

#用来显示结果
plt.plot(np.array(result))
plt.ylabel("bestvalue")
plt.xlabel("t")
plt.show()

输出为·:

最佳路线 [ 0 21 30 17  2 16 20 41  6  1 29 22 19 49 28 15 45 43 33 34 48 31 44 18
 40  7  8  9 42 32 50 10 51 13 12 46 25 26 27 11 24  3  5 14  4 23 47 37
 36 39 38 35]

在这里插入图片描述

案例二

在这里插入图片描述
代码为:

# from __future__ import division
import numpy as np
import matplotlib.pyplot as plt
import math


# 定义目标函数
def aimFunction(x):
    y = x ** 3 - 60 * x ** 2 - 4 * x + 6 # 目标函数值
    return y

x = [i / 10 for i in range(1000)] # x随机生成100个
y = [0 for i in range(1000)]
for i in range(1000):
    y[i] = aimFunction(x[i]) # y[i]对应的值是目标函数的值

plt.plot(x, y, color='red') # 画图

# 模拟退火
def main():
    T = 1000  # 初始温度
    Tmin = 10  # 温度最小值
    # 定义域[0,100]内的最小值
    x = np.random.uniform(low=0, high=100)  # x从一个均匀分布[low,high)中随机采样
    k = 50  # 内循环次数为50
    y = 0  # 目标函数值初始为0
    t = 0  # 时间初始为0
    while T >= Tmin:
        for i in range(k):
            y = aimFunction(x)#  计算 y值
        # 通过变换函数在 x 的邻域中生成一个新的 x
            xNew = x + np.random.uniform(low=-0.055, high=0.055) * T
        # 新的 x满足条件则生成新的 y
            if (0 <= xNew and xNew <= 100):
                yNew = aimFunction(xNew)
                if yNew - y < 0:
                    x = xNew # 如果新的 y减去原来的y 值小于0则x就等于新的 x
                else:
                    # Metropolis准则
                    p = math.exp((y - yNew) / T)
                    r = np.random.uniform(low=0, high=1)
                    if r < p:
                        x = xNew
        t += 1
        # print(t)
        T = 1000 / (1 + t)  # 降温函数,也可使用T=0.9T
    return x

print('x的最小值:', main())
print('对应y在定义域[0,100]内的最小值:', aimFunction(main()))

输出为:

x的最小值: 40.077681414057
对应y在定义域[0,100]内的最小值: -32153.948192769945

在这里插入图片描述
对算法中添加温度差系数delta,并对最优解点进行特殊标记以便于我们在图上看:

import numpy as np
import matplotlib.pyplot as plt

# 目标函数
def inputfun(x):
    return x ** 3 - 60 * x ** 2 - 4 * x + 6 

initT = 1000 # 初始温度
minT = 10    # 温度下限
iterL = 1000 # 每个T值的迭代次数
delta = 0.95 # 温度衰减系数
k = 1

initx = 10*(2*np.random.rand()-1)# 初始X
nowt = initT # 现在温度=初始温度
print("初始解:",initx)


# 初始数据量 0到100 分为50份并且是均匀分
xx = np.linspace(0,100,50)
yy = inputfun(xx)
plt.figure()
plt.plot(xx,yy)  # 绘制函数
plt.plot(initx,inputfun(initx),'ob')  # 初始解位置

# 模拟退火算法寻找最小值过程
while nowt > minT:
    for i in np.arange(1,iterL,1):  # 迭代iterL次
        funVal = inputfun(initx)
        xnew = initx+(2*np.random.rand()-1)  # 产生新解
        if xnew>=0 and xnew<=100:
            funnew = inputfun(xnew)
            res = funnew-funVal  # 函数增量
            if res<0:
                initx = xnew
            else:
                p = np.exp(-(res) / (k*nowt))
                if np.random.rand() < p:
                    initx = xnew
    plt.plot(initx, inputfun(initx), 'oy')  # 不断优化的过程
#    print initx
#    print nowt
    nowt = nowt*delta  # 缓慢降低温度

print("最优解:",initx)
print("最优值:",inputfun(initx))
plt.plot(initx,inputfun(initx),'*r')
plt.show()
plt.show()

输出为:

初始解: 3.0849740486399635
最优解: 39.72706078370775
最优值: -32148.458827016995

在这里插入图片描述


总结

模拟退火算法的应用很广泛,可以高效地求解NP完全问题,如旅行商问题(Travelling Salesman Problem,简记为TSP)、最大截问题(Max Cut Problem)、0-1背包问题(Zero One Knapsack Problem)、图着色问题(Graph Colouring Problem)等等,但其参数难以控制,不能保证一次就收敛到最优值,一般需要多次尝试才能获得(大部分情况下还是会陷入局部最优值)。观察模拟退火算法的过程,具有以下主要优势

  • 迭代搜索效率高,并且可以并行化;
  • 算法中有一定概率接受比当前解较差的解,因此一定程度上可以跳出局部最优;
  • 算法求得的解与初始解状态S无关,因此有一定的鲁棒性;
  • 具有渐近收敛性,已在理论上被证明是一种以概率 i 收敛于全局最优解的全局优化算法。
  • 模拟退火算法与初始值无关,算法求得的解与初始解状态S(是算法迭代的起点)无关。
  • 1
    点赞
  • 3
    收藏
    觉得还不错? 一键收藏
  • 打赏
    打赏
  • 0
    评论

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包

打赏作者

小白只对大佬的文章感兴趣

你的鼓励将是我创作的最大动力

¥1 ¥2 ¥4 ¥6 ¥10 ¥20
扫码支付:¥1
获取中
扫码支付

您的余额不足,请更换扫码支付或充值

打赏作者

实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值