基于遺傳算法的TSP算法

文章目錄

• 一、理論基礎
• 二、案例背景
• • 1，問題描述
  • 2，解決思路和步驟
  • • (1).算法流程
    • (2).遺傳算法實作
• 三、MATLAB程式實作
• • (1).種群初始化
  • (2).适應度函數
  • (3).選擇操作
  • (4).交叉操作
  • (5).變異操作
  • (6).進化逆轉操作
  • (7).畫路線軌迹圖
  • (8).遺傳算法主函數
  • (9).結果分析
• 四、遺傳算法的改進
• • 1. 使用精英政策
  • 2. 使用進化逆轉操作
• 五、算法的局限性
• 六、參考文獻

一、理論基礎

TSP(traveling salesman problem,旅行商問題)是典型的NP完全問題，即其最壞情況下的時間複雜度随着問題規模的增大按指數方式增長，到目前為止還未找到一個多項式時間的有效算法。

TSP問題可描述為：已知 n n n個城市互相之間的距離，某一旅行商從某個城市出發通路每個城市有且僅有一次，最後回到出發城市，如何安排才使其所走路線距離最短。簡言之，就是尋找一條最短的周遊 n n n個城市的路徑，或者說搜尋自然子集 X = { 1 , 2 , ⋯   , n } X=\{1,2,\dotsm,n\} X={1,2,⋯,n}( X X X的元素表示對 n n n個城市的編号)的一個排列 π ( X ) = { V 1 , V 2 , ⋯   , V n } π(X)=\{V_1,V_2,\dotsm,V_n\} π(X)={V1​,V2​,⋯,Vn​}，使得 T d = ∑ i = 1 n + 1 d ( V i , V i + 1 ) + d ( V n , V 1 ) T_d=\sum_{i=1}^{n+1} {d(V_i,V_{i+1})}+d(V_n,V_1) Td​=i=1∑n+1​d(Vi​,Vi+1​)+d(Vn​,V1​)取得最小值，其中 d ( V i , V i + 1 ) d(V_i,V_{i+1}) d(Vi​,Vi+1​)表示城市 V i V_i Vi​到城市 V i + 1 V_{i+1} Vi+1​的距離。

二、案例背景

1，問題描述

本案例以14個城市為例，假定14個城市的位置坐标如表1所列。尋找出一條最短的周遊14個城市的路徑。

表1 14個城市的位置坐标

2，解決思路和步驟

(1).算法流程

遺傳算法TSP問題的流程圖如圖1所示。

圖1 遺傳算法TSP問題求解的流程圖

(2).遺傳算法實作

<1>編碼

采用整數排列編碼方法。對于 n n n個城市的TSP問題，染色體分為 n n n段，其中每一段為對應城市的編号，對10個城市的TSP問題 1 , 2 , 3 , 4 , 5 , 6 , 7 , 8 , 9 , 10 {1,2,3,4,5,6,7,8,9,10} 1,2,3,4,5,6,7,8,9,10，則 ∣ 1 ∣ 10 ∣ 2 ∣ 4 ∣ 5 ∣ 6 ∣ 8 ∣ 7 ∣ 9 ∣ 3 |1|10|2|4|5|6|8|7|9|3 ∣1∣10∣2∣4∣5∣6∣8∣7∣9∣3就是一個合法的染色體。

<2>種群初始化

在完成染色體編碼以後，必須産生一個初始種群作為起始解，是以首先需要決定初始化種群的數目。初始化種群的數目一般根據經驗得到，一般情況下種群的數量視城市規模的大小而定，其取值在50~200之間浮動。

<3>适應度函數

設 k 1 ∣ k 2 ∣ ⋯ ∣ k i ∣ ⋯ ∣ k n ∣ k_1|k_2|\dotsm|k_i|\dotsm|k_n| k1​∣k2​∣⋯∣ki​∣⋯∣kn​∣為一個采用整數編碼的染色體， D k i k j D_{k_i k_j} Dki​kj​​為城市 k i k_i ki​到城市 k j k_j kj​的距離，則該個體的适應度為 f i t n e s s = 1 ∑ i = 1 n − 1 D k i k j + D k n k 1 fitness=\frac {1}{ \displaystyle\sum_{i=1}^{n-1} {D_{k_i k_j}}+D_{k_nk_1}} fitness=i=1∑n−1​Dki​kj​​+Dkn​k1​​1​即适應度函數為恰好走遍 n n n個城市再回到出發城市的距離的倒數。優化的目标就是選擇适應度函數值盡可能大的染色體，适應度函數值越大的染色體越優質，反之越劣質。

