模拟退火算法

模拟退火算法

一、算法简介

模拟退火算法(Simulated Annealing,SA)是一种模拟物理退火的过程而设计的随机优化算法,结合爬山法和随机行走算法,同时避免算法进入局部最优,早期用于组合优化,后来发展成一种通用的优化算法。它的基本思想最早在1953年就被Metropolis提出,但直到1983年Kirkpatrick等人才设计出真正意义上的模拟退火算法并进行应用。

该算法采用类似于物理退火的过程,先在一个高温状态下(相当于算法随机搜索),然后逐渐退火,在每个温度下(相当于算法的每一次状态转移),徐徐冷却(相当于算法局部搜索),最终达到物理基态(相当于算法找到最优解)。

高温过程——增强粒子的热运动,使其偏离平衡位置,目的是消除系统原先可能存在的非均匀态;

等温过程——退火过程中要让温度慢慢降低,在每一个温度下要达到热平衡状态,对于与环境换热而温度不变的封闭系统满足自由能较少定律,系统状态的自发变化总是朝自由能减少的方向进行,当自由能达到最小时,系统达到平衡态;

冷却过程——使粒子热运动减弱并渐趋有序,系统能量逐渐下降,从而得到低能的晶体结构。当液体凝固为固体的晶态时退火过程完成。

因此模拟退火算法从某一高温出发,在高温状态下计算初始解,然后以预设的邻域函数产生一个扰动量,从而得到新的状态,即模拟粒子的无序运动,比较新旧状态下的能量,即目标函数的解。如果新状态的能量小于旧状态,则状态发生转化;如果新状态的能量大于旧状态,则以Metropolis接受准则发生转化。当状态稳定后,便可以看作达到了当前状态的最优解,便可以开始降温,在下一个温度继续迭代,最终达到低温的稳定状态,便得到了模拟退火算法产生的结果

该算法的关键点如下:

1、对固体退火过程的模拟;
2、采用Metropolis接受准则;
3、用冷却进度表控制算法进程,使算法在多项式时间里给出一个近似解。
固体退火过程是SAA的物理背景;Metropolis接受准则使算法跳离局部最优 “险井”;而冷却进度表的合理选择是算法应用的前提。

1、固体退火过程

固体退火过程的数学表述:在温度T,分子停留在状态r满足Boltzmann概率分布

img

温度低时能量低的微观状态概率大,温度趋于零时,固体几乎处于概率最大、能量最小的基态。在同一个温度T,选定两个能量E1<E2,有:

img

因此,在同一个温度,分子停留在能量小的状态的概率比停留在能量大的状态的概率要大。

2、Metropolis准则

Metropolis准则(1953)——以概率接受新状态,固体在恒定温度下达到热平衡的过程可以用Monte Carlo方法(计算机随机模拟方法)加以模拟。

Monte Carlo模拟退火过程:蒙特卡罗(Monte Carlo)方法,或称计算机随机模拟方法,是一种基于“随机数”的计算方法。

在这里插入图片描述

可以看出,高温下可接受与当前状态能差较大的状态为重要状态,而在低温下只能接受与当前状态能差较小的新状态为重要状态。在温度趋于零时,就不再接受Ej>Ei的新状态j了.如此反复,达到系统在此温度下的热平衡。这个过程称作Metropolis抽样过程。上面这种方法被称为重要性采样法,该方法以概率形式接受新状态,避免局部最优。

3、冷却进度表

在Metropolis抽样过程中温度T缓慢的降低,模拟退火过程就是通过T参数的变化使状态收敛于最小能量处。因而,T参数的选择对于算法最后的结果有很大影响。初始温度和终止温度设置的过低或过高都会延长搜索时间。降温步骤太快,往往会漏掉全局最优点,使算法收敛至局部最优点。降温步骤太慢,则会大大延长搜索全局最优点的计算时间,从而难以实际应用。因此选择一个合适的降温函数非常重要,冷却进度表就是指的从某一高温状态T向低温状态冷却时的降温函数。

对物理退火过程和模拟退火算法进行一个简单的对比:

img

二、算法流程

有了上面对该算法基础的理解后,下面对算法的流程进行介绍:

img

