模拟退火算法举例及其matlab实现

已知敌方100个目标的经度、纬度如表1所示。
表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 
-274-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 
我方有一个基地，经度和纬度为（70,40）。假设我方飞机的速度为1000公里/小时。
我方派一架飞机从基地出发，侦察完敌方所有目标，再返回原来的基地。在敌方每一目
标点的侦察时间不计，求该架飞机所花费的时间（假设我方飞机巡航时间可以充分长）。
这是一个旅行商问题。我们依次给基地编号为1，敌方目标依次编号为2，3，…，
101，最后我方基地再重复编号为102（这样便于程序中计算）。距离矩阵 102 102
) (
× = ij d D ，
其中ij d 表示表示 j i, 两点的距离， 102 , , 2 , 1 , L = j i ，这里D为实对称矩阵。则问题是
求一个从点1出发，走遍所有中间点，到达点102的一个最短路径。
上面问题中给定的是地理坐标（经度和纬度），我们必须求两点间的实际距离。设
B A, 两点的地理坐标分别为 ) , (
1 1
y x ， ) , (
2 2
y x ，过 B A, 两点的大圆的劣弧长即为两点
的实际距离。以地心为坐标原点O，以赤道平面为XOY平面，以0度经线圈所在的平
面为XOZ平面建立三维直角坐标系。则 B A, 两点的直角坐标分别为：
) sin , cos sin , cos cos (
1 1 1 1 1
y R y x R y x R A
) sin , cos sin , cos cos (
2 2 2 2 2
y R y x R y x R B
其中 6370 = R 为地球半径。
B A, 两点的实际距离
⎟
⎟
⎠
⎞
⎜
⎜
⎝
⎛
⋅
⋅
=
OB
OB
R d
OA
OA
arccos ，
化简得
] sin sin cos cos ) ( arccos[cos
2 1 2 1 2 1
y y y y x x R d + − = 。
求解的模拟退火算法描述如下：
（1）解空间
解空间S可表为{ 102 , 101 , , 2 , 1 L }的所有固定起点和终点的循环排列集合，即
} 102 , } 101 , , 3 , 2 { ) , , ( , 1 | ) , , {(
102 101 2 1 102 1 = = = π π π π π π 的循环排列 为 L L L S
其中每一个循环排列表示侦察100个目标的一个回路， j
i = π 表示在第i次侦察j点，
初始解可选为 ) 102 , , 2 , 1 ( L ，本文中我们使用Monte Carlo方法求得一个较好的初始解。
（2）目标函数
此时的目标函数为侦察所有目标的路径长度或称代价函数。我们要求
∑
=
+
=
101
1
2 1 1 102
) , , , ( min
i
i i
d f
π π π π π L
而一次迭代由下列三步构成：
（3）新解的产生
①2变换法 
-275-任选序号  v u, （ v u< ）交换u与v之间的顺序，此时的新路径为：
102 1 1 1 1 1 π π π π π π π π L L L + + + − v u u v v u
②3变换法
任选序号 v u, 和w，将u和v之间的路径插到w之后，对应的新路径为（设
w v u < < ）
102 1 1 1 1 π π π π π π π π L L L L + + − w v u w v u
（4）代价函数差
对于2变换法，路径差可表示为
) ( ) (
1 1 1 1 + − + −
+ − + = Δ
v v u u v u v u
d d d d f
π π π π π π π π
（5）接受准则
⎩
⎨
⎧
≥ Δ Δ −
< Δ
=
0 ) / exp(
0   1
f T f
f
P
如果 0 < Δf ，则接受新的路径。否则，以概率 ) / exp( T f Δ −   接受新的路径，即若
) / exp( T f Δ − 大于0到1之间的随机数则接受。
（6）降温
利用选定的降温系数α进行降温即： T T α ← ，得到新的温度，这里我们取
999 . 0 = α 。
（7）结束条件
用选定的终止温度
30
10
−
= e ，判断退火过程是否结束。若 e T< ，算法结束，输出
当前状态。
我们编写如下的matlab程序如下：
clc,clear 
load sj.txt %加载敌方100个目标的数据，数据按照表格中的位置保存在纯文本
文件sj.txt中
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; 
-276-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