模拟退火算法及其改进

第4卷第2期2007年4月

工程地球物理学报

CHIN ESE JO U RN A L O F EN GI NEERIN G G EOP HY SICS

V ol 14,N o 12Apr 1,2007

文章编号:1672)7940(2007)02)0135)06

模拟退火算法及其改进

蒋龙聪,刘江平

(中国地质大学地球物理与空间信息学院,武汉430074)

作者简介:蒋龙聪(1983)),男,硕士研究生,现在主要从事地震数据处理和反演理论方法研究。E -mail:longcja@https://www.wendangku.net/doc/1712623804.html,

刘江平(1957)),男,教授,博士生导师,主要从事地震勘探的科研与教学工作。E -mail:liujp@https://www.wendangku.net/doc/1712623804.html,

摘 要:借鉴遗传算法中的非均匀变异思想,用非均匀变异策略对当前模型扰动产生新的模型,对传统的模

拟退火算法提出了改进,通过多峰值函数数值优化测试结果表明,该算法在高温的时候能够进行大范围的搜索,随着温度的降低,逐渐缩小解的搜索范围,大大加快了收敛速度,证实了该改进算法的有效性和高效性。

关键词:模拟退火算法;非均匀变异;数值最优化;反演中图分类号:P631

文献标识码:A

收稿日期:2006)12)07

Revised Simulated Annealing Algorithm

Jiang Longcong,Liu Jiangping

(I ns titute of Geop hy sics and Geomatics ,China Univer sity of Geosciences ,W uhan 430074,China)

Abstract:Based on the idea o f non -uniform mutation in genetic alg orithm ,w e pr esent a novel rev ised simulated annealing (RSA),w hich used the no n -unifo rm m utatio n to generate a new model from current m odel.T ested by so me numerical functions,RSA can search in the large

ar ea for the solutions in high temperature.With the low er ing of the tem perature,the area of searching the so lutions w ill be g radually reduced and conv erg ence w ill speed up.So the re -sults prove the effectiveness o f RSA.

Key words:simulate annealing;non -uniform m utatio n;numerical o ptimal;inversion

1 引 言

人类对地球内部物理性质(包括速度、密度、电导率、温度等)以及矿产资源分布的了解,大多来自地表地质和地球物理、地球化学资料的反演和解释[1]。反演方法可以分为线性反演和非线性反演两种,线性反演已成为一套科学的反演理论,然而,绝大部分地球物理问题都是非线性的,并且实践表明,线性反演方法有容易陷入局部极值和依赖于初始值等缺点。因此,地球物理学者们不

断的尝试开发非线性反演方法,比如人工神经网

络[2]、小波多尺度反演[3]、模拟退火算法[4]等。

模拟退火算法是近年发展起来的全局最优化算法,其主要优点是;不用求目标函数的偏导数及解大型矩阵方程组,即能找到一个全局最优解,而且易于加入约束条件,编写程序简单。目前此法已开始用于解决非线性地球物理反问题,如波形反演、静校正、叠前偏移速度分析等非线性反演中,并取得了较好的效果。

然而,由于模拟退火法是建立在随机搜寻方法的基础上,要达到一定的精度要求,每一模型参

数的离散点必须足够大。每次迭代必须进行多次目标函数计算,因而在处理实际资料时计算效率不高,影响着它的广泛应用。为了提高模拟退火算法的计算效率,出现了许多改进的方法,如采用依赖温度的Cauchy或似Cauchy分布代替常规模拟退火方法中的高斯分布产生新模型,基于Cauchy分布产生新模型的优点是,高温状态下可进行大范围的搜索,低温状态下只在当前模型附近进行搜索,而且由于似Cauchy分布有一个平坦的/尾巴0,使其易于迅速跳出局部极值区。王山山等[5]将Cauchy分布和Gibbs分布结合起来产生新模型,使得反演过程更加合理,加快了反演过程的收敛,并且提高了算法的抗干扰能力。王凌等[6]详细研究函数优化中基于Cauchy分布的状态发生器SGC(State generator based on Cauchy distribution)和基于Gaussian分布的状态发生器SGG(State generator based o n Gauss-i an distribution)对SA(Simulated annealing)算法性能的影响。对分布机制的研究表明,SGC有利于大范围搜索和脱离极小区域,而SGG较适合于局部搜索。对不同复杂度的典型问题的仿真表明,优化简单单极小问题时SGC的优化效率优于基于SGG,优化复杂多极小或存在平坦区的简单问题时SGC的优化度和鲁棒性均优于SGG,进而利用对尺度参数的/退温0控制,提出了SGC 的改进策略,较大程度上提高了算法的优化度。

