模擬退火演算法舉例及其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