[數學建模(二)模擬退火法與旅行商問題]

1.旅行商問題

        旅行商問題(Traveling Salesman Problem,TSP)，是由愛爾蘭數學家Sir William Rowan Hamilton和英國數學家Thomas Penyngton Kirkman在19世紀提出的數學問題。它的描述是這樣的：一名商人要到若干城市去推銷商品，已知城市個數和各城市間的路程（或旅費），要求找到一條從城市1出發，經過所有城市且每個城市只能訪問一次，最后回到城市1的路線，使總的路程（或旅費）最小。現已證明它屬於NP(Non-deterministic Polynomial---非確定多項式)難題。

       若設城市數目為n時，那么組合路徑數則為(n-1)!。當城市數量較小時，通過枚舉法就可以找出最短的路徑。隨着城市數目的不斷增大，組合路線數將呈指數級數規律急劇增長，以至達到無法計算的地步。此時，很難找到一個多項式時間復雜度的算法來求解該問題。

      TSP是一個典型的組合優化問題，是諸多領域內出現的多種復雜問題的集中概括和簡化形式，並且已成為各種啟發式的搜索、優化算法的間接比較標准。

 

 

圖1 當城市數量較少時，直接通過枚舉求解

2.模擬退火法

2.1算法

     模擬退火算法由KirkPatrick於1982提出，他將退火思想引入到組合優化領域,提出一種求解大規模組合優化問題的方法，對於NP-完全組合優化問題尤其有效。（化學上的退火方法講述比較復雜，有興趣的可以自行查找）

模擬退火算法可以分解為解空間、目標函數和初始解3部分。其基本思想是：

(1)初始化：初始溫度T(充分大)，初始解狀態s(是算法迭代的起點)，每個T值的迭代次數L；

(2)對k=1，……，L做第(3)至第6步；

(3)產生新解s′；

(4)計算增量cost=cost(s′)-cost(s)，其中cost(s)為評價函數；

(5)若t<0則接受s′作為新的當前解，否則以概率exp(-t′/T)接受s′作為新的當前解；

(6)如果滿足終止條件則輸出當前解作為最優解，結束程序。終止條件通常取為連續若干個新解都沒有被接受時終止算法；

(7)T逐漸減少，且T趨於0，然后轉第2步運算。

2.2 參數選擇

（1）溫度T的初始值設置。

    溫度T的初始值設置是影響模擬退火算法全局搜索性能的重要因素之一。初始溫度高，則搜索到全局最優解的可能性大，但因此要花費大量的計算時間；反之，則可節約計算時間，但全局搜索性能可能受到影響。實際應用過程中，初始溫度一般需要依據實驗結果進行若干次調整。

（2）溫度衰減函數的選取。

      衰減函數用於控制溫度的退火速度，一個常用的函數為：

               T(t+1)=aT(t) 

     式中a是一個非常接近於1的常數，t為降溫的次數。

（3）馬爾可夫鏈長度L的選取。

通常的原則是：在衰減參數T的衰減函數已選定的前提下，L的選取應遵循在控制參數的每一取值上都能恢復准平衡的原則。

3. TSP算法實現

