MATLAB 模擬退火SA
阿新 • • 發佈:2020-09-10
模擬退火首先從某個初始候選解開始,當溫度大於0時執行迴圈。
在迴圈中,通過隨機擾動產生一個新的解,然後求得新解和原解之間的能量差,如果差小於0,則採用新解作為當前解。
如果差大於0,則採用一個當前溫度與能量差成比例的概率來選擇是否接受新解。溫度越低,接受的概率越小,差值越大,同樣接受概率越小。
是否接受的概率用此公式計算:p=exp(-ΔE/T)。這裡ΔE為新解與原解的差,T為當前的溫度。
由於溫度隨迭代次數逐漸降低,因此獲得一個較差的解的概率較小。
典型的模擬退火演算法還使用了蒙特卡洛迴圈,在溫度降低之前,通過多次迭代來找到當前溫度下比較好的解。
這裡使用模擬退火解旅行商問題,因為這個問題本身是一個NP難問題,所以也就求不到最優解,不過應該可以求得一個比較好的解,然後再手工優化。
具體到程式中,這裡的隨機擾動就是隨機置換兩個城市的位置,能量就是旅行商路線的總長度。
演算法結果如下:
初始旅行商路線:
最終求得的旅行商路線:
每次迭代求得的旅行距離:
matlab程式碼如下:
main.m
1 clear all;close all;clc
2
3 n=20; %城市個數
4 temperature=100*n; %初始溫度
5 iter=100; %內部蒙特卡洛迴圈迭代次數
6
7 %隨機初始化城市座標
8 city=struct([]);
9 for i=1:n
10 city(i).x=floor(1 +100*rand());
11 city(i).y=floor(1+100*rand());
12 end
13
14 l=1; %統計迭代次數
15 len(l)=computer_tour(city,n); %每次迭代後的路線長度
16 netplot(city,n); %初始旅行路線
17
18 while temperature>0.001 %停止迭代溫度
19
20 for i=1:iter %多次迭代擾動,一種蒙特卡洛方法,溫度降低之前多次實驗
21 len1=computer_tour(city,n); %計算原路線總距離
22 tmp_city=perturb_tour(city,n); %產生隨機擾動
23 len2=computer_tour(tmp_city,n); %計算新路線總距離
24
25 delta_e=len2-len1; %新老距離的差值,相當於能量
26 if delta_e<0 %新路線好於舊路線,用新路線代替舊路線
27 city=tmp_city;
28 else %溫度越低,越不太可能接受新解;新老距離差值越大,越不太可能接受新解
29 if exp(-delta_e/temperature)>rand() %以概率選擇是否接受新解
30 city=tmp_city; %可能得到較差的解
31 end
32 end
33 end
34 l=l+1;
35 len(l)=computer_tour(city,n); %計算新路線距離
36 temperature=temperature*0.99; %溫度不斷下降
37
38 end
39 figure;
40 netplot(city,n); %最終旅行路線
41
42 figure;
43 plot(len)
computer_tour.m
1 function len=computer_tour(city,n) %計算路線總長度,每個城市只計算和下家城市之間的距離。
2 len=0;
3 for i=1:n-1
4 len=len+sqrt((city(i).x-city(i+1).x)^2+(city(i).y-city(i+1).y)^2);
5 end
6 len=len+sqrt((city(n).x-city(1).x)^2+(city(n).y-city(1).y)^2);
7 end
perturb_tour.m
1 function city=perturb_tour(city,n)
2
3 %隨機置換兩個不同的城市的座標
4 %產生隨機擾動
5 p1=floor(1+n*rand());
6 p2=floor(1+n*rand());
7
8 while p1==p2
9 p1=floor(1+n*rand());
10 p2=floor(1+n*rand());
11 end
12
13 tmp=city(p1);
14 city(p1)=city(p2);
15 city(p2)=tmp;
16
17 end
netplot.m
1 function netplot(city,n) %連線各城市,將路線畫出來
2 hold on;
3 for i=1:n-1
4 plot(city(i).x,city(i).y,'r*');
5 line([city(i).x city(i+1).x],[city(i).y city(i+1).y]); %只連線當前城市和下家城市
6 end
7
8 plot(city(n).x,city(n).y,'r*');
9 line([city(n).x city(1).x],[city(n).y city(1).y]); %最後一家城市連線第一家城市
10 hold off;
11 end