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

應(yīng)用GA和PSO算法求解10城市TSP問題1問題描述旅行團(tuán)方案近期在城市A、B、C、D、E、F、G、H、I和J共10個城市間進(jìn)行一次周游旅行，為了盡量節(jié)省旅行開支，希望能找到一條里程數(shù)最少或相對較少的旅行路線。問題1，給定10個城市之間的公路里程如表1所示，并要求使用GA算法求解優(yōu)化問題。問題2，與問題1數(shù)據(jù)相同，要求使用PSO算法求解優(yōu)化問題。表1城市位置坐標(biāo)〔單位：km〕橫坐標(biāo)縱坐標(biāo)城市A40 44.39城市B24.3914.63城市C17.0722.93城市D22.9376.1城市E51.7194.14城市F87.3265.36城市G68.7852.19城市H84.8836.09城市I66.8325.36城市J61.9526.342使用GA算法求解2.1算法描述 (1)編碼和適應(yīng)度函數(shù) 分別用1-10表示城市A-J，然后采用自然數(shù)編碼方式為TSP問題編碼，例如，旅程〔16289105734〕表示從城市A出發(fā)分別經(jīng)過了F-B-H-I-J-E-G-C-D的一次旅行。每一個問題的解及算法中的個體都可以計算相應(yīng)的距離。那么種群中的最小距離和最大距離也相應(yīng)的可以確定。選擇種群個數(shù)為50。根據(jù)種群中個體的距離并考慮使用自適應(yīng)的標(biāo)定方法，定義如下的適應(yīng)度函數(shù)，使用此適應(yīng)度函數(shù)的后面的乘方次數(shù)可以調(diào)整來改變淘汰速度。這里選擇2，表示淘汰速度比擬適中。 (2)定義算子選擇算子，根據(jù)適應(yīng)度函數(shù)的大小進(jìn)行選擇。計算每個個體被選中的概率，以各個個體所分配到的概率值作為其遺傳到下一代的概率，基于這些概率用賭盤選擇法來產(chǎn)生下一代群體。 交叉算子，采用局部映射交叉(PartiallyMappedCrossover,PMX)方法實(shí)現(xiàn)算法交叉。首先選取選需要交叉的區(qū)間段，然后確定區(qū)間段的映射關(guān)系，接下來交換區(qū)間段的遺傳信息，此時未換局部的遺傳信息會出現(xiàn)不合法的情況，因此根據(jù)將未換局部按映射關(guān)系進(jìn)行交換。交叉率為0.9。 變異算子，把一個染色體中的兩個基因的交換作為變異算法。在算法中隨機(jī)的產(chǎn)生一個染色體中需要交換的兩個基因的位置，將這兩個基因進(jìn)行交換來實(shí)現(xiàn)變異。變異率為0.2。 (3)算法流程 根據(jù)上述的遺傳算子的定義，并設(shè)置最大的迭代次數(shù)為1000，將GA算法流程表達(dá)如下， (i)隨機(jī)生成初始種群。 (ii)按照本節(jié)(1)中的公式計算各個個體適應(yīng)度的值。 (iii)判斷是否到達(dá)了最大的迭代次數(shù)。假設(shè)是，算法結(jié)束，輸出計算結(jié)果；假設(shè)否，轉(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.2GA算法計算結(jié)果 使用Matlab編程實(shí)現(xiàn)2.1中的算法，計算結(jié)果如下。 圖1GA算法過程的距離值變化 圖2GA算法求解的10城市TSP問題的結(jié)果 最后的結(jié)果編碼為〔89102314567〕，解為H-I-J-B-C-A-D-E-F-G，距離為269.0671。 從計算結(jié)果可以看出，算法迭代到300步時已經(jīng)收斂到了局部的最優(yōu)值。經(jīng)過大量的計算發(fā)現(xiàn)至多迭代到300步，算法收斂到局部最優(yōu)值。經(jīng)過進(jìn)一步的驗(yàn)證發(fā)現(xiàn)，這個局部最優(yōu)解也是全局最優(yōu)解。3使用PSO算法求解3.1算法描述 (1)TSP問題的交換子和交換序設(shè)n個節(jié)點(diǎn)的TSP問題的解序列為s=(ai)，i=1…n。定義交換子SO(i1,i2)為交換解S中的點(diǎn)ai1和ai2，那么S’=S+SO(i1,i2)為解S經(jīng)算子SO(i1,i2)操作后的新解。一個或多個交換子的有序隊(duì)列就是交換序，記作SS，SS=(SO1,SO2,…SON)。其中，SO1,…,SON等是交換子，之間的順序是有意義的,意味著所有的交換子依次作用于某解上。

假設(shè)干個交換序可以合并成一個新的交換序，定義為兩個交換序的合并算子。設(shè)兩個交換序SS1和SS2按先后順序作用于解S上，得到新解S’。假設(shè)另外有一個交換序SS’作用于同一解S上，能夠得到形同的解S’，可定義

和屬于同一等價集，在交換序等價集中，擁有最少交換子的交換序稱為該等價集的根本交換序。TSP問題的速度和位置更新算式根據(jù)(1)中的交換子和交換序的定義，可以根據(jù)根本的PSO算法速度更新算式和位置更新算式，重新定義如下的速度和位置更新算式， 式中，、為[0,1]區(qū)間的隨機(jī)數(shù)。為粒子與個體極值的交換序，以概率保存；為粒子與全局極值的交換序，以概率保存。粒子的位置按照交換序進(jìn)行更新。為慣性權(quán)重。 (3)算法流程 按照本節(jié)中的有關(guān)定義給出算法流程如下， (i)初始化粒子群，給每一個粒子一個初始解和隨機(jī)的交換序。 (ii)判斷是否到達(dá)最大迭代次數(shù)1000。假設(shè)是，算法結(jié)束，輸出結(jié)果；假設(shè)不是，轉(zhuǎn)到(iii)。 (iii)根據(jù)粒子當(dāng)前位置計算下一個新解： (a)計算A=，A是一個根本交換序，表示A作用于得到； (b)計算B=，B是一個根本交換序； (c)按照(2)中的公式更新速度和位置。 (d)如果得到了更好的個體位置，更新。 (iv)如果得到了更好的群體位置，更新。3.2PSO算法的計算結(jié)果 編程實(shí)現(xiàn)3.1中的算法，計算結(jié)果如下。計算程序見附錄。 圖3PSO算法過程的距離值變化 圖4PSO算法求解的10城市TSP問題的結(jié)果 最后的結(jié)果編碼為，解為A-D-E-F-G-H-I-J-B-C，距離為269.0671。 從計算結(jié)果可以看出，算法迭代到200步時已經(jīng)收斂到了局部的最優(yōu)值。這個局部最優(yōu)解也是全局最優(yōu)解。從收斂的速度的平均意義上而言，PSO算法與GA算法比收斂的更快。附錄%GA算法的代碼%GA算法程序.data=load('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);fori=1:nfarm(i,:)=randperm(N);endR=farm(1,:);len=zeros(n,1);fitness=zeros(n,1);counter=0;%GA算法迭代whilecounter<Cfori=1:nlen(i,1)=myLength(D,farm(i,:));endmaxlen=max(len);minlen=min(len);fitness=len;fori=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;%計算適應(yīng)度函數(shù)nn=0;fori=1:niffitness(i,1)>=randnn=nn+1;FARM(nn,:)=farm(i,:);endendFARM=FARM(1:nn,:);%FARM(nn*N)[aa,bb]=size(FARM);%(aa=nn)%交叉FARM2=FARM;fori=1:2:aaifPc>rand&&i<aaA=FARM(i,:);B=FARM(i+1,:);L=length(a);ifL<=10%確定交叉寬度W=9;elseif((L/10)-floor(L/10))>=rand&&L>10W=ceil(L/10)+8;elseW=floor(L/10)+8;endp=unidrnd(L-W+1);%隨機(jī)選擇交叉范圍fori=1:Wx=find(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));endFARM(i,:)=A;FARM(i+1,:)=B;endend%變異FARM2=FARM;fori=1:aaifPm>=randFARM(i,:)=mutate(FARM(i,:));endend%群體更新FARM2=zeros(n-aa+1,N);ifn-aa>=1fori=1:n-aaFARM2(i,:)=randperm(N);endendFARM=[R;FARM;FARM2];[aa,bb]=size(FARM);ifaa>nFARM=FARM(1:n,:);endfarm=FARM;clearFARMRR(counter+1)=myLength(D,R)%統(tǒng)計迭代次數(shù)。counter=counter+1;end%結(jié)果圖形顯示figureholdonplot([a(R(1),1),a(R(10),1)],[a(R(1),2),a(R(10),2)],'ms-','LineWidth',2,'MarkerEdgeColor','k','MarkerFaceColor','g');scatter(a(:,1),a(:,2),'ms')gridonholdonfori=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','MarkerFaceColor','g')holdonend%結(jié)果輸出RRlength%其他函數(shù)functiona=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;functionlen=myLength(D,p)[N,NN]=size(D);len=D(p(1,N),p(1,1));fori=1:(N-1)len=len+D(p(1,i),p(1,i+1));end%PSO算法代碼%%PSO算法代碼data=load('pos10.txt')cityCoor=[data(:,2)data(:,3)];n=size(cityCoor,1);cityDist=zeros(n,n);fori=1:nforj=1:nifi~=jcityDist(i,j)=((cityCoor(i,1)-cityCoor(j,1))^2+...(cityCoor(i,2)-cityCoor(j,2))^2)^0.5;endcityDist(j,i)=cityDist(i,j);endendindividual=zeros(indiNumber,n);%初始化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);forN=1:nMax%更新最優(yōu)值indiFit=fitness(individual,cityCoor,cityDist);fori=1:indiNumberifindiFit(i)<recordPbest(i)recordPbest(i)=indiFit(i);tourPbest(i,:)=individual(i,:);endifindiFit(i)<recordGbestrecordGbest=indiFit(i);tourGbest=individual(i,:);endend[value,index]=min(recordPbest);recordGbest(N)=recordPbest(index);%%更新每一個粒子的位置fori=1:indiNumber%個體最優(yōu)值更新速度c1=unidrnd(n-1);c2=unidrnd(n-1);whilec1==c2c1=round(rand*(n-2))+1;c2=round(rand*(n-2))+1;endchb1=min(c1,c2);chb2=max(c1,c2);cros=tourPbest(i,chb1:chb2);ncros=size(cros,2);forj=1:ncrosfork=1:nifxnew1(i,k)==cros(j)xnew1(i,k)=0;fort=1:n-ktemp=xnew1(i,k+t-1);xnew1(i,k+t-1)=xnew1(i,k+t);xnew1(i,k+t)=temp;endendendendxnew1(i,n-ncros+1:n)=cros;dist=0;forj=1:n-1dist=dist+cityDist(xnew1(i,j),xnew1(i,j+1));enddist=dist+cityDist(xnew1(i,1),xnew1(i,n));ifindiFit(i)>distindividual(i,:)=xnew1(i,:);end%全體最優(yōu)值更新速度c1=round(rand*(n-2))+1;c2=round(rand*(n-2))+1;whilec1==c2c1=round(rand*(n-2))+1;c2=round(rand*(n-2))+1;endchb1=min(c1,c2);chb2=max(c1,c2);cros=tourGbest(chb1:chb2);ncros=size(cros,2);forj=1:ncrosfork=1:nifxnew1(i,k)==cros(j)xnew1(i,k)=0;fort=1:n-ktemp=xnew1(i,k+t-1);xnew1(i,k+t-1)=xnew1(i,k+t);xnew1(i,k+t)=temp;endendendendxnew1(i,n-ncros+1:n)=cros;dist=0;forj=1:n-1dist=dist+cityDist(xnew1(i,j),xnew1(i,j+1));enddist=dist+cityDist(xnew1(i,1),xnew1(i,n));ifindiFit(i)>distindividual(i,:)=xnew1(i,:);end%%慣性項(xiàng)c1=round(rand*(n-1))+1;c2=round(ran