应用GA和PSO算法求解10城市TSP问题

1、精选优质文档-倾情为你奉上应用GA和PSO算法求解10城市TSP问题1 问题描述旅行团计划近期在城市A、B、C、D、E、F、G、H、I和J共10个城市间进行一次周游旅行，为了尽量节省旅行开支，希望能找到一条里程数最少或相对较少的旅行路线。问题1，给定10个城市之间的公路里程如表1所示，并要求使用GA算法求解优化问题。问题2，与问题1数据相同，要求使用PSO算法求解优化问题。表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) 编码和适应度函数分别用1-10表示城市A-J，然后采用自然数编码方式为TSP问题编码，例如，旅程（1 6 2 8 9 10 5 7 3 4）表示从城市A出发分别经过了F-B-H-I-J-E-G-C-D的一次旅行。每一个问题的解及算法中的个体都可以计算相应的距离。那么种群中的最小距离和最大距离也相应的可以确定。选择种群个数为50。根据种群中个体的距离并考虑使用自适应的标定方法，定义如下的适应度函数，使用此适应度函数的后面的乘方次数可以调整来改变淘汰速度。这里选择2，表示淘汰速度比较

3、适中。(2) 定义算子选择算子，根据适应度函数的大小进行选择。计算每个个体被选中的概率，以各个个体所分配到的概率值作为其遗传到下一代的概率，基于这些概率用赌盘选择法来产生下一代群体。交叉算子，采用部分映射交叉(Partially Mapped Crossover, PMX)方法实现算法交叉。首先选取选需要交叉的区间段，然后确定区间段的映射关系，接下来交换区间段的遗传信息，此时未换部分的遗传信息会出现不合法的情况，因此根据将未换部分按映射关系进行交换。交叉率为0.9。变异算子，把一个染色体中的两个基因的交换作为变异算法。在算法中随机的产生一个染色体中需要交换的两个基因的位置，将这两个基因进行交换

4、来实现变异。变异率为0.2。(3) 算法流程根据上述的遗传算子的定义，并设置最大的迭代次数为1000，将GA算法流程叙述如下，(i) 随机生成初始种群。(ii) 按照本节(1)中的公式计算各个个体适应度的值。(iii) 判断是否达到了最大的迭代次数。若是，算法结束，输出计算结果；若否，转到(iv)。(iv) 按照本节(2)中的选择算子进行选择操作。(v) 选择交叉宽度并随机确定交叉切点，按照本节(2)中的交叉算子进行交叉操作。(vi) 按照本节(2)中的变异算子进行变异操作。(vii) 将遗传算子产生的种群作为新的种群，转到(ii)。2.2 GA算法计算结果使用Matlab编程实现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步时已经收敛到了局部的最优值。经过大量的计算发现至多迭代到300步，算法收敛到局部最优值。经过进一步的验证发现，这个局部最优解也是全局最优解。3 使用PSO算法求解3.1 算法描述(1)TSP问题的交换子和交换序设n个节点的TSP问题的解序列为s=(ai)，i=1n。定义交换子SO(i1,i2)为交换解S中的点ai1和ai2，则S=S+SO

6、(i1,i2)为解S经算子SO(i1,i2)操作后的新解。 一个或多个交换子的有序队列就是交换序，记作SS，SS=(SO1,SO2,SON)。其中，SO1,SON等是交换子，之间的顺序是有意义的, 意味着所有的交换子依次作用于某解上。若干个交换序可以合并成一个新的交换序，定义为两个交换序的合并算子。设两个交换序SS1和SS2按先后顺序作用于解S上，得到新解S。假设另外有一个交换序SS作用于同一解S上，能够得到形同的解S，可定义和属于同一等价集，在交换序等价集中，拥有最少交换子的交换序称为该等价集的基本交换序。(2) TSP问题的速度和位置更新算式根据(1)中的交换子和交换序的定义，可以根据基本

7、的PSO算法速度更新算式和位置更新算式，重新定义如下的速度和位置更新算式，式中，、为0,1区间的随机数。为粒子与个体极值的交换序，以概率保留；为粒子与全局极值的交换序，以概率保留。粒子的位置按照交换序进行更新。为惯性权重。(3) 算法流程按照本节中的有关定义给出算法流程如下，(i) 初始化粒子群，给每一个粒子一个初始解和随机的交换序。(ii) 判断是否达到最大迭代次数1000。若是，算法结束，输出结果；若不是，转到(iii)。(iii) 根据粒子当前位置计算下一个新解：(a) 计算A=，A是一个基本交换序，表示A作用于得到；(b) 计算B=，B是一个基本交换序；(c) 按照(2)中的公式更新速

8、度和位置。(d) 如果得到了更好的个体位置，更新。(iv) 如果得到了更好的群体位置，更新。3.2 PSO算法的计算结果编程实现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步时已经收敛到了局部的最优值。这个局部最优解也是全局最优解。从收敛的速度的平均意义上而言，PSO算法与GA算法比收敛的更快。附录% GA算法的代码% GA算法程序. data = l

9、oad('pos10.txt');a=data(:,2) data(:,3); n=50; % 种群数目C=1000; % 迭代次数 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; %计算适应度函数 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)% 统计迭代次数。 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% 其他函数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; % 迭代次数indiNumber=50; % 粒子个数% 初始化个体和群体最优值indiFit=fitness(individual,cityCoor,cityDist);value,index=min(indiFit);tourPbest=individual; tourGbest=individual(index,:) ; recordPbest=inf*ones(1,indiNumber); recordGbest=indiFit(index); xnew1=individual;% 循环寻找最优路径L_best=zeros(1,nMax);for N

18、=1:nMax % 更新最优值 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 % 个体最优值更新速度 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、 % 全体最优值更新速度 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=