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

1、应用GA和PSO算法求解10城市TSP问题应用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）横坐标纵坐标城市A40 44.39城市B24.3914.63城市C17.0722.93城市D22.9376.1城市E51.7194.14城市F87.3265.36城市G68.7

2、852.19城市H84.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。根据种群中个体的距离并考虑使用自适应的标定方法，定义如下的适应度函数， 使用此适应度函数的后面的乘方次数可以调整来改变淘汰速度。这里

3、选择2，表示淘汰速度比较适中。 (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 G

5、A算法计算结果 使用Matlab编程实现2.1中的算法，计算结果如下。 图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。定义交换子

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

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

8、用于得到； (b) 计算B=，B是一个基本交换序； (c) 按照(2)中的公式更新速度和位置。 (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步时已经收敛到了局部的最优值。这个局部最优解也是全局最优解。从收敛的速度的平均意义上

9、而言，PSO算法与GA算法比收敛的更快。附录% GA算法的代码% GA算法程序. data = load(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 cou

10、nter=rand nn=nn+1; FARM(nn,:)=farm(i,:); end end FARM=FARM(1:nn,:); % FARM (nn*N) aa,bb=size(FARM);%(aa=nn) % 交叉 FARM2=FARM; for i=1:2:aa if Pcrand&iaa A=FARM(i,:); B=FARM(i+1,:); L=length(a);if L=rand&L10 W=ceil(L/10)+8;else W=floor(L/10)+8;endp=unidrnd(L-W+1);%随机选择交叉范围for i=1:W x=find(a=b(1,p+i-1)

11、; 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=1:n-aa FARM2(i,:)=randpe

12、rm(N); end end FARM=R;FARM;FARM2; aa,bb=size(FARM); if aan 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,MarkerEdgeColor,k,MarkerFaceColor,g);scatter(a(:,1)

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

14、),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); cityDist=zeros(n,n); for i=1:n for j=1:n if i=j cityDi

15、st(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);% 初始化nMax=1000; % 迭代次数indiNumber=50; % 粒子个数% 初始化个体和群体最优值indiFit=fitness(individual,cityCoor,cityDist);value,index=min(indiFit);tourPbest=individual; tourGb

16、est=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:indiNumber if indiFit(i)recordPbest(i) recordPbest(i)=indiFit(i); tourPbest(i,:)=indiv

17、idual(i,:); end if indiFit(i)dist individual(i,:)=xnew1(i,:); end % 全体最优值更新速度 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(

18、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

19、(i,:); end % 惯性项 c1=round(rand*(n-1)+1; c2=round(rand*(n-1)+1; while c1=c2 c1=round(rand*(n-2)+1; c2=round(rand*(n-2)+1; end temp=xnew1(i,c1); xnew1(i,c1)=xnew1(i,c2); xnew1(i,c2)=temp; 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 i

20、ndiFit(i)dist individual(i,:)=xnew1(i,:); end end value,index=min(indiFit); L_best(N)=indiFit(index); tourGbest=individual(index,:); end% 其他函数function dist=dist(x,D)n=size(x,2);dist=0;for i=1:n-1 dist=dist+D(x(i),x(i+1);enddist=dist+D(x(1),x(n);function indiFit=fitness(x,cityCoor,cityDist)m=size(x,1);n=size(cityCoor,1);indiFit=zeros(m,1);for i=1:m for j=1:n-1 indiFit(i)=indiFit(i)+cityDist(x(i,j),x(i,j+1); end indiFit(i)=indiFit(i)+cityDist(x(i,1),x(i,n);end