3.1 TSP算法描述

   （1）TSP問題的解空間和初始解

    TSP的解空間S是遍訪每個城市恰好一次的所有回路，是所有城市排列的集合。TSP問題的解空間S可表示為{1,2,…,n}的所有排列的集合，即S = {(c₁,c₂,…,c_n) | ((c₁,c₂,…,c_n)為{1,2,…,n}的排列）}，其中每一個排列S_i表示遍訪n個城市的一個路徑，c_i= j表示在第i次訪問城市j。模擬退火算法的最優解與初始狀態無關，故初始解為隨機函數生成一個{1,2,…,n}的隨機排列作為S₀。

（2）目標函數

    TSP問題的目標函數即為訪問所有城市的路徑總長度，也可稱為代價函數：

                                              

    現在TSP問題的求解就是通過模擬退火算法求出目標函數C(c₁,c₂,…,c_n)的最小值，相應地，s^*= (c^*₁,c^*₂,…,c^*_n)即為TSP問題的最優解。

   （3）新解產生

新解的產生對問題的求解非常重要。新解可通過分別或者交替用以下2種方法產生：

①二變換法：任選序號u,v(設uvn)，交換u和v之間的訪問順序，若交換前的解為s_i= (c₁,c₂,…,c_u,…,c_v,…,c_n),交換后的路徑為新路徑，即：

②三變換法：任選序號u,v和ω(u≤vω)，將u和v之間的路徑插到ω之后訪問，若交換前的解為s_i= (c₁,c₂,…,c_u,…,c_v,…,c_ω,…,c_n)，交換后的路徑為的新路徑為：

（4）目標函數差

計算變換前的解和變換后目標函數的差值：

（5）Metropolis接受准則

根據目標函數的差值和概率exp(-ΔC′/T)接受s_i′作為新的當前解s_i，接受准則：

                            

3.2 TSP算法流程

    根據以上對TSP的算法描述，可以寫出用模擬退火算法解TSP問題的流程圖2 所示：

圖 2  TSP的模擬退火流程

4.MATLAB實現

clear
clc

coordinates = load('data.txt');  %加載數據

%設置參數
a = 0.99;   % 溫度衰減函數的參數
t0 = 97; tf = 3; t = t0;
Markov_length = 10000;  % Markov鏈長度
amount = size(coordinates,1);   % 城市的數目

% 通過向量化的方法計算距離矩陣
dist_matrix = zeros(amount, amount);   
coor_x_tmp1 = coordinates(:,1) * ones(1,amount);
coor_x_tmp2 = coor_x_tmp1';
coor_y_tmp1 = coordinates(:,2) * ones(1,amount);
coor_y_tmp2 = coor_y_tmp1';
dist_matrix = sqrt((coor_x_tmp1-coor_x_tmp2).^2 + (coor_y_tmp1-coor_y_tmp2).^2);  %存儲各個城市之間距離的矩陣

sol_new = 1:amount;         % 產生初始解
% sol_new是每次產生的新解；sol_current是當前解；sol_best是冷卻中的最好解；
E_current = inf;E_best = inf;       % E_current是當前解對應的回路距離；
% E_new是新解的回路距離；
% E_best是最優解的
sol_current = sol_new; sol_best = sol_new;
p = 1;

while t>=tf
    for r=1:Markov_length       % Markov鏈長度
        % 產生隨機擾動
        if (rand < 0.5) % 隨機決定是進行兩交換還是三交換
            % 兩交換
            ind1 = 0; ind2 = 0;
            while (ind1 == ind2)
                ind1 = ceil(rand.*amount);
                ind2 = ceil(rand.*amount);
            end
            tmp1 = sol_new(ind1);
            sol_new(ind1) = sol_new(ind2);
            sol_new(ind2) = tmp1;
        else
            % 三交換
            ind1 = 0; ind2 = 0; ind3 = 0;
            while (ind1 == ind2) || (ind1 == ind3)|| (ind2 == ind3) || (abs(ind1-ind2) == 1)
                ind1 = ceil(rand.*amount);
                ind2 = ceil(rand.*amount);
                ind3 = ceil(rand.*amount);
            end
            tmp1 = ind1;tmp2 = ind2;tmp3 = ind3;
            % 確保ind1 < ind2 < ind3
            if (ind1 < ind2) && (ind2 < ind3)
                ;
            elseif (ind1 < ind3) && (ind3 < ind2)
                ind2 = tmp3;ind3 = tmp2;
            elseif (ind2 < ind1) && (ind1 < ind3)
                ind1 = tmp2;ind2 = tmp1;
            elseif (ind2 < ind3) && (ind3 < ind1)
                ind1 = tmp2;ind2 = tmp3; ind3 = tmp1;
            elseif (ind3 < ind1) && (ind1 < ind2)
                ind1 = tmp3;ind2 = tmp1; ind3 = tmp2;
            elseif (ind3 < ind2) && (ind2 < ind1)
                ind1 = tmp3;ind2 = tmp2; ind3 = tmp1;
            end

            tmplist1 = sol_new((ind1+1):(ind2-1));
            sol_new((ind1+1):(ind1+ind3-ind2+1)) = sol_new((ind2):(ind3));
            sol_new((ind1+ind3-ind2+2):ind3) = tmplist1;
        end

        %檢查是否滿足約束

        % 計算目標函數值（即內能）
        E_new = 0;
        for i = 1 : (amount-1)
            E_new = E_new +  dist_matrix(sol_new(i),sol_new(i+1));
        end
        % 再算上從最后一個城市到第一個城市的距離
        E_new = E_new + dist_matrix(sol_new(amount),sol_new(1));

        if E_new < E_current
            E_current = E_new;
            sol_current = sol_new;
            if E_new < E_best
                % 把冷卻過程中最好的解保存下來
                E_best = E_new;
                sol_best = sol_new;
            end
        else
            % 若新解的目標函數值小於當前解的，
            % 則僅以一定概率接受新解
            if rand < exp(-(E_new-E_current)./t)
                E_current = E_new;
                sol_current = sol_new;
            else
                sol_new = sol_current;
            end
        end
    end
    t=t.*a;     % 控制參數t（溫度）減少為原來的a倍
end

for i=1:length(coordinates)
    plot(coordinates(i,1),coordinates(i,2),'r*');
    hold on;
end;

x=coordinates([sol_best sol_best(1)],1);
y=coordinates([sol_best sol_best(1)],2);
plot(x,y);
disp('最優解為：')
disp(sol_best)
disp('最短距離：')
disp(E_best)

 以上代碼來自：http://blog.csdn.net/zhangzhengyi03539/article/details/46673545

txt數據：

565 575
25 185
345 750
945 685
845 655
880 660
25 230
525 1000
580 1175
650 1130
1605 620
1220 580
1465 200
1530 5
845 680
725 370
145 665
415 635
510 875
560 365
300 465
520 585
480 415
835 625
975 580
1215 245
1320 315
1250 400
660 180
410 250
420 555
575 665
1150 1160
700 580
685 595
685 610
770 610
795 645
720 635
760 650
475 960
95 260
875 920
700 500
555 815
830 485
1170 65
830 610
605 625
595 360
1340 725
1740 245