应用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）横坐标纵坐标城市A4044.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)编码和适应度函数分别用1-10表示城市A-J，然后采用自然数编码方式为TSP问题编码，例如，旅程（16289105734）表示从城市A出发分别经过了F-B-H-I-J-E-G-C-D的一次旅行。每一个问题的解及算法中的个体都可以计算相应的距离。那么种群中的最小距离和最大距离也相应的可以确定。选择种群个数为50。根据种群中个体的距离并考虑使用自适应的标定方法，定义如下的适应度函数，2i)-x1()(种群最小距离种群最大距离种群最小距离距离ixf使用此适应度函数的后面的乘方次数可以调整来改变淘汰速度。这里选择2，表示淘汰速度比较适中。(2)定义算子选择算子，根据适应度函数的大小进行选择。计算每个个体被选中的概率Niiiixfxfxp1)()()(，以各个个体所分配到的概率值作为其遗传到下一代的概率，基于这些概率用赌盘选择法来产生下一代群体。交叉算子，采用部分映射交叉(PartiallyMappedCrossover,PMX)方法实现算法交叉。首先选取选需要交叉的区间段，然后确定区间段的映射关系，接下来交换区间段的遗传信息，此时未换部分的遗传信息会出现不合法的情况，因此根据将未换部分按映射关系进行交换。交叉率为0.9。变异算子，把一个染色体中的两个基因的交换作为变异算法。在算法中随机的产生一个染色体中需要交换的两个基因的位置，将这两个基因进行交换来实现变异。变异率为0.2。(3)算法流程根据上述的遗传算子的定义，并设置最大的迭代次数为1000，将GA算法流程叙述如下，(i)随机生成初始种群。(ii)按照本节(1)中的公式计算各个个体适应度的值。(iii)判断是否达到了最大的迭代次数。若是，算法结束，输出计算结果；若否，转到(iv)。(iv)按照本节(2)中的选择算子进行选择操作。(v)选择交叉宽度并随机确定交叉切点，按照本节(2)中的交叉算子进行交叉操作。(vi)按照本节(2)中的变异算子进行变异操作。(vii)将遗传算子产生的种群作为新的种群，转到(ii)。2.2GA算法计算结果使用Matlab编程实现2.1中的算法，计算结果如下。01002003004005006007008009001000260280300320340360380GA算法过程迭代次数距离值图1GA算法过程的距离值变化102030405060708090102030405060708090100kmkmGA算法求解10个城市TSP问题规划路径规划路径图2GA算法求解的10城市TSP问题的结果最后的结果编码为（89102314567），解为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=1…n。定义交换子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’，可定义21'SSSSSS'SS和21SSSS属于同一等价集，在交换序等价集中，拥有最少交换子的交换序称为该等价集的基本交换序。(2)TSP问题的速度和位置更新算式根据(1)中的交换子和交换序的定义，可以根据基本的PSO算法速度更新算式和位置更新算式，重新定义如下的速度和位置更新算式，ididididgdididididvxxxpxpvv)()(式中，、为[0,1]区间的随机数。ididxp为粒子与个体极值的交换序，以概率保留；idgdxp为粒子与全局极值的交换序，以概率保留。粒子的位置按照交换序idv进行更新。为惯性权重。(3)算法流程按照本节中的有关定义给出算法流程如下，(i)初始化粒子群，给每一个粒子一个初始解idx和随机的交换序idv。(ii)判断是否达到最大迭代次数1000。若是，算法结束，输出结果；若不是，转到(iii)。(iii)根据粒子当前位置计算下一个新解：(a)计算A=ididxp，A是一个基本交换序，表示A作用于idx得到idp；(b)计算B=idgdxp，B是一个基本交换序；(c)按照(2)中的公式更新速度和位置。(d)如果得到了更好的个体位置，更新idp。(iv)如果得到了更好的群体位置，更新gdp。3.2PSO算法的计算结果编程实现3.1中的算法，计算结果如下。计算程序见附录。01002003004005006007008009001000260280300320340360380PSO算法过程迭代次数距离值图3PSO算法过程的距离值变化102030405060708090102030405060708090100kmkm规划路径规划路径图4PSO算法求解的10城市TSP问题的结果最后的结果编码为（14567891023），解为A-D-E-F-G-H-I-J-B-C，距离为269.0671。从计算结果可以看出，算法迭代到200步时已经收敛到了局部的最优值。这个局部最优解也是全局最优解。从收敛的速度的平均意义上而言，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);fori=1:nfarm(i,:)=randperm(N);endR=farm(1,:);len=zeros(n,1);fitness=zeros(n,1);counter=0;%GA算法迭代whilecounterCfori=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;%计算适应度函数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:aaifPcrand&&iaaA=FARM(i,:);B=FARM(i+1,:);L=length(a);ifL=10%确定交叉宽度W=9;elseif((L/10)-floor(L/10))=rand&&L10W=ceil(L/10)+8;elseW=floor(L/10)+8;endp=unidrnd(L-W+1);%随机选择交叉范围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);ifaanFARM=FARM(1:n,:);endfarm=FARM;clearFARMRR(counter+1)=myLength(D,R)%统计迭代次数。counter=counter+1;end%结果图形显示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%结果输出RRlength%其他函数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;%迭代次数indiNumber=50;%粒子个数%%初始化个体和群体最优值indiFit=fitness(individual,cityCoor,c