模拟退火算法

模拟退火算法

(SimulatedAnnealing)1983，Kirkpatrick成功应用于组合优化问题起源：1953，Metropolis，固体退火过程SA是基于MonteCarlo迭代求解策略的一种随机寻优算法，理论上是一个全局最优算法。其出发点是基于物理中固体物质的退火过程与一般组合优化问题之间的相似性。统计力学表明材料中粒子的不同结构对应于粒子的不同能量水平。在高温条件下，粒子的能量较高，可以自由运动和重新排列。在低温条件下，粒子能量较低。如果从高温开始，非常缓慢地降温（这个过程被称为退火），粒子就可以在每个温度下达到热平衡。当系统完全被冷却时，最终形成处于低能状态的晶体。统计力学的研究表明，在温度T，进行了充分转换之后，材料将达到热平衡。此时分子停留在状态r的概率满足Boltzmann分布：其中，E(r)为状态r的能量，为玻尔兹曼常数，为分子能量的一个随机变量，D表示状态空间集合.在同一个温度，分子停留在能量小状态的概率比停留在能量大状态的概率大。当温度相当高时，每个状态的概率基本相同，接近平均值当温度越低时，能量越低的状态的概率值越高。当温度趋于0时，分子停留在最低能量状态的概率趋于1。类比：组合优化问题（可行）解最优解费用函数金属物体状态能量最低的状态能量模拟退火算法基本流程：Step1：随机产生一个初始解，令，并计算目标函数值若则否则若时Step3：若满足停止条件，终止计算，输出最优解；否则，回到step2.Step2：若在该温度达到内循环停止条件，则到step3；否则，从当前最优解的邻域中随机选一计算新的目标函数值并计算目标函数值的增量模拟退火算法的数学模型：马尔可夫链（可达性，渐近不依赖起点，分布稳定性，收敛到最优解）如果温度下降十分缓慢，而在每个温度都有足够多次的状态转移，使之在每一个温度下达到热平衡，则全局最优解将以概率1被找到。因此可以说模拟退火算法能找到全局最优解。模拟退火算法中，包含一个内循环和一个外循环，内循环表示在同一个温度时，在一些状态随机搜索。外循环主要包括温度下降变化，迭代步数的增加和停止条件。模拟退火算法实现的技术问题1）解的形式和邻域结构解的表现形式直接决定于邻域的构造例：的解可采用0-1编码TSP问题可采用n个城市的一个排序表示问题的一个解，可通过城市间不同位置交换构造邻域例：同样的解的表达式和邻域结构用于背包问题和车间作业问题时，由于背包问题有包容积约束的限制和车间作业的死锁现象，无法保证每个邻域中邻居都是可行解。处理这个问题有两类方法：一是用罚函数将不可行解可行化；另一个方法是研究解和邻域结构，这一方法要求对问题有相当深入的了解，难度较大。2）状态转移概率（接受概率）p状态转移概率是指从一个状态xold（一个可行解）向另一个状态xnew（另一个可行解）的转移概率通俗的理解是接受一个新解为当前解的概率它与当前的温度参数T有关，随温度下降而减小一般采用Metropolis准则3）初始温度T0实际计算中，可以选K=10，20，100…等试验值K充分大的数实验表明，初温越大，获得高质量解的几率越大，但花费的计算时间将增加；初温过低则使算法过早陷入局部最优点。因此，初温的确定应折中考虑优化质量和优化效果。4）冷却进度T(t)冷却进度表是指从某一高温状态T0向低温状态冷却时的降温管理表。经典模拟退火算法的降温方式为：快速模拟退火算法的降温方式为：常用的两种简单下降方式：其中其中K为算法温度下降的总次数需综合考虑解的性能和算法速度5）内循环终止准则或称Metropolis抽样稳定准则，用于决定在各温度下产生候选解的数目。常用的方法：固定长度：在每一温度迭代相同的步数，步数的选取同问题的大小有关，通常采用与邻域大小直接相关的规则。由接受和拒绝的比率来控制迭代步数：随着温度的下降，将同一温度的迭代步数增加。6）外循环终止准则即算法终止准则，常用的包括：零度法：设置终止温度的阀值：小正数循环总数控制法：如设置外循环迭代次数基于不改进规则的控制法：算法搜索到的最优值连续若干步保持不变2）优点：高效、健壮、通用、灵活模拟退火算法实验性能分析1）模拟退火法与局部搜索算法的差异3）不足：返回一个高质量近似解的时间花费较多，当问题规模不可避免增大时，难于承受的运行时间将使算法丧失可行性。4）模拟退火算法理论上具有全局最优性，但实际应用中的模拟退火算法是一个启发式算法，它有诸多的参数需要调整。解决这个矛盾的方法主要通过大量的数值模拟计算，从中选择比较好的参数配置。模拟退火算法的matlab实现利用模拟退火算法求函数极小：无约束或有界约束x=simulannealbnd(fun,x0)x=simulannealbnd(fun,x0,lb,ub)x=simulannealbnd(fun,x0,lb,ub,options)[x,fval]=simulannealbnd(...)[x,fval,exitflag]=simulannealbnd(...)[x,fval,exitflag,output]=simulannealbnd(...)dejong5fcnx0=[00];[x,fval]=simulannealbnd(@dejong5fcn,x0)Optimizationterminated:changeinbestfunctionvaluelessthanoptions.TolFun.x=0.0216-31.9955fval=2.9821Optimizationterminated:changeinbestfunctionvaluelessthanoptions.TolFun.x=-15.9669-31.9749fval=1.9920exitflag=1output=iterations:1608funccount:1621message:[1x80char]randstate:[625x1uint32]randnstate:[2x1double]problemtype:'unconstrained'temperature:[2x1double]totaltime:0.8268[x,fval,exitflag,output]=simulannealbnd(@dejong5fcn,[0,0])x0=[00];lb=[-64-64];ub=[6464];[x,fval]=simulannealbnd(@dejong5fcn,x0,lb,ub)Optimizationterminated:changeinbestfunctionvaluelessthanoptions.TolFun.x=-32.0169-31.9879fval=0.9980fun=@(x)3*sin(x(1))+exp(x(2));x=simulannealbnd(fun,[1;1],[00])Optimizationterminated:changeinbestfunctionvaluelessthanoptions.TolFun.x=896.92340.0000options=saoptimset('PlotFcns',{@saplotbestf,@saplottemperature,@saplotf,@saplotstopping});simulannealbnd(@dejong5fcn,x0,lb,ub,options);options=saoptimset('InitialTemperature',[30050]);options=saoptimset(options,'TemperatureFcn',@temperaturefast);options=saoptimset(options,'ReannealInterval',50);options=saoptimset(options,'Display','iter','DisplayInterval',400);options=saoptimset(options,'TolFun',1e-5);1）加载数据表129个城市的坐标城市序号12345678910X坐标1150.0630.040.0750.0750.01030.01650.01490.0790.0710.0Y坐标1760.01660.02090.01100.02030.02070.0650.01630.02260.01310.0城市序号11121314151617181920X坐标840.01170.0970.0510.0750.01280.0230.0460.01040.0590.0Y坐标550.02300.01340.0700.0900.01200.0590.0860.0950.01390.0城市序号212223242526272829X坐标

