應用GA和PSO算法求解10城市TSP問題

1、精選優(yōu)質文檔-傾情為你奉上應用GA和PSO算法求解10城市TSP問題1 問題描述旅行團計劃近期在城市A、B、C、D、E、F、G、H、I和J共10個城市間進行一次周游旅行，為了盡量節(jié)省旅行開支，希望能找到一條里程數(shù)最少或相對較少的旅行路線。問題1，給定10個城市之間的公路里程如表1所示，并要求使用GA算法求解優(yōu)化問題。問題2，與問題1數(shù)據(jù)相同，要求使用PSO算法求解優(yōu)化問題。表1 城市位置坐標（單位：km）橫坐標縱坐標城市A4044.39城市B24.3914.63城市C17.0722.93城市D22.9376.1城市E51.7194.14城市F87.3265.36城市G68.7852.19城市H

2、84.8836.09城市I66.8325.36城市J61.9526.342 使用GA算法求解2.1 算法描述(1) 編碼和適應度函數(shù)分別用1-10表示城市A-J，然后采用自然數(shù)編碼方式為TSP問題編碼，例如，旅程（1 6 2 8 9 10 5 7 3 4）表示從城市A出發(fā)分別經(jīng)過了F-B-H-I-J-E-G-C-D的一次旅行。每一個問題的解及算法中的個體都可以計算相應的距離。那么種群中的最小距離和最大距離也相應的可以確定。選擇種群個數(shù)為50。根據(jù)種群中個體的距離并考慮使用自適應的標定方法，定義如下的適應度函數(shù)，使用此適應度函數(shù)的后面的乘方次數(shù)可以調整來改變淘汰速度。這里選擇2，表示淘汰速度比較

3、適中。(2) 定義算子選擇算子，根據(jù)適應度函數(shù)的大小進行選擇。計算每個個體被選中的概率，以各個個體所分配到的概率值作為其遺傳到下一代的概率，基于這些概率用賭盤選擇法來產(chǎn)生下一代群體。交叉算子，采用部分映射交叉(Partially Mapped Crossover, PMX)方法實現(xiàn)算法交叉。首先選取選需要交叉的區(qū)間段，然后確定區(qū)間段的映射關系，接下來交換區(qū)間段的遺傳信息，此時未換部分的遺傳信息會出現(xiàn)不合法的情況，因此根據(jù)將未換部分按映射關系進行交換。交叉率為0.9。變異算子，把一個染色體中的兩個基因的交換作為變異算法。在算法中隨機的產(chǎn)生一個染色體中需要交換的兩個基因的位置，將這兩個基因進行交換

4、來實現(xiàn)變異。變異率為0.2。(3) 算法流程根據(jù)上述的遺傳算子的定義，并設置最大的迭代次數(shù)為1000，將GA算法流程敘述如下，(i) 隨機生成初始種群。(ii) 按照本節(jié)(1)中的公式計算各個個體適應度的值。(iii) 判斷是否達到了最大的迭代次數(shù)。若是，算法結束，輸出計算結果；若否，轉到(iv)。(iv) 按照本節(jié)(2)中的選擇算子進行選擇操作。(v) 選擇交叉寬度并隨機確定交叉切點，按照本節(jié)(2)中的交叉算子進行交叉操作。(vi) 按照本節(jié)(2)中的變異算子進行變異操作。(vii) 將遺傳算子產(chǎn)生的種群作為新的種群，轉到(ii)。2.2 GA算法計算結果使用Matlab編程實現(xiàn)2.1中的算

5、法，計算結果如下。圖1 GA算法過程的距離值變化圖2 GA算法求解的10城市TSP問題的結果最后的結果編碼為（8 9 10 2 3 1 4 5 6 7），解為H-I-J-B-C-A-D-E-F-G，距離為269.0671。從計算結果可以看出，算法迭代到300步時已經(jīng)收斂到了局部的最優(yōu)值。經(jīng)過大量的計算發(fā)現(xiàn)至多迭代到300步，算法收斂到局部最優(yōu)值。經(jīng)過進一步的驗證發(fā)現(xiàn)，這個局部最優(yōu)解也是全局最優(yōu)解。3 使用PSO算法求解3.1 算法描述(1)TSP問題的交換子和交換序設n個節(jié)點的TSP問題的解序列為s=(ai)，i=1n。定義交換子SO(i1,i2)為交換解S中的點ai1和ai2，則S=S+SO