<4>選擇操作

選擇操作即從舊群體中以一定機率選擇個體到新群體中，個體被選中的機率跟适應度值有關，個體适應度值越大，被選中的機率越大。

<5>交叉操作

采用部分映射雜交，确定交叉操作的父代，将父代樣本兩兩分組，每組重複以下過程(假定城市數為10)：

①産生兩個 [ 1 , 10 ] [1,10] [1,10]區間内的随機整數 r 1 r_1 r1​和 r 2 r_2 r2​，确定兩個位置，對兩位置的中間資料進行交叉，如 r 1 = 4 , r 2 = 7 r_1=4,r_2=7 r1​=4,r2​=7

9 5 1 ∣ 3 7 4 2 ∣ 10 8 6 9\quad5\quad1|3\quad7\quad4\quad2|10\quad8\quad6\quad 951∣3742∣1086 10 5 4 ∣ 6 3 8 7 ∣ 2 1 9 10\quad5\quad4|6\quad3\quad8\quad7|2\quad1\quad9\quad 1054∣6387∣219交叉為 9 5 1 ∣ 6 3 8 7 ∣ 10 ∗ ∗ 9\quad5\quad1|6\quad3\quad8\quad7|10\quad*\quad*\quad 951∣6387∣10∗∗ 10 5 ∗ ∣ 3 7 4 2 ∣ ∗ 1 9 10\quad5\quad*|3\quad7\quad4\quad2|*\quad1\quad9\quad 105∗∣3742∣∗19②交叉後，同一個個體中有重複的城市編号，不重複的數字保留，有沖突的數字(帶 ∗ * ∗位置)采用部分映射的方法消除沖突，即利用中間段的對應關系進行映射。結果為 9 5 1 ∣ 6 3 8 7 ∣ 10 4 2 9\quad5\quad1|6\quad3\quad8\quad7|10\quad4\quad2\quad 951∣6387∣1042 10 5 8 ∣ 3 7 4 2 ∣ 6 1 9 10\quad5\quad8|3\quad7\quad4\quad2|6\quad1\quad9\quad 1058∣3742∣619<6>變異操作

變異政策采取随機選取兩個點，将其對換位置。産生兩個 [ 1 , 10 ] [1,10] [1,10]範圍内的随機整數 r 1 r_1 r1​和 r 2 r_2 r2​，确定兩個位置，将其對換位置，如 r 1 = 4 , r 2 = 7 r_1=4,r_2=7 r1​=4,r2​=7 9 5 1 ∣ 6 ∣ 3 8 ∣ 7 ∣ 10 4 2 9\quad5\quad1|6|\quad3\quad8\quad|7|10\quad4\quad2\quad 951∣6∣38∣7∣1042變異後為 9 5 1 ∣ 7 ∣ 3 8 ∣ 6 ∣ 10 4 2 9\quad5\quad1|7|\quad3\quad8\quad|6|10\quad4\quad2\quad 951∣7∣38∣6∣1042<7>進化逆轉操作

為改善遺傳算法的局部搜尋能力，在選擇、交叉、變異之後引進連續多次的進化逆轉操作。這裡的“進化”是指逆轉算子的單方向性，即隻有經過逆轉後，适應度值有提高的才接受下來，否則逆轉無效。

産生兩個 [ 1 , 10 ] [1,10] [1,10]區間内的随機整數 r 1 r_1 r1​和 r 2 r_2 r2​，确定兩個位置，将其對換位置，如 r 1 = 4 , r 2 = 7 r_1=4,r_2=7 r1​=4,r2​=7 9 5 1 ∣ 7 3 8 ∣ 6 10 4 2 9\quad5\quad1|7\quad3\quad8|6\quad10\quad4\quad2\quad 951∣738∣61042進化逆轉後為 9 5 1 ∣ 8 3 7 ∣ 6 10 4 2 9\quad5\quad1|8\quad3\quad7|6\quad10\quad4\quad2\quad 951∣837∣61042對每個個體進行交叉變異，然後代入适應度函數進行評估， x x x選擇出适應度值大的個體進行下一代的交叉和變異以及進化逆轉操作。循環操作：判斷是否滿足設定的最大遺傳代數MAXGEN ，不滿足則跳入适應度值的計算；否則，結束遺傳操作。

