遺傳演算法 求解旅行商 TSP 問題，matlab程式碼
阿新 • • 發佈：2019-02-11
學習啟發式演算法時，旅行商問題是一個經典的例子。其中，遺傳演算法可以用來求解該問題。遺傳演算法是一種進化演算法，由於其啟發式演算法的屬性，並不能保證得到最優解。求解效果與初始種群選取，編碼方法，選擇方法，交叉變異規則有關。
上課時，老師不知從哪裡找了一個非常粗糙的程式，自己將不少錯誤修正，增加了一些註釋方便理解，並增加了一些程式碼使程式更完美。該程式碼可以動態顯示每一代的路線，非常直觀！
顯示效果：
程式碼：
function ga_TSP % mainly amended by Chen Zhen, 2012~2016 CityNum=30; %you chan choose 10, 30, 50, 75 [dislist,Clist]=tsp(CityNum); inn=30; %初始種群大小 gnmax=500; %最大代數 pc=0.8; %交叉概率 pm=0.8; %變異概率 %產生初始種群 s=zeros(inn,CityNum); for i=1:inn s(i,:)=randperm(CityNum); end [~,p]=objf(s,dislist); gn=1; ymean=zeros(gn,1); ymax=zeros(gn,1); xmax=zeros(inn,CityNum); scnew=zeros(inn,CityNum); smnew=zeros(inn,CityNum); while gn<gnmax+1 for j=1:2:inn seln=sel(p); %選擇操作 scro=cro(s,seln,pc); %交叉操作 scnew(j,:)=scro(1,:); scnew(j+1,:)=scro(2,:); smnew(j,:)=mut(scnew(j,:),pm); %變異操作 smnew(j+1,:)=mut(scnew(j+1,:),pm); end s=smnew; %產生了新的種群 [f,p]=objf(s,dislist); %計算新種群的適應度 %記錄當前代最好和平均的適應度 [fmax,nmax]=max(f); ymean(gn)=1000/mean(f); ymax(gn)=1000/fmax; %記錄當前代的最佳個體 x=s(nmax,:); xmax(gn,:)=x; drawTSP(Clist,x,ymax(gn),gn,0); gn=gn+1; end [min_ymax,index]=min(ymax); drawTSP(Clist,xmax(index,:),min_ymax,index,1); figure(2); plot(ymax,'r'); hold on; plot(ymean,'b');grid; title('搜尋過程'); legend('最優解','平均解'); fprintf('遺傳演算法得到的最短距離:%.2f\n',min_ymax); fprintf('遺傳演算法得到的最短路線'); disp(xmax(index,:)); end %------------------------------------------------ %計算所有種群的適應度 function [f,p]=objf(s,dislist) inn=size(s,1); %讀取種群大小 f=zeros(inn,1); for i=1:inn f(i)=CalDist(dislist,s(i,:)); %計算函式值，即適應度 end f=1000./f'; %取距離倒數 %根據個體的適應度計算其被選擇的概率 fsum=0; for i=1:inn fsum=fsum+f(i)^15;% 讓適應度越好的個體被選擇概率越高 end ps=zeros(inn,1); for i=1:inn ps(i)=f(i)^15/fsum; end %計算累積概率 p=zeros(inn,1); p(1)=ps(1); for i=2:inn p(i)=p(i-1)+ps(i); end p=p'; end %-------------------------------------------------- %根據變異概率判斷是否變異 function pcc=pro(pc) test(1:100)=0; l=round(100*pc); test(1:l)=1; n=round(rand*99)+1; pcc=test(n); end %-------------------------------------------------- %“選擇”操作 function seln=sel(p) seln=zeros(2,1); %從種群中選擇兩個個體，最好不要兩次選擇同一個個體 for i=1:2 r=rand; %產生一個隨機數 prand=p-r; j=1; while prand(j)<0 j=j+1; end seln(i)=j; %選中個體的序號 if i==2&&j==seln(i-1) %%若相同就再選一次 r=rand; %產生一個隨機數 prand=p-r; j=1; while prand(j)<0 j=j+1; end end end end %------------------------------------------------ %“交叉”操作 function scro=cro(s,seln,pc) bn=size(s,2); pcc=pro(pc); %根據交叉概率決定是否進行交叉操作，1則是，0則否 scro(1,:)=s(seln(1),:); scro(2,:)=s(seln(2),:); if pcc==1 c1=round(rand*(bn-2))+1; %在[1,bn-1]範圍內隨機產生一個交叉位 c2=round(rand*(bn-2))+1; chb1=min(c1,c2); chb2=max(c1,c2); middle=scro(1,chb1+1:chb2); scro(1,chb1+1:chb2)=scro(2,chb1+1:chb2); scro(2,chb1+1:chb2)=middle; for i=1:chb1 %似乎有問題 while find(scro(1,chb1+1:chb2)==scro(1,i)) zhi=find(scro(1,chb1+1:chb2)==scro(1,i)); y=scro(2,chb1+zhi); scro(1,i)=y; end while find(scro(2,chb1+1:chb2)==scro(2,i)) zhi=find(scro(2,chb1+1:chb2)==scro(2,i)); y=scro(1,chb1+zhi); scro(2,i)=y; end end for i=chb2+1:bn while find(scro(1,1:chb2)==scro(1,i)) zhi=logical(scro(1,1:chb2)==scro(1,i)); y=scro(2,zhi); scro(1,i)=y; end while find(scro(2,1:chb2)==scro(2,i)) zhi=logical(scro(2,1:chb2)==scro(2,i)); y=scro(1,zhi); scro(2,i)=y; end end end end %-------------------------------------------------- %“變異”操作 function snnew=mut(snew,pm) bn=size(snew,2); snnew=snew; pmm=pro(pm); %根據變異概率決定是否進行變異操作，1則是，0則否 if pmm==1 c1=round(rand*(bn-2))+1; %在[1,bn-1]範圍內隨機產生一個變異位 c2=round(rand*(bn-2))+1; chb1=min(c1,c2); chb2=max(c1,c2); x=snew(chb1+1:chb2); snnew(chb1+1:chb2)=fliplr(x); end end %------------------------------------------------ %城市位置座標 function [DLn,cityn]=tsp(n) DLn=zeros(n,n); if n==10 city10=[0.4 0.4439;0.2439 0.1463;0.1707 0.2293;0.2293 0.761;0.5171 0.9414; 0.8732 0.6536;0.6878 0.5219;0.8488 0.3609;0.6683 0.2536;0.6195 0.2634];%10 cities d'=2.691 for i=1:10 for j=1:10 DLn(i,j)=((city10(i,1)-city10(j,1))^2+(city10(i,2)-city10(j,2))^2)^0.5; end end cityn=city10; end if n==30 city30=[41 94;37 84;54 67;25 62;7 64;2 99;68 58;71 44;54 62;83 69;64 60;18 54;22 60; 83 46;91 38;25 38;24 42;58 69;71 71;74 78;87 76;18 40;13 40;82 7;62 32;58 35;45 21;41 26;44 35;4 50];%30 cities d'=423.741 by D B Fogel for i=1:30 for j=1:30 DLn(i,j)=((city30(i,1)-city30(j,1))^2+(city30(i,2)-city30(j,2))^2)^0.5; end end cityn=city30; end if n==50 city50=[31 32;32 39;40 30;37 69;27 68;37 52;38 46;31 62;30 48;21 47;25 55;16 57; 17 63;42 41;17 33;25 32;5 64;8 52;12 42;7 38;5 25; 10 77;45 35;42 57;32 22; 27 23;56 37;52 41;49 49;58 48;57 58;39 10;46 10;59 15;51 21;48 28;52 33; 58 27;61 33;62 63;20 26;5 6;13 13;21 10;30 15;36 16;62 42;63 69;52 64;43 67];%50 cities d'=427.855 by D B Fogel for i=1:50 for j=1:50 DLn(i,j)=((city50(i,1)-city50(j,1))^2+(city50(i,2)-city50(j,2))^2)^0.5; end end cityn=city50; end if n==75 city75=[48 21;52 26;55 50;50 50;41 46;51 42;55 45;38 33;33 34;45 35;40 37;50 30; 55 34;54 38;26 13;15 5;21 48;29 39;33 44;15 19;16 19;12 17;50 40;22 53;21 36; 20 30;26 29;40 20;36 26;62 48;67 41;62 35;65 27;62 24;55 20;35 51;30 50; 45 42;21 45;36 6;6 25;11 28;26 59;30 60;22 22;27 24;30 20;35 16;54 10;50 15; 44 13;35 60;40 60;40 66;31 76;47 66;50 70;57 72;55 65;2 38;7 43;9 56;15 56; 10 70;17 64;55 57;62 57;70 64;64 4;59 5;50 4;60 15;66 14;66 8;43 26];%75 cities d'=549.18 by D B Fogel for i=1:75 for j=1:75 DLn(i,j)=((city75(i,1)-city75(j,1))^2+(city75(i,2)-city75(j,2))^2)^0.5; end end cityn=city75; end end %------------------------------------------------ %適應度函式 function F=CalDist(dislist,s) DistanV=0; n=size(s,2); for i=1:(n-1) DistanV=DistanV+dislist(s(i),s(i+1)); end DistanV=DistanV+dislist(s(n),s(1)); F=DistanV; end %------------------------------------------------ %畫圖 function drawTSP(Clist,BSF,bsf,p,f) CityNum=size(Clist,1); for i=1:CityNum-1 plot([Clist(BSF(i),1),Clist(BSF(i+1),1)],[Clist(BSF(i),2),Clist(BSF(i+1),2)],'ms-','LineWidth',2,'MarkerEdgeColor','k','MarkerFaceColor','g'); text(Clist(BSF(i),1),Clist(BSF(i),2),[' ',int2str(BSF(i))]); text(Clist(BSF(i+1),1),Clist(BSF(i+1),2),[' ',int2str(BSF(i+1))]); hold on; end plot([Clist(BSF(CityNum),1),Clist(BSF(1),1)],[Clist(BSF(CityNum),2),Clist(BSF(1),2)],'ms-','LineWidth',2,'MarkerEdgeColor','k','MarkerFaceColor','g'); title([num2str(CityNum),'城市TSP']); if f==0&&CityNum~=10 text(5,5,['第 ',int2str(p),' 代',' 最短距離為 ',num2str(bsf)]); else text(5,5,['最終搜尋結果：最短距離 ',num2str(bsf),'， 在第 ',num2str(p),' 代達到']); end if CityNum==10 if f==0 text(0,0,['第 ',int2str(p),' 代',' 最短距離為 ',num2str(bsf)]); else text(0,0,['最終搜尋結果：最短距離 ',num2str(bsf),'， 在第 ',num2str(p),' 代達到']); end end hold off; pause(0.05); end