1. 程式人生 > 實用技巧 >MATLAB 模擬退火SA

MATLAB 模擬退火SA

模擬退火首先從某個初始候選解開始,當溫度大於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