模拟退火算法学习及实验分析

1、模拟退火算法学习及试验分析清华大学计算机系 李军 2007-4大纲n1. 介绍n2. Six-hump camel back function 试验n3. Schwefels function 试验对比n4. 试验总结n5. 结论与未来工作n6. 参考1.1 优化问题介绍n描述: Find the values of a vector that minimize a scalar valued loss function L(). : the domain of allowable values for a vector 注注: loss function 也被称为也被称为: performa

2、nce measure, cost function, objective function,fitness function, or criterion etc.1.2 模拟退火算法介绍n用于解决优化问题的一种启发式算法,理论上是一个全局最优算法.n以一定概率跳出局部极值区域从而增大了找到全局极值的概率.n伪码描述:Simulated-Annealing()Create initial solution Srepeatfor i=1 to iteration-length doGenerate a random transition from S to SiIf ( C(S) random0

3、,1) ) thenS=SiReduce Temperature tuntil ( no change in C(S) )C(S): Cost or Loss function of Solution S.2. Six-hump camel back function 试验nThe 2-D Six-hump camel back function is a global optimization test function. n global minimum: f(x1,x2)=-1.0316; (x1,x2)=(-0.0898,0.7126), (0.0898,-0.7126).n注: 在简

4、单问题上成立的结论才有可能推广到更复杂的问题上.22222121412121)44()31 . 24(),(xxxxxxxxxf 331x221x2.1 函数值分布图2.2 函数值分布底部区域局部 2 .3 与旅行商问题对比优化问题术语Six Camel function problemTSP(Traveling salesman problem)损失(代价)函数值f(x1,x2)路径长度初始解(x1,x2 ) 坐标初始城市排列邻域结构任选一维加上一个符合正态分布的随机数e.g., 2-opt(2个元素交换) 或者 inversion, translation, and switching (

5、部分反转,移动与交换的混合)等全局最优解f(x1,x2)最小路径最短局部最优解f(x1,x2)较小路径相对较短邻域大小增大或减少随机数的大小相对上一个解变动较大. (e.g., 增大部分反转的比例,减少移动和交换的比例)2.4 主要试验参数设置n CoolSched: (0.8*T) %温度下降速率0.8n Generator: 生成邻域: 从 x,y中随机地选择一个再加上一个随机数R , R = randn/100; (rand符合标准正态分布N(0,1)的伪随机数, 意味着 随机变量落入-1,1内的概率是68.26%, 落入-2,2内的概率是 95.44% 落入-3,3内的概率是 99.7

6、2%, 除以100以后也就是大概范围在e-3量级的小数,也就是大约以95%的概率处于-0.02,0.02)nInitTemp: 1 %起始温度nMaxTries: 300 %同一温度下的最大迭代次数nStopTemp: 1e-8 %终止温度n 2.5 初始解的位置对最终解的影响 (R=randn/100)xySix hump camel back function -3-2-10123-1.5-1-0.500.511.5-20246810-4-3-2-101234-3-2-10123effect of initial position -1-0.500.511.52x1,x2分别在区间-3,3

7、, -2,2 步长0.5进行grid 式的初始值设置,然后用模拟退火求最小值的结果( 每一点对应运行一次模拟退火算法之后得到的解)Six-hump camel back function 极值分布的等高线图n 可以看出, 在当前参数设置情况下,初始解的位置与局部极值的区域基本是一一对应的.也就是从初始解的位置出发,通过邻域搜索,往往落入最近的局部极值区域.2.6 R=randn/100的3维效果图-4-2024-202-1-0.500.511.52 -1-0.500.511.522.7 R=randn/10 与 rand/1 的3维效果图-4-2024-202-1-0.500.511.52 -

8、1-0.500.511.52-4-2024-202-1-0.500.511.52 -1.0316-1.0315-1.0314-1.0313-1.0312-1.0311-1.031-1.0309-1.0308-1.03072.8 R=randn*10 与 rand*40 的3维效果图-4-2024-202-1-0.500.511.52 -1.03-1.02-1.01-1-0.99-0.98-0.97-0.96-0.95-4-2024-202-1-0.500.511.52 -1-0.9-0.8-0.7-0.6-0.5-0.42.9 R=randn/100,randn/10,randn/1 的数据比

9、较最终解平均值最终解标准差sd总计迭代次数平均值总迭代次数标准差sdRandn/100 0.0048 1.2351 6731.6 627.6Randn/10-0.6502 0.5410 10705 1082.7 Randn/1-1.0315 0.000111412 1515.3Randn*10-1.02310.0122 11331.71707.3Randn*40-0.9689 0.1001 11703.21422.2n可见,随着邻域的范围的增大, 跳出局部极小区域,最终进入全局极小区域的概率越来越大, 但是代价是总的迭代次数增加.n但是随着邻域范围的增大,会出现所谓的在极值附近来回”振荡”而无

10、法落入极值点的现象. 可以推测,随着邻域范围的进一步增大及其随机特性,模拟退火算法将退化到随机寻找极值并进行记录极值的算法.n注: 当 randn*10, rand*40的时候超出我们限定的搜索范围-3,3,-2,2的概率增大很多,与前3个试验在某种程度上不具备可比性, 尽管如此,仍然具有启发性的意义.2.10 更改降温速率后的运行结果(改为0.95)-4-2024-202-1-0.500.511.52 -1-0.500.511.52-4-2024-202-1-0.500.511.52 -1-0.500.511.52-4-2024-202-1-0.500.511.52 -1.0316-1.03

11、14-1.0312-1.031-1.0308-1.0306-1.0304-1.0302-1.03-1.0298Mean_fvaleStd_fvaluetotalStd_totalRandn/100 0.0922 1.2630 24843 2378.2Randn/10-0.8224 0.3579 39835 2431.3 Randn/1-1.0314 0.000272 40637 3986.5 n可见,随着邻域的范围的增大, 效果与前一个试验类似,区别在于 由于速率放缓, 经历了更多的迭代次数才达到最终解.n而且从中间的图可以看出，进入全局极值的点的初始位置分布较散,这是因为随着迭代次数的增多,

12、以及邻域结构的随机向外伸展性质,由初始位置导致陷入临近局部极值的可能性降低,进入全局极值区域的可能性增大.2.11 其他参数尝试n起始温度(提高到1000)n每个温度时的最大迭代次数(提高到3000)n限于时间,更多参数和变化形式没有进行尝试.n得到的试验结论与前面基本相同,只是在初始解的临近位置周围略有微小的变化.n所以, 在一个合适的参数设置情况下,对寻找最优解起主要作用的重要因素有2个: q初始解的起始位置q邻域解的构造结构2.12 与Native Brute Force Search比较Mean_fvalueMean_total_iterationRandn/1000.00486731

13、.6 Randn/10-0.650210705 Rand/1-1.031511412 nx1,x2以 s的步长在-3,3,-2,2的区间内进行 Brute Force Search, 判断搜索过程中的最小值并记录下来.n在解的搜索空间范围较小时,暴力搜索的优势非常明显,只需很少的迭代次数就能达到与模拟退火相当的程度,但是随着解空间以几何倍数增大,需要的迭代次数也迅速增大.( 本试验中达到同样的精度是模拟退火的21倍.)FvalueTotal_iterationS=1035S=0.5-0.75117S=0.01-1.0315241001降温速率0.8时模拟退火试验结果 暴力搜索的试验结果3. 2

14、-d Schwefels function 试验对比nFunction defination:) |sin()(1iniixxxf 500500 ix50005000 ix3.1 初始解位置(step=40)与最终解值分布图-5000500-5000500-1200-1000-800-600-400-2000 -1100-1000-900-800-700-600-500-400-300-200-100n 观察到了在 Six hump camle function 中同样的试验效果. 区别在于由于解空间的较多个极值区域的存在,迭代次数变化不像只有6个极值区域时为达到全局极值而明显地增加了总的迭代

15、次数.n注1: 由于定义域范围增大,解空间增大,所以将邻域范围设大,从/1到*100.n注2: 由于邻域范围放大之后,处于自定义域边界的初始解生成的第一次邻域解以95%的概率落在-700,700范围之内, 也就得到500,500之外函数极值,所以第3个图中的最小值明显比第1,2个图中的函数值要小,这说明在原定义域外存在更小的极值点n注3: 因为画图的原因,3个图的纵轴的坐标是不一致的, 数值依次降低.也就是跳出了原来的局部极小区域. Randn/1, T=0.8T Randn*40 , T=0.8T Randn*100, T=0.8T-5000500-5000500-3500-3000-250

16、0-2000-1500 -3200-3000-2800-2600-2400-2200-2000-1800mean_fvalue fvalue_std mean_total_itertotal_std Randn/1-525.49 228.411444 881.4 Randn*40-646.67 160.312169 1569.5 Randn*100-2330.90 287.5 11789 1760.5 -5000500-5000500-1400-1200-1000-800-600-400-200 -1200-1100-1000-900-800-700-600-500-4003.2 与Nativ

17、e Brute Force Search比较Mean_fvalueMean_total_iterationRand/1-525.49 11,444 Rand*40-646.67 12,169 Rand*100-2330.90 11,789 nx1,x2以 s的步长在-500,500,-500,500的区间内进行 Brute Force Search, 判断搜索过程中的最小值并记录下来.n同上一个试验结论相同, 在以非常粗的粒度进行解的搜索时,暴力搜索的优势非常明显,只需很少的迭代次数就能达到与模拟退火相当的程度,但是随着解空间以几何倍数增大,需要的迭代次数也迅速增大.( 本试验中达到接近的精度

18、所需迭代次数是模拟退火的573倍.), 进而计算时间也大大增加.FvalueTotal_iterationS=100-730.35121S=40-837.73 676S=1 -500,500S=1 -700,700S=1 -800,800S=1 -1300,1300-837.97-1357.81430.1-2593.11,002,001 1,962,8012,563,2016,765,201 降温速率0.8时模拟退火试验平均结果 暴力搜索的试验结果3.3 3-d Schwefels function-5000500-5000500-5000500 x13-d schwefel function

19、 values distributionx2 x3-1600-1400-1200-1000-800-600nx1,x2,x3以 80的步长在-500,500,-500,500,-500,500的区间内进行模拟退火求极小值, randn*40, 降温速率0.8.n由图中可以看出,较大的解多集中在中间的球形区域,外围的值则要小一些, 这是由于变量以0为中心的区域是相对于外围空间的局部极小值区域,所以从中间区域出发通过模拟退火找到的解也相对于外围要大一些. 得到的结论仍同前面的试验结果相同.n注: 这个图同前几个图不同,Z轴代表变量x3, 颜色的变化代表找到的解的值的大小.4. 试验总结n在2个to

20、y problem 上得到的结论:q在其他参数合适的情况下,初始解的位置与邻域结构是2个重要因素.q初始解的位置与邻域结构的设计之间对于找到最优解的问题而言存在依赖关系.q随着邻域的范围的增大, 跳出局部极小区域,最终进入全局极小区域的概率越来越大, 同时也可能会产生振荡,无法落入极值区域.q模拟退火本质上也是一种 暴力搜索,只不过是利用随机的性质和局部极值区域连续(也取决于求解问题中邻域的结构)的特性避免了大量计算.进而在较短时间内从某种概率上逼近局部最优值和全局最优值.4.2 试验总结-500-400-300-200-1000100200300400500-500-400-300-200-

21、1000100200300400500 x1function value1-d schwefel function (x2=0)q问题本身解空间的性质也会影响求解.-2.5-2-1.5-1-0.500.511.522.53.9583.95853.9593.95953.963.96053.9613.96153.9623.96253.963x 104x1function value1-d Six-hump camel back function (x2=10)-2.5-2-1.5-1-0.500.511.522.5-2-10123456x1function value1-d Six-hump ca

22、mel back function (x2=10)q如果全局最优解不在连续的空间上,(isolated point)那么将非常难以找到.5.1 结论与未来工作1n本试验结论能否推广到更高维空间的问题上,对于不同性质的问题是否有同样的结论还有待更多的试验及数据分析. 例如更改邻域结构,使用tsp问题来分析,使用更高维的实际问题等.n现实问题是复杂的,不同问题的解空间的性质也是仍需探索的.n对问题本身建模,以及构造合适的邻域结构,是解决问题的关键因素.再加上某种程度的初始解的grid式的粗遍历并应用模拟退火或其他类似算法, 或许可以解决很多问题.5.2 结论与未来工作2n对工程实践的启发性意义q先

23、进行粗粒度的搜索, 找到一个较优的解,然后在这个解的周围,利用设计合适的邻域结构,进而找到比这个解更优的解.n理论中通过迭代无穷次达到平稳分布,其实也就是是相当于遍历了所有的解空间.n某种程度的”brute force”(less bugs,no human miss)是解决当前组合优化问题的一种有效途径.qE.g., 1997年深蓝(或许是当时最快的计算机)打败世界国际象棋冠军是暴力搜索的结果. -许丰雄6. 参考n1. matlab code url: http:/ Minimizes a function url: http:/ Steven Skiena The Algorithm Design Manual n4. James C. Spall Introduction to Stochastic Search and Optimizationn5. 刑文训刑文训 谢金星谢金星 现代优化计算方法现代优化计算方法附录1. TSP的邻域的一种构造方法示例nThe Inversion, translation, and switching operations are selected randomly with probabiliti