纪晨等[7]在模拟退火算法中引入均匀设计方法,Basu等[8]提出用试验方法确定临界温度,算法由稍高于临界温度开始,Press等[9]采用单纯形法与模拟退火算法相结合的综合算法,在不同程度上提高了模拟退火法的计算效率,然而在实际应用中还有待于提高。

随着多维分形(M ultifractals)理论与应用研究日益受到重视,人们展开了关于这一理论的统计学研究。T sallis给出了广义Bo ltzmann-Gibbs 统计理论及相应的Gibbs分布。基于这一理论, Penna提出了新的模拟退火方法,即由广义Gibbs分布给出新的接收概率计算表达式,用于求解推销员问题,表明这种方法在较高的温度就能收敛到全局极值,具有较高的计算效率。张霖斌等以广义Boltzmann-Gibbs统计理论为基础,采用非常快速模拟退火法中依赖于温度的似Cauchy分布产生新的扰动模型,提出了快速模拟退火算法FSA(Fast Simulated Annealing),进一步提高了模拟退火方法的计算效率[10]。

笔者借鉴遗传算法中的非均匀变异思想[11],用非均匀变异策略产生新的扰动模型,对传统的模拟退火算法进行了改进(Rev ised simulated an-nealing,RSA),该算法在高温的时候能够进行大范围的搜索,随着温度的降低,逐渐缩小解的搜索范围,能够大大加快收敛速度,通过对多峰值函数数值优化,测试结果表明该改进算法的有效性和高效性。

2模拟退火原理及其改进

211模拟退火原理

模拟退火算法(Sim ulated Annealing)源于统计物理学,据统计热力学,物体中的每个分子的状态服从Gibbs分布,即:

Q(r i)=

ex p(-

E(r i)

kT)

E M

j=1

exp(-

E(r j)

kT)

(1)

式中:E(r i)为第i个分子的能量函数;r i为第i个分子所处的状态;k为波耳兹曼常数,T表示温度;Q(r i)为第i个分子的概率密度,为了方便起见令k=1。

模拟退火算法最早的思想是由M etro po lis 在1953年提出的,由Kir kpatrick于1983年成功地应用在组合优化问题中,目前已经应用到各门学科中以解决非线性系统中的优化问题。SA法是局部搜索算法的扩展,它不同于局部搜索之处是以一定的概率选择领域中的最优值状态。理论上已经证明,它是一个全局最优算法,并且以概率1接近最优值。

算法源于对实际固体退火过程的模拟,即先将固体加温至充分高,再逐渐冷却。加温时,固体内部粒子变为无序状态,内能增大;而逐渐降温时,粒子趋于有序,在每个温度都达到平衡态,最后在常温时达到基态,内能减为最小。因此,算法实际上是将优化问题类比为退火过程中能量的最低状态,也就是温度达到最低点时,概率分布中具有最大概率的状态。

模拟退火算法的具体步骤如下:

136工程地球物理学报(Chinese Jour nal of Engineer ing Geo phy sics)第4卷

1)给定模型每一个参数变化范围,在这个范围内随机选择一个初始模型m0,并计算相应的目标函数值E(m0);

2)对当前模型进行扰动产生一个新模型m,计算相应的目标函数值E(m),得到

$E=E(m)-E(m0)

3)若$E<0,则新模型m被接受;若$E>0,则新模型m按概率P=ex p(-$E/T)进行接受, T为温度。当模型被接受时,置m0=m;

4)在温度T下,重复一定次数的扰动和接受过程,即重复步骤2)、3);

5)缓慢降低温度T;

6)重复步骤2)、5),直至收敛条件满足为止。

为了避免最好的解在优化过程中被忽略掉,可以稍做改进,即在整个搜索过程中随时记下最好的值,因为该法的特点决定了最后的优化值并不一定就是最优的值,而只能是较优值。

212算法的改进

21211模型扰动

模拟退火算法中新模型的产生是对当前模型进行扰动得到的,这个扰动是用随机函数控制的。Ing her于1989年提出了速度非常快的模拟退火算法V FSA(Very Fast Simulated annealing),在该算法中,采用了依赖于温度的似Cauchy分布产生新的模型[5],即

x c i=x i+y i(x i max-x i min)

y i=T k sg n(N-0.5)

@[(1+1/T k)|2N-1|-1](2)式中,x i为当前模型参数,x c i为扰动后的模型参数,N为(0,1)区间上的随机数,x i min,x i max为x i的取值范围,sg n为符号函数,y i称之为扰动因子。

笔者分析了快速模拟退火算法的改进扰动策略,提出了一种强化局部搜索能力的算法,该算法借鉴了遗传算法中的非均匀变异思想[11],用非均匀变异策略对当前模型参数扰动产生新的模型参数,即:

