运筹优化 | 模拟退火算法解决旅行商问题详解_我方有一个基地,经度纬度坐标为(70,40)

一、题目

        我方有一个基地,经度和纬度为(70,40)。假设我方飞机的速度为 400km/h。我方派一架飞机从基地出发,侦察完所有目标,再返回原来的基地。在每一目标点的侦察时间不计.求该架飞机所花费的时间(假设我方飞机巡航时间可以充分长)。经纬度数据如下图，可复制数据见文末。

二、题目解析

         对于实际距离、交换法、路径差的理解如下： 

代码如下：

clear;clc
sj=xlsread("sj.xlsx");% 加载excel数据 
 
x=sj(:,1);% 第一列为经度
y=sj(:,2);% 第二列为纬度
dl=[70,40];% 起点坐标
sj=[dl;sj;dl];
% 我们设置第一行为起点，最后一行还是起点，从起点出发遍历100个坐标，又回到起点，共102个坐标
 
sj=sj*pi/180;% 角度制转化为弧度制
d=zeros(102,102);% 初始化距离矩阵
% a向量和b向量夹角的余弦值等于 a向量点乘b向量 除以 a向量的模乘上b向量的模
% 那么a向量与b向量的夹角等于arccos（a向量点乘b向量 除以 a向量的模乘上b向量的模）
% 由于扇形弧长等于夹角乘半径
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=[];% 初始化路径矩阵，这个不能用zeros去初始化
long=inf;% 长度初始化
 
rand('state',sum(clock));% 用年月日时分秒的和来设置随机流
% 产生较好的初始解：path数组和long
% 其判断的根据在于long尽可能接近最短巡航长度，从而在模拟退火过程中更容易找到全局最优解
for j=1:10000% 模拟1万次 
    path0=randperm(102);% 把1到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
 
T=1;% 初始温度
e=0.1^30;% 终止温度
L=20000;% 迭代次数
at=0.999;% 温度衰减系数
 
for k=1:L% 退火过程
    c=2+floor(100*rand(1,2));
    % rand生成0-1之间的随机数，乘上100并且向下取整是避免出现小数以及大于100的情况
    % 由于MATLAB下标从1开始，并且下面的路径差里面的下标还有减1，所以加2是为了避免出现小于等于0的情况
    c=sort(c);
    % 生成两个下标，让它们交换
    u=c(1);
    v=c(2);
    % 路径差
    df=d(path(u -1),path(v)) +d(path(u),path(v +1)) - d(path(u -1),path(u))-d(path(v),path(v+1));
    if df<0% 路径差小于0代表新路径长度小于原路径长度，那么接受
        path=[path(1:u -1),path(v:-1:u),path(v +1:102)];
        long =long +df;
    elseif rand<=exp(-df/T)% 如果不满足if的条件，那么执行elseif：生成的随机数小于概率
        path=[path(1:u -1),path(v:-1:u),path(v +1:102)];% 那就接受
    end
    T=T*at;% 此时更新温度（退火）
    if T<e% 温度到达终止温度，结束退火过程
        break;
    end
end
 
path,long % 输出路径数组和最短巡航长度
time=long/1000% 输出最短巡航时间
xx=sj(path,1);
yy=sj(path,2);
plot(xx,yy,'-*')

结果如下：

此次模拟的最短路径长度为： 4.5687e+04 km 

最短时间为：114.218 h

经纬度数据如下：