6、(i1,i2)為解S經(jīng)算子SO(i1,i2)操作后的新解。 一個或多個交換子的有序隊列就是交換序，記作SS，SS=(SO1,SO2,SON)。其中，SO1,SON等是交換子，之間的順序是有意義的, 意味著所有的交換子依次作用于某解上。若干個交換序可以合并成一個新的交換序，定義為兩個交換序的合并算子。設兩個交換序SS1和SS2按先后順序作用于解S上，得到新解S。假設另外有一個交換序SS作用于同一解S上，能夠得到形同的解S，可定義和屬于同一等價集，在交換序等價集中，擁有最少交換子的交換序稱為該等價集的基本交換序。(2) TSP問題的速度和位置更新算式根據(jù)(1)中的交換子和交換序的定義，可以根據(jù)基本

7、的PSO算法速度更新算式和位置更新算式，重新定義如下的速度和位置更新算式，式中，、為0,1區(qū)間的隨機數(shù)。為粒子與個體極值的交換序，以概率保留；為粒子與全局極值的交換序，以概率保留。粒子的位置按照交換序進行更新。為慣性權重。(3) 算法流程按照本節(jié)中的有關定義給出算法流程如下，(i) 初始化粒子群，給每一個粒子一個初始解和隨機的交換序。(ii) 判斷是否達到最大迭代次數(shù)1000。若是，算法結束，輸出結果；若不是，轉到(iii)。(iii) 根據(jù)粒子當前位置計算下一個新解：(a) 計算A=，A是一個基本交換序，表示A作用于得到；(b) 計算B=，B是一個基本交換序；(c) 按照(2)中的公式更新速

8、度和位置。(d) 如果得到了更好的個體位置，更新。(iv) 如果得到了更好的群體位置，更新。3.2 PSO算法的計算結果編程實現(xiàn)3.1中的算法，計算結果如下。計算程序見附錄。圖3 PSO算法過程的距離值變化圖4 PSO算法求解的10城市TSP問題的結果最后的結果編碼為（1 4 5 6 7 8 9 10 2 3），解為A-D-E-F-G-H-I-J-B-C，距離為269.0671。從計算結果可以看出，算法迭代到200步時已經(jīng)收斂到了局部的最優(yōu)值。這個局部最優(yōu)解也是全局最優(yōu)解。從收斂的速度的平均意義上而言，PSO算法與GA算法比收斂的更快。附錄% GA算法的代碼% GA算法程序. data = l

