應(yīng)用GA和PSO算法求解10城市TSP問(wèn)題

1、精選優(yōu)質(zhì)文檔-傾情為你奉上應(yīng)用GA和PSO算法求解10城市TSP問(wèn)題1 問(wèn)題描述旅行團(tuán)計(jì)劃近期在城市A、B、C、D、E、F、G、H、I和J共10個(gè)城市間進(jìn)行一次周游旅行，為了盡量節(jié)省旅行開(kāi)支，希望能找到一條里程數(shù)最少或相對(duì)較少的旅行路線。問(wèn)題1，給定10個(gè)城市之間的公路里程如表1所示，并要求使用GA算法求解優(yōu)化問(wèn)題。問(wèn)題2，與問(wèn)題1數(shù)據(jù)相同，要求使用PSO算法求解優(yōu)化問(wèn)題。表1 城市位置坐標(biāo)（單位：km）橫坐標(biāo)縱坐標(biāo)城市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) 編碼和適應(yīng)度函數(shù)分別用1-10表示城市A-J，然后采用自然數(shù)編碼方式為T(mén)SP問(wèn)題編碼，例如，旅程（1 6 2 8 9 10 5 7 3 4）表示從城市A出發(fā)分別經(jīng)過(guò)了F-B-H-I-J-E-G-C-D的一次旅行。每一個(gè)問(wèn)題的解及算法中的個(gè)體都可以計(jì)算相應(yīng)的距離。那么種群中的最小距離和最大距離也相應(yīng)的可以確定。選擇種群個(gè)數(shù)為50。根據(jù)種群中個(gè)體的距離并考慮使用自適應(yīng)的標(biāo)定方法，定義如下的適應(yīng)度函數(shù)，使用此適應(yīng)度函數(shù)的后面的乘方次數(shù)可以調(diào)整來(lái)改變淘汰速度。這里選擇2，表示淘汰速度比較

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

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

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

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

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

8、度和位置。(d) 如果得到了更好的個(gè)體位置，更新。(iv) 如果得到了更好的群體位置，更新。3.2 PSO算法的計(jì)算結(jié)果編程實(shí)現(xiàn)3.1中的算法，計(jì)算結(jié)果如下。計(jì)算程序見(jiàn)附錄。圖3 PSO算法過(guò)程的距離值變化圖4 PSO算法求解的10城市TSP問(wèn)題的結(jié)果最后的結(jié)果編碼為（1 4 5 6 7 8 9 10 2 3），解為A-D-E-F-G-H-I-J-B-C，距離為269.0671。從計(jì)算結(jié)果可以看出，算法迭代到200步時(shí)已經(jīng)收斂到了局部的最優(yōu)值。這個(gè)局部最優(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; %計(jì)算適應(yīng)度函數(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);%隨機(jī)選擇交叉范圍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)計(jì)迭代次數(shù)。 counter=counter+1 ; end % 結(jié)果圖形顯示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 % 結(jié)果輸出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; % 粒子個(gè)數(shù)% 初始化個(gè)體和群體最優(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); % 更新每一個(gè)粒子的位置 for i=1:indiNumber % 個(gè)體最優(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 % 慣性項(xiàng) c1=round(rand*(n-1)+1; c2=