三、MATLAB程式實作

(1).種群初始化

種群初始化函數InitPop的代碼：

function Chrom = InitPop(NIND,N)
%% 初始化種群
% 輸入：
% NIND：種群大小
% N：   個體染色體長度（這裡為城市的個數）  
% 輸出：
% 初始種群
Chrom = zeros(NIND, N);         % 用于存儲種群
for i = 1:NIND
    Chrom(i, :) = randperm(N);    % 随機生成初始種群
end
           

(2).适應度函數

求種群個體的适應度函數Fitness的代碼：

function FitnV = Fitness(len)
%% 适應度函數     
% 輸入：
% 個體的長度（TSP的距離）
% 輸出：
% 個體的适應度值
FitnV = 1./len;
           

(3).選擇操作

選擇操作函數Select的代碼：

function SelCh = Select(Chrom,FitnV,GGAP)
%% 選擇操作
% 輸入
% Chrom 種群
% FitnV 适應度值
% GGAP：選擇機率
% 輸出
% SelCh  被選擇的個體
NIND = size(Chrom,1);                         % 種群個數
NSel = max(floor(NIND*GGAP+.5), 2);   % 被選中個體的數目
ChrIx = Sus(FitnV,NSel);                 % 被選中的個體索引号
SelCh = Chrom(ChrIx, :);                  % 被選中的個體
           

其中，函數Sus的代碼為：

function NewChrIx = Sus(FitnV,Nsel)
% 輸入:
% FitnV  個體的适應度值
% Nsel   被選擇個體的數目
% 輸出:
% NewChrIx  被選擇個體的索引号
% Identify the population size (Nind)
   [Nind,ans] = size(FitnV);
