模拟退火算法TSP问题matlab

TSP问题数据，52个城市的坐标：城市编号X坐标Y坐标城市编号X坐标Y坐标城市编号X坐标Y坐标156557519510875377706102251852056036538795645334575021300465397206354945685225205854076065058456552348041541475960688066024835625429526072523025975580438759208525100026121524544700500958011752713203154555581510650113028125040046830485111605620296601804711706512122058030410250488306101314652003142055549605625141530532575665505953601584568033115011605113407251672537034700580521740245171456653568559553184156353668561052clearclca=0.99;%温度衰减函数的参数t0=97;tf=3;t=t0;Markov_length=10000;%Markov链长度coordinates=[1565.0575.0;225.0185.0;3345.0750.0;4945.0685.0;5845.0655.0;6880.0660.0;725.0230.0;8525.01000.0;9580.01175.0;10650.01130.0;111605.0620.0;121220.0580.0;131465.0200.0;141530.05.0;15845.0680.0;16725.0370.0;17145.0665.0;18415.0635.0;19510.0875.0;20560.0365.0;21300.0465.0;22520.0585.0;23480.0415.0;24835.0625.0;25975.0580.0;261215.0245.0;271320.0315.0;281250.0400.0;29660.0180.0;30410.0250.0;31420.0555.0;32575.0665.0;331150.01160.0;34700.0580.0;35685.0595.0;36685.0610.0;37770.0610.0;38795.0645.0;39720.0635.0;40760.0650.0;41475.0960.0;4295.0260.0;43875.0920.0;44700.0500.0;45555.0815.0;46830.0485.0;471170.065.0;48830.0610.0;49605.0625.0;50595.0360.0;511340.0725.0;521740.0245.0;];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)