一隻菜雞的遺傳演算法求解TSP問題嘗試

身為一個剛入坑不久的菜雞，剛剛學運籌學不久，自從看了《迷茫的旅行商》[1]後對TSP問題的興趣都很大，雖然網上有很多解法，也有很多前輩做出來，深知在逼乎中各路大佬都有解決過，但還是忍不住想把自己幾天的成果記錄下來，作為一個小小的紀念。各位大佬見到莫笑。

這次求解的一個TSP問題是一個只有十個城市的小型的TSP問題，對於新手來說比較友好，大佬可以直接跳過。主要是為像我一樣剛剛入坑的小白提供一個入門的基礎，讓小白們對遺傳演算法求解TSP問題有一個比較基礎的認識。其中如果有一些錯誤的地方，還請各位指出。

話不多說，我們直接進入正題。這次我們選取的題目如下：

題目的數據來自百度文庫[2]，是一個已經解決的問題的數據，只是題目的原話換了而已。在實際求解之前，我們先來回顧一下遺傳演算法的步驟：

• 種群中的每個個體表示一個旅行方案 ，每個旅行方案用1~N的某種排列表示，用來表示一次TSP環行的路線順序。
• 每個方案的總費用為
• 適應度函數取
• 採用最優個體遺傳策略（我也不記得具體叫啥了反正就是把最好的那個個體遺傳到下一代）
• 輪盤賭策略的意思就是隨機生成一個0到1的數，像一個骰子一樣，然後選擇最小的比這個數大的累計概率的個體進行交叉操作。
• 交叉操作就是在選取的父代個體中隨機選取交叉位點，交換這兩個個體交叉位點之間的序列，得到兩個序列，再通過衝突調整，得到新的子代個體，然後加入到新的子代種群中。我的程序中為了控制種群的大小，選擇把兩個子代個體中更好的那個加入到新子代中，而較差的那個則捨棄不用。重複若干次，直到達到流程圖中的要求為止。
• 衝突調整我參考了[3]中的衝突調整方法。

寫出的MATLAB程序代碼如下（演算法的一些具體細節在程序的注釋中都已寫清）：

注意：這裡的程序代碼是好幾個m文件的程序放一起了的，一定要注意分別建立相應的m文件之後再運行main.m（已經有好多個小夥伴還沒弄懂怎麼用matlab就直接複製上去了，然後就出錯……每個m文件前邊我都有注釋，注釋的格式是：%xxx.m其中xxx是m文件名）

