Matlab随笔之模拟退火算法

浏览数：45 / 时间：2015年06月08日

问题描述：

我方有一个基地，经度和纬度为（ 70,40）。假设我方飞机的速度为 1000 公里/小时。
我方派一架飞机从基地出发，侦察完敌方所有目标，再返回原来的基地。在敌方每一目
标点的侦察时间不计，求该架飞机所花费的时间（假设我方飞机巡航时间可以充分长）。
这是一个旅行商问题。我们依次给基地编号为 1，敌方目标依次编号为 2， 3，…，
101，最后我方基地再重复编号为 102（这样便于程序中计算）。距离矩阵 D = ( dij )102×102 ，
其中
dij 表示表示 i, j 两点的距离， i, j = 1,2,L ,102 ，这里 D 为实对称矩阵。则问题是
求一个从点 1 出发，走遍所有中间点，到达点 102 的一个最短路径。

算法如下：

%模拟退火算法之旅行商算法 
%ij.txt存放的是各点经纬度 
clc;clear 
load C:\Users\LittleLala\Desktop\ij.txt;%加载目标数据ij.txt存放25x8矩阵 
x=ij(:,1:2:8);%x为25x4矩阵 
x=x(:);%x化为10x1矩阵 
y=ij(:,2:2:8);%y为25x4矩阵 
y=y(:);%y化为100x1矩阵 
sj=[x y];%100x2矩阵 
d1=[70,40];%起始点 
sj=[d1;sj;d1];%首尾添加上出发点和终点 
sj=sj*pi/180;%角度化为弧度 
%计算距离 
d=zeros(102);%先建一个102x102的空的距离矩阵，再把对应距离填进去 
%A(R cos x1 cos y1 , R sin x1 cos y1 , R sin y1 ) 
%B(R cos x2 cos y2 , R sin x2 cos y2 , R sin y2 ) 
%其中 R = 6370 为地球半径。 
%A, B 两点的实际距离: 
%d = R arccos[cos(x1 − x2 ) cos y1 cos y2 + sin y1 sin y2 ] 
for i=1:101 
    for j=i+1:102 
        d(i,j)=6370*acos(cos(sj(i,1)-sj(j,1))*cos(sj(i,2))*cos(sj(j,2))... 
            +sin(sj(i,2))*sin(sj(j,2))); 
    end 
end 
%上述的d为上三角矩阵，根据对称性，化为对称矩阵,d矩阵包含各点间的距离信息 
d=d+d‘; 
S0=[]; 
Sum=inf; 
rand(‘state‘,sum(clock));%作用是定义一个随时间变化的初值 
%matlab里面的随机生成函数基本都是以rand为基函数通过函数关系式得到， 
%比如 normrnd，unidrnd等，你每次重启matlab后运行已编好的含随机数生成的函数你将得到相同的结果， 
%比如我的电脑上重启matlab运行unidrnd(100),每次的值都是82，这是因为rand函数的初值都一样， 
%所以为了避免上述问题经常在程序前运行或加命令rand(‘state‘,sum(clock))，这样重启matlab，运行随机数生成值就不同了。 
for i=1:1000%循环1000次，尽可能让所有排序情况都出现一遍 
    s=[1,1+randperm(100),102];%s为中间2到101的排列组合 
    temp=0; 
    for j=1:101 
        temp=temp+d(s(j),s(j+1)); 
    end 
    if temp<Sum 
       Sum=temp; 
       S0=s;%保存当前最小值的排序情况 
    end 
end 
e=0.1^30;%选定的终止温度 e = 10 −30 ，判断退火过程是否结束 
at=0.999;%选定的降温系数α进行降温即：T =αT 
L=20000; 
T=1; 
%退火过程 
for k=1:L 
c=2+floor(100*rand(1,2));%floor函数其功能是“向下取整”，或者说“向下舍入”，即取不大于x的最大整数 
%ceil函数向上取整 
c=sort(c);% [a,b]=sort(X)是按列从小到大排序，而 [a,b]=sort(X,2)是按行，b为排序情况 
c1=c(1);c2=c(2); 
%计算代价函数值 
df=d(S0(c1-1),S0(c2))+d(S0(c1),S0(c2+1))-d(S0(c1-1),S0(c1))-d(S0(c2),S0(c2+1)); 
if df<0%接受准则 
    S0=[S0(1:c1-1),S0(c2:-1:c1),S0(c2+1:102)]; 
    Sum=Sum+df; 
elseif exp(-df/T)>rand(1)%以概率exp(−df/T)接受新的路径%注意时elseif而不是else if 
    S0=[S0(1:c1-1),S0(c2:-1:c1),S0(c2+1:102)]; 
    Sum=Sum+df; 
end; 
T=T*at;%退火 
if T<e%达到终止温度 
    break; 
end; 
end; 
%输出路径和路径长度 
S0,Sum 
%画出其中一个巡航路径 
plot(sj(S0,1)/pi*180,sj(S0,2)/pi*180); 
hold on 
plot(sj(S0,1)/pi*180,sj(S0,2)/pi*180,‘rx‘); 
axis([-5,75,-5,45]);