可以看出该算法有内外两次循环,内循环用来模拟同一温度下多次的状态转移,也称为抽样稳定准则(为什么称为“抽样稳定”准则,我的理解是等温过程在邻域选择下一次状态转移,一般需要在邻域内采样多次,进行多次状态转移,进而达到这一温度下的稳定状态),外循环作为算法终止准则。
另外其中有3个比较重要的函数,分别为状态产生函数(邻域函数),状态接受函数,退温函数。同时另外还需注意到一个关键部分就是初温、初始解的初始化。下面对上述内容展开介绍:

1、状态产生函数(邻域函数)

设计的出发点:尽可能保证产生的候选解遍布全部解空间。该函数包括两部分的内容:产生候选解的方式、候选解产生的概率分布。前者决定由当前解产生候选解的方式,后者决定在当前解产生的候选解中选择不同状态的概率。

候选解的产生方式由问题的性质决定,通常在当前状态的邻域结构内以一定概率方式产生,而邻域函数和概率方式可以多样化设计,其中概率分布可以是均匀分布、正态分布、指数分布、柯西分布等。

2、状态接受函数

状态接受函数的目的是尽可能接受优化解,状态接受函数一般以概率的方式给出,不同接受函数的差别主要在于接受概率的形式不同,设计状态接受概率,应遵循的原则:

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

下面为标准形式下的状态接受概率,通过将该概率与[0,1]之间的随机数进行比较,从而决定是否接受该候选解:

img

3、初温

实验表明:初温值只要选择充分大,获得高质量解的概率就大!但花费计算时间增加。
因此,初温的选择要足够高。初温的确定应折衷考虑优化质量和优化效率,有以下两种方式:

  • 均匀抽样一组状态,选各状态目标值的方差

  • 利用经验公式
    a.随机产生一组状态,确定两两状态间的最大目标值差,然后依据差值,利用一定的函数确定初温。例如

    img

    b img

4、温度更新函数

温度更新函数即温度的下降方式,用于在外循环中修改温度值,衰减量“以小为宜”

实验表明降温速度越慢,获得高质量解的几率越大,但花费的计算时间将同时增加。因此可以使其温度高时下降的慢些,温度低时下降的快些。主要有以下几种方式:
a.img

Kirkpatrick首先提出对于lambd可选0.95,Johnson,Bonomi及Lutton采用[0.5,0.99].

b. Nahar及Skiscim等人把[0,t0 ]划分成K个小区间,温度更新函数为

img

5、Metropolis抽样稳定准则

用于决定在各温度下产生候选解的数目(Lk ).
具体应与问题规模成比例。实验表明高温时迭代次数越多越好,低温时迭代次数可以适当减少。
在这里引入下马尔可夫链

img

img

马氏链模型:

img

时齐算法的收敛性

img

结论1 时齐模拟退火算法对应的有限状态马氏链存在平稳分布。

结论2 当温度趋于0时,马氏链以概率1收敛到最优状态集,而收敛到非最优状态的概率为0。

实现途径:通过各温度下各状态序列无限长得以实现!

非时齐算法的收敛性:

img

对于非时齐算法

img

SA算法从某一个初始状态开始后,每一步状态转移均是在当前状态的邻域中随机产生新状态,然后以一定概率进行接受的。接受概率仅依赖于新状态和当前状态,并由温度加以控制。因此,SAA对应一个马氏链。

对于时齐算法,若固定每一温度,算法均计算马氏链的变化直至平稳分布,然后下降温度。
对于非时齐算法,若无需各温度下算法均达到平稳分布,但温度需按照一定的速率下降。或称非平稳马氏链算法。

  • 马尔可夫链这部分我目前的理解还不是很深刻,后面有了认识后会加以补充。先记住下面的结论:
    非时齐SAA:每个温度下只产生一个或少量候选解。
    时齐算法——常用Metropolis抽样稳定准则包括:
  • 检验目标函数值的均值是否稳定
    连续若干步目标函数值的变化较小
    按一定的步数抽样

6、算法终止准则

用于决定算法何时结束。
理论上要求温度终值趋于零,通常的做法包括:

  • 设置终止温度的阀值
  • 设置外循环迭代次数(6-50)
  • 算法搜索到的最优值连续若干步保持不变
  • 检验系统熵是否稳定(暂未理解)