%sumlong.mfunction s=sumlong(x,d)%這個函數是求方案所走總距離的函數，返回的是該方案對應的位點序列%x是方案對應的位點序列%waylong是各點間的距離s=d(x(1),x(end));for i=1:length(x)-1 s=s+d(x(i),x(i+1));%main.m%這是遺傳演算法的主程序clcclear% d=[0 6 15 19 23 9 11 20 24 10% 6 0 27 30 32 15 16 23 29 18% 15 27 0 4 7 20 21 33 40 19% 19 30 4 0 3 17 15 39 35 25% 23 32 7 3 0 22 23 41 37 28% 9 15 20 17 22 0 2 18 26 9% 11 16 21 15 23 2 0 15 24 10% 20 23 33 39 41 18 15 0 6 23% 24 29 40 35 37 26 24 6 0 30% 10 18 19 25 28 9 10 23 30 0];%以上是另一個問題的距離矩陣d=[0 7 4 5 8 6 12 13 11 18 7 0 3 10 9 14 5 14 17 17 4 3 0 5 9 10 21 8 27 12 5 10 5 0 14 9 10 9 23 16 8 9 9 14 0 7 8 7 20 19 6 14 10 9 7 0 13 5 25 13 12 5 21 10 8 13 0 23 21 18 13 14 8 9 7 5 23 0 18 12 11 17 27 23 20 25 21 18 0 16 18 17 12 16 19 13 18 12 16 0];n=input(種群大小n=); %確定初始種群大小age=input(迭代次數=); %確定遺傳代數r=input(交叉率=); %交叉率v=input(變異率=); %變異率L=length(d(1,:)); %獲得城市數s=[]; %初始化種群矩陣for i=1:L d(i,i)=1000000000000000000; %這個操作可要可不要，作用是防止取到無效的路徑endfor i=1:n s=[s;randperm(L)]; %隨機產生n個個體，形成初始種群 end% % for i=1:n% [s(i,:),~]=modifycircle(s(i,:),L,d);% end[numway,wayside]=size(s); %獲得種群個數和城市數，其實有點重複可以改一下，但是我懶得改了，如果你要改的話請把之後相應的變數也要做好處理，建議是不要改，因為它們對結果沒有影響mins=[]; %初始化最優值矩陣，這個是用來記錄每次迭代生成的種群中最優值的for i=1:age [s,mins(i)]=GABY2(s,d,r,v); %對種群不斷進行迭代，同時記錄每次迭代出來的最優值endfitness=[];for i=1:numway fitness=[fitness;1/sumlong(s(i,:),d)]; %計算最終得到的種群中各個個體的適應度end[maxfitness,indexmax]=max(fitness); %找出最優個體的位置和相應的最大適應度%以下打注釋部分的語句是用來繪製路徑圖的，如果你的TSP問題是使用各點之間的二維坐標表示，且相應的兩點間的距離使用兩點坐標決定的歐氏距離的話，可以使用下面的語句生成最優解的路徑圖% figure;% l=length(s(indexmax,:));% plot([x(s(indexmax,1)),x(s(indexmax,end))],[y(s(indexmax,1)),y(s(indexmax,end))]);% hold on;% for j=1:l-1% plot([x(s(indexmax,j)),x(s(indexmax,j+1))],[y(s(indexmax,j)),y(s(indexmax,j+1))]);% hold on;% end% xlabel(x);% ylabel(y);% title(路徑);%輸出最優路徑及相應的耗費disp(最優路徑為:);best_way=s(indexmax,:)disp(最優路徑的耗費：);min_cost=sumlong(s(indexmax,:),d)%繪製最優值變化趨勢圖figure; plot(1:age,mins)xlabel(迭代次數);ylabel(每次迭代的種群中的最少耗費);title(最優值變化趨勢)%這個jiancha只是我當初為了檢查結果正確與否而設置的，因為調整的時候老是出錯jiancha=sort(s(indexmax,:))%GABY2.m%這是遺傳演算法的迭代函數function [A,mins]=GABY2(s,d,r,v)%s是種群樣本矩陣%d是兩點間距離矩陣%r是每次迭代中通過交叉取代群體成員的比例，建議取值較大，在0.9以上a1=[];[m,t]=size(s);fitness=[];for i=1:mfitness=[fitness;1/sumlong(s(i,:),d)]; %計算各個個體的適應度endf=sum(fitness(:));p=fitness./f; %各方案被選中的概率[~,maxA]=max(p);mins=sumlong(s(maxA,:),d); %記錄當前種群中最優個體的費用sump=[];for i=1:length(p)sump=[sump,sum(p(1:i))]; %累積概率endn=round((1-r)*m); %被保留下來的最優樣本數[~,l]=sort(p);a1=[a1;s(l(length(p)-n+1:end),:)]; %適應度最高的(1-r)n個樣本被保留傳遞至下一代tnews=s(l(1:m-n),:);[g,numcity]=size(tnews); %從剩下的個體中使用輪盤賭策略挑選父代個體進行交叉變異操作i=1;a=[];while i<=gfm=rand(1,2); %隨機產生兩個0-1間的數，作為累積概率的取值point=sort(round(1+(numcity-3)*rand(1,2))); %隨機選取交叉基因位點 while point(1)==point(2)%|| fpoint==1 || mpoint==1 || fpoint==k ||mpoint==k point=sort(round(1+(numcity-3)*rand(1,2))); %防止交叉位點重複 end father=[]; mother=[]; fathernum=1; mothernum=1; while fm(1)>sump(fathernum) fathernum=fathernum+1; %根據隨機產生的概率選擇父代個體 end father=s(fathernum,:); while fm(2)>sump(mothernum) mothernum=mothernum+1; %根據隨機產生的概率選擇父代個體 end mother=s(mothernum,:);%這個while循環是防止產生的父親和母親重複 while mother==father fm=rand(1,2); father=[]; mother=[]; fathernum=1; mothernum=1; while fm(1)>sump(fathernum) fathernum=fathernum+1; end father=s(fathernum,:); while fm(2)>sump(mothernum) mothernum=mothernum+1; end mother=s(mothernum,:); end %初始化子代個體 children1=[]; children2=[]; flag1=1; flag2=1; %交叉操作，間交叉位點中的基因進行調換 children1=[father(1,1:point(1)),mother(1,(point(1)+1):(point(2)-1)),father(1,point(2):end)]; children2=[mother(1,1:point(1)),father(1,(point(1)+1):(point(2)-1)),mother(1,point(2):end)]; %以下兩個while循環是為了處理產生的兩個子代個體的基因，作用是防止產生的子代個體中的基因出現重複 while flag1>0 flag1=0; for j=1:length(children1) for k=1:length(children1) if j~=k && children1(j)==children1(k) if children1(j)~=children2(j) x=children1(k); children1(k)=children2(j); children2(j)=x; flag1=1; else x=children1(j); children1(j)=children2(k); children2(k)=x; flag1=1; end end end end end while flag2>0 flag2=0; for j=1:length(children2) for k=1:length(children2) if j~=k && children2(j)==children2(k) if children1(j)~=children2(j) x=children2(k); children2(k)=children1(j); children1(j)=x; flag2=1; else x=children2(j); children2(j)=children1(k); children1(k)=x; flag2=1; end end end end end %選擇適應度更高的個體作為子代個體加入子代中 if sumlong(children1,d)<=sumlong(children2,d) a=[a;children1]; else a=[a;children2]; end i=i+1;end%在生成的子代中，隨機選取vn個個體進行變異操作xc=round((g-1)*rand(1,fix(m*v))+1);for i=1:length(xc)% [a(xc(i),:),~]=modifycircle(a(xc(i),:),t,d);%這是使用LK改圈法對所需變異個體進行變異操作 a(xc(i),:)=randperm(t);%這是重新打亂子代的基因順序，從而生成一個新的子代enda1=[a1;a];%輸出迭代後的子代種群A=a1;