9、oad('pos10.txt');a=data(:,2) data(:,3); n=50; % 種群數(shù)目C=1000; % 迭代次數(shù) Pc=0.9; % 交叉概率Pm=0.2; % 變異概率 % GA算法初始化. N,NN=size(D); farm=zeros(n,N); for i=1:n farm(i,:)=randperm(N); end R=farm(1,:); len=zeros(n,1); fitness=zeros(n,1); counter=0; % GA算法迭代 while counter<C for i=1:n len(i,1)=myLength(D

10、,farm(i,:); end maxlen=max(len); minlen=min(len); fitness=len;for i=1:length(len) fitness(i,1)=(1-(len(i,1)-minlen)/(maxlen-minlen+0.0001).2;end; rr=find(len=minlen); R=farm(rr(1,1),:); FARM=farm; %計算適應度函數(shù) nn=0; for i=1:n if fitness(i,1)>=rand nn=nn+1; FARM(nn,:)=farm(i,:); end end FARM=FARM(1:nn

11、,:); % FARM (nn*N) aa,bb=size(FARM);%(aa=nn) % 交叉 FARM2=FARM; for i=1:2:aa if Pc>rand&&i<aa A=FARM(i,:); B=FARM(i+1,:); L=length(a);if L<=10 %確定交叉寬度 W=9;elseif (L/10)-floor(L/10)>=rand&&L>10 W=ceil(L/10)+8;else W=floor(L/10)+8;endp=unidrnd(L-W+1);%隨機選擇交叉范圍for i=1:W x=f

12、ind(a=b(1,p+i-1); y=find(b=a(1,p+i-1); a(1,p+i-1),b(1,p+i-1)=exchange(a(1,p+i-1),b(1,p+i-1); a(1,x),b(1,y)=exchange(a(1,x),b(1,y); end FARM(i,:)=A; FARM(i+1,:)=B; end end % 變異 FARM2=FARM; for i=1:aa if Pm>=rand FARM(i,:)=mutate(FARM(i,:); end end % 群體更新 FARM2=zeros(n-aa+1,N); if n-aa>=1 for i=

13、1:n-aa FARM2(i,:)=randperm(N); end end FARM=R;FARM;FARM2; aa,bb=size(FARM); if aa>n FARM=FARM(1:n,:); end farm=FARM; clear FARM RR(counter+1)=myLength(D,R)% 統(tǒng)計迭代次數(shù)。 counter=counter+1 ; end % 結果圖形顯示figurehold onplot(a(R(1),1),a(R(10),1),a(R(1),2),a(R(10),2),'ms-','LineWidth',2,'

14、;MarkerEdgeColor','k','MarkerFaceColor','g');scatter(a(:,1),a(:,2),'ms')grid onhold onfor i=2:length(R) x0=a(R(i-1),1); y0=a(R(i-1),2); x1=a(R(i),1); y1=a(R(i),2); xx=x0,x1; yy=y0,y1; plot(xx,yy,'LineWidth',4,'MarkerEdgeColor','k','Mark

15、erFaceColor','g') hold onend % 結果輸出RRlength% 其他函數(shù)function a=mutate(a)L=length(a);rray=randperm(L);a(rray(1),a(rray(2)=exchange(a(rray(1),a(rray(2);function x,y=exchange(x,y)temp=x;x=y;y=temp;function len=myLength(D,p)N,NN=size(D);len=D(p(1,N),p(1,1);for i=1:(N-1) len=len+D(p(1,i),p(1,i+1

16、);end% PSO算法代碼% PSO算法代碼data=load('pos10.txt')cityCoor=data(:,2) data(:,3); n=size(cityCoor,1); cityDist=zeros(n,n); for i=1:n for j=1:n if i=j cityDist(i,j)=(cityCoor(i,1)-cityCoor(j,1)2+. (cityCoor(i,2)-cityCoor(j,2)2)0.5; end cityDist(j,i)=cityDist(i,j); endendindividual=zeros(indiNumber,n

17、);% 初始化nMax=1000; % 迭代次數(shù)indiNumber=50; % 粒子個數(shù)% 初始化個體和群體最優(yōu)值indiFit=fitness(individual,cityCoor,cityDist);value,index=min(indiFit);tourPbest=individual; tourGbest=individual(index,:) ; recordPbest=inf*ones(1,indiNumber); recordGbest=indiFit(index); xnew1=individual;% 循環(huán)尋找最優(yōu)路徑L_best=zeros(1,nMax);for N

18、=1:nMax % 更新最優(yōu)值 indiFit=fitness(individual,cityCoor,cityDist); for i=1:indiNumber if indiFit(i)<recordPbest(i) recordPbest(i)=indiFit(i); tourPbest(i,:)=individual(i,:); end if indiFit(i)<recordGbest recordGbest=indiFit(i); tourGbest=individual(i,:); end end value,index=min(recordPbest); recor

19、dGbest(N)=recordPbest(index); % 更新每一個粒子的位置 for i=1:indiNumber % 個體最優(yōu)值更新速度 c1=unidrnd(n-1); c2=unidrnd(n-1); while c1=c2 c1=round(rand*(n-2)+1; c2=round(rand*(n-2)+1; end chb1=min(c1,c2); chb2=max(c1,c2); cros=tourPbest(i,chb1:chb2); ncros=size(cros,2); for j=1:ncros for k=1:n if xnew1(i,k)=cros(j) x

20、new1(i,k)=0; for t=1:n-k temp=xnew1(i,k+t-1); xnew1(i,k+t-1)=xnew1(i,k+t); xnew1(i,k+t)=temp; end end end end xnew1(i,n-ncros+1:n)=cros; dist=0; for j=1:n-1 dist=dist+cityDist(xnew1(i,j),xnew1(i,j+1); end dist=dist+cityDist(xnew1(i,1),xnew1(i,n); if indiFit(i)>dist individual(i,:)=xnew1(i,:); end

21、 % 全體最優(yōu)值更新速度 c1=round(rand*(n-2)+1; c2=round(rand*(n-2)+1; while c1=c2 c1=round(rand*(n-2)+1; c2=round(rand*(n-2)+1; end chb1=min(c1,c2); chb2=max(c1,c2); cros=tourGbest(chb1:chb2); ncros=size(cros,2); for j=1:ncros for k=1:n if xnew1(i,k)=cros(j) xnew1(i,k)=0; for t=1:n-k temp=xnew1(i,k+t-1); xnew1(i,k+t-1)=xnew1(i,k+t); xnew1(i,k+t)=temp; end end end end xnew1(i,n-ncros+1:n)=cros; dist=0; for j=1:n-1 dist=dist+cityDist(xnew1(i,j),xnew1(i,j+1); end dist=dist+cityDist(xnew1(i,1),xnew1(i,n); if indiFit(i)>dist individual(i,:)=xnew1(i,:); end % 慣性項 c1=round(rand*(n-1)+1; c2=