应用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城市H84.8836.09

2、城市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。根据种群中个体的距离并考虑使用自适应的标定方法，定义如下的适应度函数，_Xi距离-种群最小距离2f (Xi) =(1 _种群最大距离-种群最小距离)使用此适应度函数的

3、后面的乘方次数可以调整来改变淘汰速度。这里选择2 ,表示淘汰速度比较适中。(2) 定义算子选择算子，根据适应度函数的大小进行选择。计算每个个体被选中的概率P(Xi )= Nf(X),以各个个体所分配到的概率值作为其遗传到下一代的概率，基' f(Xi) i日丁这些概率用赌盘选择法来产生下一代群体。交义算子，采用部分映射交义(Partially Mapped Crossover, PMX)方法实现算 法交义。首先选取选需要交义的区间段， 然后确定区问段的映射关系，接下来交 换区问段的遗传信息，此时未换部分的遗传信息会出现不合法的情况， 因此根据 将未换部分按映射关系进行交换。交义率为 0.

4、9。变异算子，把一个染色体中的两个基因的交换作为变异算法。在算法中随机 的产生一个染色体中需要交换的两个基因的位置，将这两个基因进行交换来实现 变异。变异率为0.2。(3) 算法流程根据上述的遗传算子的定义，并设置最大的迭代次数为1000,将GA算法流 程叙述如下，(i) 随机生成初始种群。(ii) 按照本节(1 )中的公式计算各个个体适应度的值。(iii) 判断是否达到了最大的迭代次数。若是，算法结束，输出计算结果；若 否，转到(iv)。(iv) 按照本节(2)中的选择算子进行选择操作。(v) 选择交义宽度并随机确定交义切点，按照本节(2)中的交义算子进行交义操作。(vi) 按照本节(2)中

5、的变异算子进行变异操作。(vii) 将遗传算子产生的种群作为新的种群，转到(ii)。2.2 GA算法计算结果使用Matlab编程实现2.1中的算法，计算结果如下GA算法过程380360340值离320距3002801,1260 01002003004005006007008009001000图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步，算法收敛到局部最优值

6、。经过进一步的验 证发现，这个局部最优解也是全局最优解。3使用PSO算法求解3.1算法描述(1) TSP问题的交换子和交换序设n个节点的TSP问题的解序列为s=(ai), i=1n。定义交换子SO(i1,i2)为 交换解S中的点ai1和ai2,则S' =S+SO(i1,i2)为解S经算子SO(i1,i2)操作后的 新解。一个或多个交换子的有序队列就是交换序，记作SS, SS=(SO1,SO2,SON)。其中，SO1,SON等是交换子，之间的顺序是有意义的，意味着所有的交换子 依次作用丁某解上。若干个交换序可以合并成一个新的交换序，定义为两个交换序的合并算子。设两个交换序SS1和SS2按

7、先后顺序作用丁解S上，得到新解S'。假设另 外有一个交换序SS'作用丁同一解S上，能够得到形同的解S',可定义SS = SS3S§SS'和SS由SS届丁同一等价集，在交换序等价集中，拥有最少交换子的交换序称为该等价集的基本交换序。(2) TSP问题的速度和位置更新算式根据(1)中的交换子和交换序的定义，可以根据基本的PSO算法速度更新算式和位置更新算式，重新定义如下的速度和位置更新算式，U =仍 (Pid Xd) E(pgd Xd) :X =Xid +Vid式中，a、P为0,1区间的随机数。Pid-Xid为粒子与个体极值的交换序，以概率 a保留；Pgd

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

9、(iV)如果得到了更好的群体位置，更新 Pgd。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 = load( 'pos10.txt&

10、#39; );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:nfarm(i,:)=randperm(N);endR=farm(1,:);len=zeros(n,1);fitness=zeros(n,1);counter=0;% GA算法迭代while counter<Cfor i=1:nlen(i,1)=myLength(D,farm(i,:);endmaxlen=max(len);minlen=min(l

11、en);fitness=len;for i=1:length(len)f itness(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:nif fitness(i,1)>=randnn=nn+1;FARM(nn,:)=farm(i,:);endendFARM=FARM(1:nn,:);% FARM (nn*N)aa,bb=size(FARM);%(aa=nn)%交叉FARM2=FARM;

12、for i=1:2:aaif Pc>rand&&i<aaA=FARM(i,:);B=FARM(i+1,:);L=length(a);if L<=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);咖机选择交叉范围for i=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

13、,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;for i=1:aaif Pm>=randFARM(i,:)=mutate(FARM(i,:);endend%群体更新FARM2=zeros(n-aa+1,N);if n-aa>=1for i=1:n-aaFARM2(i,:)=randperm(N);endendFARM=R;FARM;FARM2;aa,bb=size(FARM);if aa>nFARM=FARM(1:

14、n,:);endfarm=FARM;clearFARMRR(counter+1)=myLength(D,R)%统计迭代次数。counter=counter+1 ;end%结果图形显ZKfigurehold onplot(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),'

15、;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' , 'MarkerFaceColor' , 'g')hold onend%结果输出RRlength%其他函数function a=mutate(a)L=length(a);rray=ran

16、dperm(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);end% PSO算法代码% PSO算法代码data=load( 'pos10.txt' )cityCoor=data(:,2) data(:,3);n=size(cityCoor,1

17、); cityDist=zeros(n,n);for i=1:nfor 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;endcityDist(j,i)=cityDist(i,j);endendindividual=zeros(indiNumber,n);%初始化nMax=1000;%迭代次数indiNumber=50;% 粒子个数%初始化个体和群体最优值indiFit=fitness(individual,cityCoor,cityDist); val

18、ue,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=1:nMax%更新最优值indiFit=fitness(individual,cityCoor,cityDist);for i=1:indiNumberif indiFit(i)<recordPbest(i) rec

19、ordPbest(i)=indiFit(i); tourPbest(i,:)=individual(i,:);endif indiFit(i)<recordGbest recordGbest=indiFit(i); tourGbest=individual(i,:);endendvalue,index=min(recordPbest);recordGbest(N)=recordPbest(index);%更新每一个粒子的位置for i=1:indiNumber%个体最优值更新速度c1=unidrnd(n-1);c2=unidrnd(n-1);while c1=c2 c1=round(ra

20、nd*(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);for j=1:ncrosfor 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;endendendendxnew1(i,n-ncros+1:n)=cros;dist=0;

21、for j=1:n-1dist=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-2)+1;c2=round(rand*(n-2)+1;while c1=c2c1=round(rand*(n-2)+1;c2=round(rand*(n-2)+1;end chb1=min(c1,c2);chb2=max(c1,c2);cro

22、s=tourGbest(chb1:chb2);ncros=size(cros,2);for j=1:ncrosfor k=1:nif xnew1(i,k)=cros(j) xnew1(i,k)=0;for t=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;for j=1:n-1dist=dist+cityDist(xnew1(i,j),xnew1(i,j+1); enddist=dist+cityDist(xnew1(i,1),xnew1(i,n);if indiFit(i)>distindividual(i,:)=xne