matlab練習程序（模擬退火SA）

模擬退火首先從某個初始候選解開始，當溫度大於0時執行循環。

在循環中，通過隨機擾動產生一個新的解，然后求得新解和原解之間的能量差，如果差小於0，則采用新解作為當前解。

如果差大於0，則采用一個當前溫度與能量差成比例的概率來選擇是否接受新解。溫度越低，接受的概率越小，差值越大，同樣接受概率越小。

是否接受的概率用此公式計算：p=exp(-ΔE/T)。這里ΔE為新解與原解的差，T為當前的溫度。

由於溫度隨迭代次數逐漸降低，因此獲得一個較差的解的概率較小。

典型的模擬退火算法還使用了蒙特卡洛循環，在溫度降低之前，通過多次迭代來找到當前溫度下比較好的解。

這里使用模擬退火解旅行商問題，因為這個問題本身是一個NP難問題，所以也就求不到最優解，不過應該可以求得一個比較好的解，然后再手工優化。

具體到程序中，這里的隨機擾動就是隨機置換兩個城市的位置，能量就是旅行商路線的總長度。

算法結果如下：

初始旅行商路線：

最終求得的旅行商路線：

每次迭代求得的旅行距離：

matlab代碼如下：

main.m

clear all;close all;clc

n=20;                   %城市個數
temperature=100*n;      %初始溫度
iter=100;               %內部蒙特卡洛循環迭代次數

%隨機初始化城市坐標
city=struct([]);
for i=1:n
    city(i).x=floor(1+100*rand()); 
    city(i).y=floor(1+100*rand());
end

l=1;                            %統計迭代次數
len(l)=computer_tour(city,n);   %每次迭代后的路線長度  
netplot(city,n);                %初始旅行路線

while temperature>0.001     %停止迭代溫度
    
    for i=1:iter     %多次迭代擾動，一種蒙特卡洛方法，溫度降低之前多次實驗
        len1=computer_tour(city,n);         %計算原路線總距離
        tmp_city=perturb_tour(city,n);      %產生隨機擾動
        len2=computer_tour(tmp_city,n);     %計算新路線總距離
        
        delta_e=len2-len1;  %新老距離的差值，相當於能量
        if delta_e<0        %新路線好於舊路線，用新路線代替舊路線
            city=tmp_city;
        else                        %溫度越低，越不太可能接受新解；新老距離差值越大，越不太可能接受新解
            if exp(-delta_e/temperature)>rand() %以概率選擇是否接受新解
                city=tmp_city;      %可能得到較差的解
            end
        end        
    end
    l=l+1;
    len(l)=computer_tour(city,n);   %計算新路線距離
    temperature=temperature*0.99;   %溫度不斷下降
  
end  
figure;
netplot(city,n);    %最終旅行路線

figure;
plot(len)  

computer_tour.m

function len=computer_tour(city,n)   %計算路線總長度，每個城市只計算和下家城市之間的距離。
    len=0;
    for i=1:n-1
        len=len+sqrt((city(i).x-city(i+1).x)^2+(city(i).y-city(i+1).y)^2);        
    end
    len=len+sqrt((city(n).x-city(1).x)^2+(city(n).y-city(1).y)^2);
end

perturb_tour.m

function city=perturb_tour(city,n)  
    
    %隨機置換兩個不同的城市的坐標
    %產生隨機擾動
    p1=floor(1+n*rand());
    p2=floor(1+n*rand());

    while p1==p2
        p1=floor(1+n*rand());
        p2=floor(1+n*rand());    
    end
    
    tmp=city(p1);
    city(p1)=city(p2);
    city(p2)=tmp;

end

netplot.m

function netplot(city,n)        %連線各城市，將路線畫出來
    hold on;
    for i=1:n-1
        plot(city(i).x,city(i).y,'r*');  
        line([city(i).x city(i+1).x],[city(i).y city(i+1).y]);  %只連線當前城市和下家城市     
    end

    plot(city(n).x,city(n).y,'r*');  
    line([city(n).x city(1).x],[city(n).y city(1).y]);     %最后一家城市連線第一家城市
    hold off;
end