% Perform stochastic universal sampling
   cumfit = cumsum(FitnV);
   trials = cumfit(Nind) / Nsel * (rand + (0:Nsel-1)');
   Mf = cumfit(:, ones(1, Nsel));
   Mt = trials(:, ones(1, Nind))';
   [NewChrIx, ans] = find(Mt < Mf & [ zeros(1, Nsel); Mf(1:Nind-1, :) ] <= Mt);
% Shuffle new population
   [ans, shuf] = sort(rand(Nsel, 1));
   NewChrIx = NewChrIx(shuf);   
% End of function
           

(4).交叉操作

交叉操作函數Recombin的代碼：

function SelCh = Recombin(SelCh,Pc)
%% 交叉操作
% 輸入
% SelCh  被選擇的個體
% Pc     交叉機率
% 輸出：
% SelCh 交叉後的個體
NSel = size(SelCh, 1);
for i = 1:2:NSel-mod(NSel, 2)
    if Pc >= rand         %  交叉機率Pc
        [SelCh(i, :), SelCh(i+1, :)] = intercross(SelCh(i, :), SelCh(i+1, :));
    end
end
function [a, b] = intercross(a, b)
%% 交叉
% 輸入：
% a和b為兩個待交叉的個體
% 輸出：
% a和b為交叉後得到的兩個個體
L = length(a);
r1 = randsrc(1, 1, [1:L]);
r2 = randsrc(1, 1, [1:L]);
if r1 ~= r2
    a0 = a; b0 = b;
    s = min([r1, r2]);
    e = max([r1, r2]);
    for i = s:e
        a1 = a; b1 = b;
        a(i) = b0(i);
        b(i) = a0(i);
        x = find(a == a(i));
        y = find(b == b(i));
        i1 = x(x ~= i);
        i2 = y(y ~= i);
        if ~isempty(i1)
            a(i1) = a1(i);
        end
        if ~isempty(i2)
            b(i2) = b1(i);
        end
    end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %% 交叉算法采用部分比對交叉
% 這部分代碼有問題
% function [a, b] = intercross(a, b)
% L = length(a);
% r1 = ceil(rand * L);
% r2 = ceil(rand * L);
% if r1 ~= r2
%     s = min([r1,r2]);
%     e = max([r1,r2]);
%     a1 = a; b1 = b;
%     a(s:e) = b1(s:e);
%     b(s:e) = a1(s:e);
%     for i = [setdiff(1:L, s:e)]
%         [tf, loc] = ismember(a(i), a(s:e));
%         if tf
%             a(i) = a1(loc+s-1);
%         end
%         [tf, loc] = ismember(b(i), b(s:e));
%         if tf
%             b(i) = b1(loc+s-1);
%         end
%     end
% end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% function [a, b] = intercross(a, b)
% L = length(a);
% r1 = ceil(rand * L);
% r2 = ceil(rand * L);
% if r1 ~= r2
%     s = min([r1, r2]);
%     e = max([r1, r2]);
%     a1 = a; b1 = b;
%     a(s:e) = b1(s:e);
%     b(s:e) = a1(s:e);
%     for i = [setdiff(1:L, s:e)]
%         a2 = a1(~ismember(a1, a));
%         b2 = b1(~ismember(b1, b));
%         if ~isempty(a2) 
%             if ~isempty(find(a(s:e) == a(i), 1))
%                a(i) = a2(1);
%             end
%         end
%         if ~isempty(b2)
%             if ~isempty( find(b(s:e) == b(i), 1))
%                 b(i) = b2(1);
%             end
%         end
%     end
% end
           

(5).變異操作

變異操作函數Mutate的代碼：

function SelCh = Mutate(SelCh, Pm)
%% 變異操作
% 輸入：
% SelCh  被選擇的個體
% Pm     變異機率
% 輸出：
% SelCh 變異後的個體
[NSel, L] = size(SelCh);
for i = 1:NSel
    if Pm >= rand
        R = randperm(L);
        SelCh(i, R(1:2)) = SelCh(i, R(2:-1:1));
    end
end

           

(6).進化逆轉操作

進化逆轉操作Reverse的代碼：

function SelCh=Reverse(SelCh,D)
%% 進化逆轉函數
% 輸入
% SelCh 被選擇的個體
% D     個城市的距離矩陣
% 輸出
% SelCh  進化逆轉後的個體
[row, col] = size(SelCh);
ObjV = PathLength(D, SelCh);        % 計算路徑長度
SelCh1 = SelCh;
for i = 1:row
    r1 = randsrc(1, 1, [1:col]);
    r2 = randsrc(1, 1, [1:col]);
    mininverse = min([r1 r2]);
    maxinverse = max([r1 r2]);
    SelCh1(i, mininverse:maxinverse) = SelCh1(i, maxinverse:-1:mininverse);
end
ObjV1 = PathLength(D, SelCh1);      % 計算路徑長度
index = ObjV1 < ObjV;
SelCh(index, :)= SelCh1(index, :);      % 取距離較短者
           

(7).畫路線軌迹圖

畫出所給路線的軌迹圖函數DrawPath的代碼：

function DrawPath(Chrom,X)
%% 畫路徑函數
%輸入
% Chrom  待畫路徑   
% X          各城市坐标位置
R = Chrom(1, :);       %一個随機解(個體)
figure;
hold on
plot([X(R, 1); X(R(1), 1)], [X(R, 2); X(R(1), 2)], 'o', 'color', [0.5,0.5,0.5])
plot(X(R(1), 1), X(R(1), 2), 'rv', 'MarkerSize', 20)
plot(X(R(end), 1), X(R(end), 2), 'rs', 'MarkerSize', 20)
text(X(R(1),1)+0.05, X(R(1),2)+0.05, ['   起點' num2str(R(1))], 'color', [1,0,0]);
text(X(R(end),1)+0.05, X(R(end),2)+0.05, ['   終點' num2str(R(end))], 'color', [1,0,0]);
for i = 1:size(X,1)
    if R(1) ~= i && R(end) ~= i
        text(X(i,1)+0.05, X(i,2)+0.05, num2str(i), 'color', [1,0,0]);
    end
end
A = [X(R, :); X(R(1), :)];
row=size(A,1);
for i=2:row
    [arrowx, arrowy] = dsxy2figxy(gca, A(i-1:i,1), A(i-1:i,2));     % 坐标轉換
    annotation('textarrow', arrowx, arrowy, 'HeadWidth', 8, 'color', [0,0,1]);
end
hold off
xlabel('橫坐标')
ylabel('縱坐标')
title('軌迹圖')
box on
           

(8).遺傳算法主函數

主函數名為GA_TSP，代碼如下：

clear
clc
close all
load CityPosition1.mat
%% 初始化參數
NIND = 100;                  % 種群大小
MAXGEN = 200;            % 最大遺傳代數
Pc = 0.9;                          % 交叉機率
Pm = 0.05;                       % 變異機率
GGAP = 0.9;                     % 代溝(Generation gap)
D = Distance(X);                 % 生成距離矩陣
N = size(D, 1);                     % (14*14)
%% 初始化種群
Chrom = InitPop(NIND, N);
%% 在二維圖上畫出所有坐标點
% figure
% plot(X(:,1),X(:,2),'o');
%% 畫出随機解的路線圖
DrawPath(Chrom(1,:), X)
%% 輸出随機解的路線和總距離
disp('初始種群中的一個随機值:')
OutputPath(Chrom(1,:));
Rlength = PathLength(D,Chrom(1,:));
disp(['總距離：', num2str(Rlength)]);
disp('~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~')
%% 優化
gen = 0;
figure;
hold on; box on
xlim([0,MAXGEN])
title('優化過程')
xlabel('代數')
ylabel('最優值')
ObjV = PathLength(D,Chrom);  %計算路線長度
preObjV = min(ObjV);
while gen < MAXGEN
    %% 計算适應度
    ObjV = PathLength(D,Chrom);  % 計算路線長度
    % fprintf('%d   %1.10f\n',gen,min(ObjV))
    line([gen-1, gen], [preObjV, min(ObjV)]);
    preObjV = min(ObjV);
    FitnV = Fitness(ObjV);
    %% 選擇
    SelCh = Select(Chrom,FitnV,GGAP);
    %% 交叉操作
    SelCh = Recombin(SelCh,Pc);
    %% 變異
    SelCh = Mutate(SelCh,Pm);
    %% 進化逆轉操作
    SelCh = Reverse(SelCh,D);
    %% 重插入子代的新種群
    Chrom = Reins(Chrom,SelCh,ObjV);
    %% 更新疊代次數
    gen = gen+1 ;
end
%% 畫出最優解的路線圖
ObjV = PathLength(D,Chrom);  %計算路線長度
[minObjV, minInd] = min(ObjV);
DrawPath(Chrom(minInd(1), :), X)
%% 輸出最優解的路線和總距離
disp('最優解:')
p = OutputPath(Chrom(minInd(1), :));
disp(['總距離：', num2str(ObjV(minInd(1)))]);
disp('-------------------------------------------------------------')
           

重插入子代得到新種群的函數Reins代碼如下：

function Chrom = Reins(Chrom, SelCh, ObjV)
%% 重插入子代的新種群
% 輸入：
% Chrom  父代的種群
% SelCh    子代種群
% ObjV     父代适應度
% 輸出：
% Chrom  組合父代與子代後得到的新種群
NIND = size(Chrom, 1);
NSel = size(SelCh, 1);
[TobjV, index] = sort(ObjV);
Chrom = [Chrom(index(1:NIND-NSel), :); SelCh];
           

(9).結果分析

優化前的一個随機路線軌迹圖如圖2所示。

圖2 随機路線圖

初始種群中的一個随機值:
10—>4—>8—>11—>9—>13—>3—>12—>14—>6—>1—>5—>2—>7—>10
總距離：70.3719
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
           

優化後的路線圖如圖3所示。

圖3 最優解路線圖

最優解:
5—>4—>3—>14—>2—>1—>10—>9—>11—>8—>13—>7—>12—>6—>5
總距離：29.3405
-------------------------------------------------------------
           

優化疊代圖如圖4所示。

圖4 遺傳算法進化過程圖

由優化圖可以看出，優化前後路徑長度得到很大改進，接近40代以後路徑長度已經保持不變了，可以認為是最優解了，總距離由優化前的70.3719變為29.3405，減為原來的41.7%。

具體細節可參考：

連結：https://pan.baidu.com/s/16vnho_Uf6O1JdoVx9L2enQ

提取碼：shgr

四、遺傳算法的改進

上述程式中，對遺傳算法做了以下兩處改進。

1. 使用精英政策

子代種群中的最優個體永遠不會比父代最優的個體差，這樣使得父代的好的個體不至于由于交叉或者變異操作而丢失。

2. 使用進化逆轉操作

同交叉算子相比較，逆轉算子能使子代繼承親代的較多資訊。

五、算法的局限性

當問題規模 n n n比較小時，得到的一般都是最優解；當規模比較大時，一般隻能得到近似解。這時可以通過增大種群大小和增加最大遺傳代數使得優化值更接近最優解。

六、參考文獻

[1] 儲理才. 基于MATLAB的遺傳算法程式設計及TSP問題求解[J]. 集美大學學報:自然科學版, 2001.

[2] 郁磊等. MATLAB智能算法30個案例分析(第2版)[M].北京航空航天大學出版社.2015年.