x c i=x i+y i(x i max-x i min)

y ic=r(1-t/N)K sgn(r-0.5)(3)其中,r为(0,1)之间的随机数,t为当前温度,N 为最大迭代次数,N与最大温度、最低温度有关, K是确定非均匀性程度的常数,也称之为形状因子。(3)式表明,在高温时期,搜索的范围比较广,随着温度的降低,逐渐缩小搜索范围,为此,分析了当K=015,1,2,5,10时相应的y ic值,指导实际计算中如何取值,并且跟(2)式y i做了对比分析(图1)。

从图1可以得出,随着迭代次数的增加(即温度的降低),y ic的值逐渐减小,也即缩小了搜索范围,加强了局部搜索的能力,常规的模拟退火算法和快速模拟退火算法都是随机进行大范围的搜索,而且,当形状因子K越大,越容易进行局部搜索,因此,在实际计算中,要合理地控制好K值,K 太大容易陷入局部极值中,M ichalew icz的建议和笔者的数值模拟表明是K取[2,5]是比较好的。随着温度的降低,新接收的模型越来越接近真实的求解模型,因此,加强局部搜索更有利于提高算法的收敛性能。

21212接收概率

对于接受概率,快速模拟退火算法采用广义Gibbs分布产生接收概率P,即

P=[1-(1-h)$E/T]1/[1-h](4)式中$E表示能量差,T表示温度,h为一实数。按照(4)式计算得到的概率进行判断是否接收新模型,该式只考虑了能量的绝对变化,然而,物体由熔化状态逐渐冷却至结晶状态所出现的相变规则,必须考虑能量的相对变化,为此顾汉明等[12]对此进行了修正,即:

P=[1-(h$E/T)(E A1/E B)]1/h(5)式中,E为新模型的目标函数,E1为当前模型的目标函数,A,B为非负实数。

21213降温方式

对于退火方案,实践表明,T取指数变化的模式比较切合退火的实质。Ingber给出的快速退火算法中的退火方式为

T(K)=T0exp(-CK1/N)(6)式中:T0为初始温度,K为迭代次数,C为给定的常数,N为待反演的参数的个数。该式子也可以改成

