模拟退火算法应用——求飞机巡航最短路径

仅作自己学习使用

---

源码：

clc, clear, close all
%% 数据准备
% sj0=load('data12_1.txt');
% 数据是100个目标的经纬度，奇数列是经度，偶数列是维度
% 数据来源，《数学建模算法》第十二章
sj0 = [
   53.7121   15.3046   51.1758    0.0322   46.3253   28.2753   30.3313    6.9348
   56.5432   21.4188   10.8198   16.2529   22.7891   23.1045   10.1584   12.4819
   20.1050   15.4562    1.9451    0.2057   26.4951   22.1221   31.4847    8.9640
   26.2418   18.1760   44.0356   13.5401   28.9836   25.9879   38.4722   20.1731
   28.2694   29.0011   32.1910    5.8699   36.4863   29.7284    0.9718   28.1477
    8.9586   24.6635   16.5618   23.6143   10.5597   15.1178   50.2111   10.2944
    8.1519    9.5325   22.1075   18.5569    0.1215   18.8726   48.2077   16.8889
   31.9499   17.6309    0.7732    0.4656   47.4134   23.7783   41.8671    3.5667
   43.5474    3.9061   53.3524   26.7256   30.8165   13.4595   27.7133    5.0706
   23.9222    7.6306   51.9612   22.8511   12.7938   15.7307    4.9568    8.3669
   21.5051   24.0909   15.2548   27.2111    6.2070    5.1442   49.2430   16.7044
   17.1168   20.0354   34.1688   22.7571    9.4402    3.9200   11.5812   14.5677
   52.1181    0.4088    9.5559   11.4219   24.4509    6.5634   26.7213   28.5667
   37.5848   16.8474   35.6619    9.9333   24.4654    3.1644    0.7775    6.9576
   14.4703   13.6368   19.8660   15.1224    3.1616    4.2428   18.5245   14.3598
   58.6849   27.1485   39.5168   16.9371   56.5089   13.7090   52.5211   15.7957
   38.4300    8.4648   51.8181   23.0159    8.9983   23.6440   50.1156   23.7816
   13.7909    1.9510   34.0574   23.3960   23.0624    8.4319   19.9857    5.7902
   40.8801   14.2978   58.8289   14.5229   18.6635    6.7436   52.8423   27.2880
   39.9494   29.5114   47.5099   24.0664   10.1121   27.2662   28.7812   27.6659
    8.0831   27.6705    9.1556   14.1304   53.7989    0.2199   33.6490    0.3980
    1.3496   16.8359   49.9816    6.0828   19.3635   17.6622   36.9545   23.0265
   15.7320   19.5697   11.5118   17.3884   44.0398   16.2635   39.7139   28.4203
    6.9909   23.1804   38.3392   19.9950   24.6543   19.6057   36.9980   24.3992
    4.1591    3.1853   40.1400   20.3030   23.9876    9.4030   41.1084   27.7149
    ];

    
x=sj0(:,[1:2:8]); x=x(:);
y=sj0(:,[2:2:8]); y=y(:);
sj=[x y]; d1=[70,40]; 
xy=[d1;sj;d1]; sj=xy*pi/180; %角度化成弧度
d=zeros(102); %距离矩阵初始化
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';
path=[]; %巡航路径记录
long=inf; %巡航总长度记录

%% 求较好的初始解
generation = 1000; %迭代次数
for j=1:generation  
    path0=[1 1+randperm(100),102]; temp=0;
    for i=1:101
        temp=temp+d(path0(i),path0(i+1));
    end
    if temp<long
        path=path0; long=temp;
    end
end

%% 新解的产生（求最短路径）
e=0.1^30;L=20000;at=0.999;T=1;
for k=1:L  %退火过程
    c=2+floor(100*rand(1,2));  %产生新解
    c=sort(c); c1=c(1);c2=c(2);
    %计算代价函数值的增量
    df=d(path(c1-1),path(c2))+d(path(c1),path(c2+1))-...
        d(path(c1-1),path(c1))-d(path(c2),path(c2+1));
    if df<0 %接受准则
        path=[path(1:c1-1),path(c2:-1:c1),path(c2+1:102)]; long=long+df;
    elseif exp(-df/T)>=rand
        path=[path(1:c1-1),path(c2:-1:c1),path(c2+1:102)]; long=long+df;
    end
    T=T*at;
    if T<e
        break;
    end
end

%% 数据可视化
path, long % 输出巡航路径及路径长度
xx=xy(path,1); yy=xy(path,2);
plot(xx,yy,'-*') %画出巡航路径

效果图：