調用程序後得到的結果：

從結果可以看出，得到的結果已經達到了線性規劃模型中所求的的最優解的下界，應視為求出最優解。所以遺傳演算法是一個比較不錯的求解TSP問題的演算法。但這也許不能看出遺傳演算法收斂的速度，那我們繼續看下面的圖。

從上圖可以看出，遺傳演算法大概在十幾次迭代時就已經收斂到最優解了，可見，遺傳演算法擁有很強的收斂性。對於大型的TSP問題，遺傳演算法有著先天的優勢。因此，在求解大型的TSP問題時，推薦使用遺傳演算法這樣的自啟發式搜索演算法。

得到的最優旅行路線為：

參考文獻

[1] 庫克(Cook,W.J.)(美)著;隋春寧(譯).迷茫的旅行商:一個無處不在的計算機演算法問題[M].北京:人民郵電出版社,2013.10

[2] 百度文庫. TSP問題及LINGO求解技巧[EB/OL].

TSP問題及LINGO求解技巧_百度文庫

[3] mylovestart.遺傳演算法解決TSP問題[EB/OL].

遺傳演算法解決TSP問題 - CSDN博客

[4] 黃厚生.求解旅行商問題的新方法研究[D].天津:天津大學，2005

[5] 米歇爾(Mitchell,T.M.)著；曾華軍等譯.機器學習[M].北京:機械工業出版社，2003.1（2016.4重印）

資料共享

本次所有的程序文件都已共享至百度網盤，僅供學習交流用，轉載請註明出處，未經允許不得用於商業運營，本人將保留法律追訴權。

百度網盤地址：

https://pan.baidu.com/s/1bpm8S8J 密碼：vfhl

謝謝大家，歡迎指出錯誤糾正。