T(K)=T0A k1/N(7)通常选择017[A[1,在本文中,用1代替1/N,A 取为0199。

在改进的算法中RSA,采用(3)式对当前模型扰动产生新的模型,新模型按照(5)进行接收,按照(7)式缓慢降温。用M atlab语言编写了本文中的算法,用于数值优化试验。

137

第2期蒋龙聪等:模拟退火算法及其改进

图1扰动因子随迭代次数变化的曲线

(a)为常规模拟退火的扰动;(b)、(c)、(d)、(e)、(f)为K=0.5,1,2,5,10

Fig.1P er turbatio n facto r via iterat ive number

3数值优化

为了验证改进的算法的有效性,选取三个多峰值函数作为测试对象,计算该函数的最小值。用改进的模拟退火算法和普通模拟退火算法对三个函数各运行20次,收敛曲线为采用局部策略之后的局部最优收敛曲线。3.1函数1:Shuber函数

f1(x1,x2)=F2k=1E5i=1i cos[(i+1)x k+i]

(8)

x k I(-10,10)该函数有全局最小值-1861731,用改进的算法(RSA)求取该函数的最小值,运行20次,每次都能找到全局最小值,而常规的模拟退火算法12次找不到全局最小值,8次

138工程地球物理学报(Chinese Jour nal of Engineer ing Geo phy sics)第4卷

找到全局最小值,常规模拟退火易陷入局部极小值。初始温度为10000,最小温度为010001,马尔科夫链长度为3,K 为5,初始值x 1、x 2同时为-10。采用了记忆方法,记忆了搜索过程中每一温度下的局部最优解,收敛过程如图2-图4

。

3.2 函数2:De Jong 函数

f 2(x 1,x 2)=100(x 21-x 2)2+(1-x 1)2

-2.048[x i [2.048(9)函数在(x 1,x 2)=(1,1)处有一个全局最小值0。用两种算法计算该函数最小值,各运行20次,改进后的模拟退火算法每次都能找到全局最小值,而普通的模拟退火算法13次找到全局最小值。初始温度为10000,最小温度为010001,马尔科夫链长度为3,K 为5,初始值x 1,x 2同时为-21048。3.3 函数3:Easom 函数

f 3(x 1,x 2)=-cos (x 1)cos (x 2)e -(x 1-P )2

-(x 1-P

)

2

-100[x i [100(10)函数在(x 1,x 2)=(P ,P )处有一个全局最小值-1.用两种算法计算该函数的最小值,各运行20次,改进后的模拟退火算法18次找到全局最小值,两

次陷入函数值0中,普通的模拟退火算法12次找到全局最小值,8次陷入函数值0中。初始温度为10000,最小温度为010001,马尔科夫链长度为3,K 为3,初始值x 1,x 2同时为-100。

随着温度的下降,函数值逐渐靠近全局最优解,此时宜采用局部搜索方法,因此改进后的模拟退火算法能够很快搜索到全局最优解,而常规的模拟退火算法由于是对解空间的解随机扰动,很容易在很长时间内,搜索不到比当前局部最优解更好的解,对三个函数的测试表明了该点。

4 结 语

1)改进的模拟退火算法相对于常规的模拟退火算法,加强了局部搜索功能,更有利于提高算法的收敛速度,数值优化结果证实改进算法的有效性和高效性。

2)对于改进的算法,K 的选择至关重要,K 太大容易陷入局部最小值,K 太小则降低了收敛速度,因此要通过测试再决定取K 的值。

3)从本质上讲,模拟退火算法属于蒙特卡罗方法,它用随机的方法代替了线性反演中的步长,相对于线性反演来讲,模拟退火算法的执行效率还是远远不及线性反演,因此,如何考虑将线性反演与非线性反演结合起来,相互利用对方的优势,并且用于实际反演中将是接下来重点研究的对象。参考文献:

[1]王家映.地球物理反演理论[M ].北京:高等教育出版

139

第2期 蒋龙聪等:模拟退火算法及其改进

社,20001

[2]徐海浪,吴小平.电阻率二维神经网络反演[J].地球物

理学报,2006,49(2):584~5891

[3]徐义贤,王家映.大地电磁多尺度反演[J].地球物理学

报,1998,41(5):704~7111

[4]姚姚.地球物理非线性反演模拟退火法的改进[J].地

球物理学报,1995,38(5):643~6501

[5]王山山,李灿平,李青仁,等.快速模拟退火地震反演

[J].地球物理学报,1995,38(增刊1):123~1341 [6]王凌,郑大钟.基于Cauchy和G aussian分布状态发

生器的模拟退火算法[J].清华大学学报(自然科学版),2000,40(9):109~1121

[7]纪晨,姚振兴.用于地球物理反演的均匀设计优化算

法[J].地球物理学报,1996,39(2):233~242.

[8]Basu A,F razer L.Rapid deter mination o f the cr itical

temperature in simulated annealing inv ersio n[J].Sci-ence,1990,249(21):1409~1412.

[9]Pr ess W H,V et terling W T,T eukolsky S A,et al.

Simulated annealing o ptimization o ver continus space [J].Comp uter s in Phy sics,1991(5):426~429. [10]张霖斌,姚振兴,纪晨,等.快速模拟退火算法及应用

[J].石油地球物理勘探,1997,32(5):654~660. [11]M icha lew icz Z.演化程序)遗传算法和数据编码的结

合[M].北京:科学出版社,2000.

[12]顾汉明,江涛,王家映.改进的模拟退火方法进行

A V O岩性参数反演[J].地球科学,1999,24(4):418

~422.

[13]钟守楠.遗传算法的收敛性与编码[J].武汉水利电力

大学学报,2000,33(1):108~112.

[14]张宏兵,尚作萍,谭胜章,等.波阻抗反演的快速模拟

退火算法[J].河海大学学报,2005,33(4):434~437.

[15]M D M art?nez,X L ana,J O lar te etc.I nv ersion o f

R ayleigh w ave phase and g ro up v elo cities by simula-ted annealing[J].Phys ics of the Ear th and Planeta-r y I nter ior s,2000,122(1):3~17.

[16]王薇,曾光明,何理.用模拟退火算法估计水质模型

参数[J].水利学报,2004,11(6):61~67.

[17]金炳尧.最新优化技术中的若干新技术[J].科技通

报,2000,16(2):119~1241

[18]邢文训,谢金星.现代优化计算方法[M].北京:清华

大学出版社,19991

[19]K irkpat rick S,Gelatt Jr C D,Vecchi M P.O pt im-i

zation by simulat ed annealing[J].Science,1983,220

(13):671~6801

[20]M antaw y A H,A bde-l M ag id,Yo ussef L,et al.A

simulated annealing a lg orithm fo r unit co mmitment [J].I EEE T r ansactions on P ow er Sy s tems,1998, 13(1):197~2051

[21]杨东,于桂荣.基于模拟退火算法预测储层参数[J].

沈阳航空工业学院学报.2005,22(1):76~77. [22]唐立山,谢云,尤矢勇,等.非数值并行算法[M].北

京:科学出版社,1994.

140工程地球物理学报(Chinese Jour nal of Engineer ing Geo phy sics)第4卷