830.0

490.01840.01260.01280.0490.01460.01260.0360.0Y坐标1770.0500.01240.01500.0790.02130.01420.01910.01980.0模拟退火算法对TSP的应用inputcities=[1150.0630.040.0750.0750.01030.01650.01490.0...790.0710.0840.01170.0970.0510.0750.01280.0230.0460.0...1040.0590.0830.0490.01840.01260.01280.0490.01460.0...1260.0360.0;1760.01660.02090.01100.02030.02070.0650.0...1630.02260.01310.0550.02300.01340.0700.0900.01200.0...590.0860.0950.01390.01770.0500.01240.01500.0790.0...2130.01420.01910.01980.0];2）计算总距离的函数该函数根据给定路径计算该路径对应总路程。functiond=distance(inputcities)%d=distance(inputcities)计算TSP问题中n个城市间某路径%的总路程%输入参数为2*n数组，其中n为城市数目，列为对应城市坐标

d=0;forn=1:length(inputcities)ifn==length(inputcities)d=d+norm(inputcities(:,n)-inputcities(:,1));elsed=d+norm(inputcities(:,n)-inputcities(:,n+1));endend3）绘制路线的函数functionf=plotcities(inputcities)%plotcities(inputcities)从inputcities参数中绘制城市的位置%输入参数为2*n数组，其中n为城市数目，列为对应城市坐标temp_1=plot(inputcities(1,:),inputcities(2,:),’b*’);set(temp_1,’erasemode’,’none’);temp_2=line(inputcities(1,:),inputcities(2,:),’Marker’,’*’);set(temp_2,’color’,’r’);x=[inputcities(1,1)inputcities(1,length(inputcities))];y=[inputcities(2,1)inputcities(2,length(inputcities))];x1=10*round(max(inputcities(1,:))/10);y1=10*round(max(inputcities(2,:))/10);ifx1==0x1=1;endify1==0y1=1;endaxis([0x10y1]);temp_3=line(x,y);set(temp_3,’color’,’r’);dist=distance(inputcities);distance_print=sprintf(...’Theroundtriplengthfor%dcitiesis%4.6funits’...,length(inputcities),dist);text(x1/15,1.05*y1,distance_print,’fontweight’,’bold’);drawnow;4）交换城市的函数%s=swapcities(inputcities,n)n个城市随机的交换返回一组m个城市s=inputcities;fori=1:ncity_1=round(length(inputcities)*rand(1));ifcity_1<1city_1=1;endcity_2=round(length(inputcities)*rand(1));ifcity_2<1city_2=1;endtemp=s(:,city_1);s(:,city_1)=s(:,city_2);s(:,city_2)=temp;end5）执行模拟退火的函数functions=simulatedannealing(inputcities,initial_temperature,...cooling_rate,threshold,numberofcitiestoswap)%s=simulatedannealing(inputcities,initial_temperature,cooling_rate,…)%通过n个城市的TSP问题的最优解返回一个新的城市构型。globaliterations;temperature=initial_temperature;initial_cities_to_swap=numberofcitiestoswap;iterations=1;complete_temperature_iterations=0;whileiterations<thresholdprevious_distance=distance(inputcities);temp_cities=swapcities(inputcities,numberofcitiestoswap);current_distance=distance(temp_cities);diff=abs(current_distance-previous_distance);ifcurrent_distance<previous_distanceinputcities=temp_cities;plotcities(inputcities);ifcomplete_temperature_iterations>=10temperature=cooling_rate*temperature;complete_temperature_iterations=0;endnumberofcitiestoswap=round(numberofcitiestoswap*exp(-diff/(iterations*temperature)));ifnumberofcitiestoswap==0numberofcitiestoswap=1;enditerations=iterations+1