x=sj(:,1:2:8);x=x(:);
y=sj(:,2:2:8);y=y(:);
sj=[x y];
d1=[70,40];
sj=[d1;sj;d1];
sj=sj*pi/180;
%距离矩阵d
d=zeros(102);
for i=1:101
for j=i+1:102
temp=cos(sj(i,1)-sj(j,1))*cos(sj(i,2))*cos(sj(j,2))+sin(sj(i,2))*sin(sj(j,2));
d(i,j)=6370*acos(temp);
end
end
d=d+d';
S0=[];Sum=inf;
rand('state',sum(clock));
for j=1:1000
S=[1 1+randperm(100),102];
temp=0;
for i=1:101
temp=temp+d(S(i),S(i+1));
end
if temp<Sum
S0=S;Sum=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(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)
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
% 表1 经度和纬度数据表
% 经度 纬度 经度 纬度 经度 纬度 经度 纬度
% 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