小结:
模拟退火算法基本要素和设定方法

img

三、算法示例

1、案例1

在这里插入图片描述

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


# define aim function
def aimFunction(x):
    y = x ** 3 - 60 * x ** 2 - 4 * x + 6
    return y


x = [i / 10 for i in range(1000)]
y = [0 for i in range(1000)]
for i in range(1000):
    y[i] = aimFunction(x[i])

plt.plot(x, y)
T = 1000  # initiate temperature
Tmin = 10  # minimum value of terperature
x = np.random.uniform(low=0, high=100)  # initiate x
k = 50  # times of internal circulation
y = 0  # initiate result
t = 0  # time
while T >= Tmin:
    for i in range(k):
        # calculate y
        y = aimFunction(x)
        # generate a new x in the neighboorhood of x by transform function
        xNew = x + np.random.uniform(low=-0.055, high=0.055) * T
        if (0 <= xNew and xNew <= 100):
            yNew = aimFunction(xNew)
            if yNew - y < 0:
                x = xNew
            else:
                # metropolis principle
                p = math.exp(-(yNew - y) / 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

print(x, aimFunction(x))

img

最小值大概在40左右,通过求导计算得到最小值为40.033

2、案例2

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

img

这是一个旅行商问题。给我方基地编号为0,目标依次编号为1,2,…,100,最后我方基地再重复编号为101(这样便于程序中计算)。距离矩阵 ,其中 表示表示 两点的距离, ,这里 为实对称矩阵。则问题是求一个从点0出发,走遍所有中间点,到达点101的一个最短路径。

img

img

img

img

img

img

img

from numpy import loadtxt,radians,sin,cos,inf,exp
from numpy import array,r_,c_,arange,savetxt
from numpy.lib.scimath import arccos
from numpy.random import shuffle,randint,rand
from matplotlib.pyplot import plot, show, rc
a=loadtxt("Pdata17_1.txt")
x=a[:,::2]. flatten(); y=a[:,1::2]. flatten()
d1=array([[70,40]]); xy=c_[x,y]
xy=r_[d1,xy,d1]; N=xy.shape[0]
t=radians(xy)  #转化为弧度
d=array([[6370*arccos(cos(t[i,0]-t[j,0])*cos(t[i,1])*cos(t[j,1])+ sin(t[i,1])*sin(t[j,1])) for i in range(N)]
   for j in range(N)]).real
savetxt('Pdata17_2.txt',c_[xy,d]) #前两列是102*2的矩阵代表102个点的位置坐标,后面是102*102的矩阵,代表的是距离矩阵
path=arange(N); L=inf
for j in range(1000):   #蒙特卡洛
	path0=arange(1,N-1); shuffle(path0)
	path0=r_[0,path0,N-1]; L0=d[0,path0[1]]  #初始化
	for i in range(1,N-1):L0+=d[path0[i],path0[i+1]]
	if L0<L: path=path0; L=L0
print(path,'\n',L)  

e=0.1**30; M=20000; at=0.999; T=1   #模拟退火
for k in range(M):
	c=randint(1,101,2); c.sort()
	c1=c[0]; c2=c[1]
	df=d[path[c1-1],path[c2]]+d[path[c1],path[c2+1]]-d[path[c1-1],path[c1]]-d[path[c2],path[c2+1]]
  
	if df<0:
    	path=r_[path[0],path[1:c1],path[c2:c1-1:-1],path[c2+1:102]]; L=L+df
	else:
    	if exp(-df/T)>=rand(1):
        	path=r_[path[0],path[1:c1],path[c2:c1-1:-1],path[c2+1:102]]
        	L=L+df
	T=T*at
	if T<e: break54 
print(path,'\n',L)  #输出巡航路径及路径长度
xx=xy[path,0]; yy=xy[path,1]; rc('font',size=16)
plot(xx,yy,'-*'); show()  #画巡航路径

img

————————————————
资料来源:CSDN博主「馋学习的身子」的原创文章
原文链接:https://blog.csdn.net/qq_38048756/article/details/109305769

评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

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

余额充值