单钻头退火算法matlab编程clearclca=0.999;%温度衰减函数的参数t0=97;tf=3;t=t0;Markov_length=2800;%Markov链长度coordinates=[];coordinates(:,1)=[];amount=size(coordinates,1);%城市的数目%通过向量化的方法计算距离矩阵dist_matrix=zeros(amount,amount);coor_x_tmp1=coordinates(:,1)*ones(1,amount);coor_x_tmp2=coor_x_tmp1';coor_y_tmp1=coordinates(:,2)*ones(1,amount);coor_y_tmp2=coor_y_tmp1';dist_matrix=sqrt((coor_x_tmp1-coor_x_tmp2).^2+...(coor_y_tmp1-coor_y_tmp2).^2);sol_new=1:amount;%产生初始解%sol_new是每次产生的新解;sol_current是当前解;sol_best是冷却中的最好解;E_current=inf;E_best=inf;%E_current是当前解对应的回路距离;%E_new是新解的回路距离;%E_best是最优解的sol_current=sol_new;sol_best=sol_new;p=1;whilet=tfforr=1:Markov_length%Markov链长度%产生随机扰动if(rand0.5)%随机决定是进行两交换还是三交换%两交换ind1=0;ind2=0;while(ind1==ind2)ind1=ceil(rand.*amount);ind2=ceil(rand.*amount);endtmp1=sol_new(ind1);sol_new(ind1)=sol_new(ind2);sol_new(ind2)=tmp1;else%三交换ind1=0;ind2=0;ind3=0;while(ind1==ind2)||(ind1==ind3)...||(ind2==ind3)||(abs(ind1-ind2)==1)ind1=ceil(rand.*amount);ind2=ceil(rand.*amount);ind3=ceil(rand.*amount);endtmp1=ind1;tmp2=ind2;tmp3=ind3;%确保ind1ind2ind3if(ind1ind2)&&(ind2ind3);elseif(ind1ind3)&&(ind3ind2)ind2=tmp3;ind3=tmp2;elseif(ind2ind1)&&(ind1ind3)ind1=tmp2;ind2=tmp1;elseif(ind2ind3)&&(ind3ind1)ind1=tmp2;ind2=tmp3;ind3=tmp1;elseif(ind3ind1)&&(ind1ind2)ind1=tmp3;ind2=tmp1;ind3=tmp2;elseif(ind3ind2)&&(ind2ind1)ind1=tmp3;ind2=tmp2;ind3=tmp1;endtmplist1=sol_new((ind1+1):(ind2-1));sol_new((ind1+1):(ind1+ind3-ind2+1))=...sol_new((ind2):(ind3));sol_new((ind1+ind3-ind2+2):ind3)=...tmplist1;end%检查是否满足约束%计算目标函数值(即内能)E_new=0;fori=1:(amount-1)E_new=E_new+...dist_matrix(sol_new(i),sol_new(i+1));end%再算上从最后一个城市到第一个城市的距离%E_new=E_new+...%dist_matrix(sol_new(amount),sol_new(1));ifE_newE_currentE_current=E_new;sol_current=sol_new;ifE_newE_best%把冷却过程中最好的解保存下来E_best=E_new;sol_best=sol_new;endelse%若新解的目标函数值小于当前解的,%则仅以一定概率接受新解ifrandexp(-(E_new-E_current)./t)E_current=E_new;sol_current=sol_new;elsesol_new=sol_current;endendendt=t.*a;%控制参数t(温度)减少为原来的a倍enddisp('最优解为:')disp(sol_best)disp('最短距离:')disp(E_best)figure()set(gcf,'Name','monituihuo-sol_best','Color','r')N=length(sol_best);scatter(coordinates(:,1),coordinates(:,2),50,'filled');holdonplot([coordinates(sol_best(1),1),coordinates(sol_best(N),1)],[coordinates(sol_best(1),2),coordinates(sol_best(N),2)])set(gca,'Color','g')holdonfori=2:Nplot([coordinates(sol_best(i-1),1),coordinates(sol_best(i),1)],[coordinates(sol_best(i-1),2),coordinates(sol_best(i),2)])holdonend