遗传算法与TSP问题的MATLAB实现

最近越来越发现很难不受外界干扰的进行学习，可能与九月份这个躁动的求职季节有关，看着师兄们每天忙着奔走于各个公司的宣讲会，心中有种莫名的心情，时常想起大学毕业时的情景：四月经历考研失败；五月忙于毕业设计；六月刚毕业答辨完就和同学离开学校奔走于武汉各招聘会；两个星期后终于将自己“卖”出去；七月做着人生的一份工作。。。时常会清晰的感觉时间的紧迫，但还是得按步就部，不能冒进，一步一个脚印，坚持学习与提高自身素质并重。

    上一个星期被导师安排做神经网络相关，在做优化神经网络时涉及到遗传算法，于是搜集资料，参照别人部分程序，初步完成遗传算法解决TSP问题。

   “旅行商问题”(Traveling Salesman Problem,TSP)可简单描述为：一位销售商从n个城市中的某一城市出发，不重复地走完其余n-1个城市并回到原出发点，在所有可能路径中求出路径长度最短的一条。

旅行商的路线可以看作是对n城市所设计的一个环形，或者是对一列n个城市的排列。由于对n个城市所有可能的遍历数目可达(n-1)!个，因此解决这个问题需要O(n!)的计算时间。而由美国密执根大学的Holland教授发展起来的遗传算法，是一种求解问题的高效并行全局搜索方法，能够解决复杂的全局优化问题，解决TSP问题也成为遗传算法界的一个目标。

     遗传算法具有广泛的应用领域，它借助于生物进化的思想和原理与计算机科学相结合，在解决实际问题中得到了很好的广泛应用。遗传算法一般由选择、交叉、变异构成。它通过不断地迭代，逐步改进当前解，直到最后搜索到最优解或满意解。算法流程图如下：

                     

遗传算法求解TSP的基本步骤

（1）       种群初始化。个体编码方法有二进制编码和实数编码，在解决TSP问题过程中个体编码方法为实数编码。对于TSP问题，实数编码为1-n的实数的随机排列，初始化的参数有种群个数M、染色体基因个数N（即城市的个数）、迭代次数C、交叉概率Pc、变异概率Pmutation。

（2）       适应度函数。在TSP问题中，对于任意两个城市之间的距离D(i,j)已知，每个染色体（即n个城市的随机排列）可计算出总距离，因此可将一个随机全排列的总距离的倒数作为适应度函数，即距离越短，适应度函数越好，满足TSP要求。

（3）       选择操作。遗传算法选择操作有轮盘赌法、锦标赛法等多种方法，本程序采用基于适应度比例的选择策略，即适应度越好的个体被选择的概率越大，同时在选择中保存适应度最高的个体。

（4）       交叉操作。遗传算法中交叉操作有多种方法。本程序中对于个体，随机选择两个个体，在对应位置交换若干个基因片段，同时保证每个个体依然是1-n的随机排列，防止进入局部收敛。

（5）       变异操作。本程序中对于变异操作，随机选取个体，同时随机选取个体的两个基因进行交换以实现变异操作。

 

以下为选择N不同时对应的遗传算法所得到的最短距离连接图：

N=8时

在N=9以下时，计算机通过计算全排列得到的最优解和遗传算法得到的结果一致，而在N大于9时，计算机需要很长时间（等了几个小时都没计算出最优解），而遗传算法则可得到次优解或满意解。

MATLAB实现程序如下：

（1）       适应度函数fit.m

function fitness=fit(len,m,maxlen,minlen)

fitness=len;

for i=1:length(len)

    fitness(i,1)=(1-(len(i,1)-minlen)/(maxlen-minlen+0.0001)).^m;

end

(2)个体距离计算函数  mylength.m

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

 

 

end

(3)交叉操作函数  cross.m

function [A,B]=cross(A,B)

L=length(A);

if L<10

    W=L;

elseif ((L/10)-floor(L/10))>=rand&&L>10

    W=ceil(L/10)+8;

else

    W=floor(L/10)+8;

end

p=unidrnd(L-W+1);

fprintf('p=%d ',p);

for i=1:W

    x=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));

end

 

end

(4)对调函数 exchange.m

function [x,y]=exchange(x,y)

temp=x;

x=y;

y=temp;

 

end

(5)变异函数 Mutation.m

function a=Mutation(A)

index1=0;index2=0;

nnper=randperm(size(A,2));

index1=nnper(1);

index2=nnper(2);

%fprintf('index1=%d ',index1);

%fprintf('index2=%d ',index2);

 

temp=0;

temp=A(index1);

A(index1)=A(index2);

A(index2)=temp;

a=A;

end

(6)连点画图函数 plot_route.m

function plot_route(a,R)

scatter(a(:,1),a(:,2),'rx');

hold on;

plot([a(R(1),1),a(R(length(R)),1)],[a(R(1),2),a(R(length(R)),2)]);

hold on;

for 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);

    hold on;

end

 

end

(7)主函数

clear;

clc;

%%%%%%%%%%%%%%%输入参数%%%%%%%%

N=50;               %%城市的个数

M=100;               %%种群的个数

C=100;               %%迭代次数

C_old=C;

m=2;                %%适应值归一化淘汰加速指数

Pc=0.4;             %%交叉概率

Pmutation=0.2;       %%变异概率

%%生成城市的坐标

pos=randn(N,2);

%%生成城市之间距离矩阵

D=zeros(N,N);

for i=1:N

    for j=i+1:N

        dis=(pos(i,1)-pos(j,1)).^2+(pos(i,2)-pos(j,2)).^2;

        D(i,j)=dis^(0.5);

        D(j,i)=D(i,j);

    end

end

%%生成初始群体

popm=zeros(M,N);

for i=1:M

    popm(i,:)=randperm(N);

end

%%随机选择一个种群

R=popm(1,:);

 

figure(1);

scatter(pos(:,1),pos(:,2),'rx');

axis([-3 3 -3 3]);

figure(2);

plot_route(pos,R);      %%画出种群各城市之间的连线

axis([-3 3 -3 3]);

%%初始化种群及其适应函数

fitness=zeros(M,1);

len=zeros(M,1);

for i=1:M

    len(i,1)=myLength(D,popm(i,:));

end

maxlen=max(len);

minlen=min(len);

fitness=fit(len,m,maxlen,minlen);

rr=find(len==minlen);

R=popm(rr(1,1),:);

for i=1:N

fprintf('%d ',R(i));

end

fprintf('\n');

fitness=fitness/sum(fitness);

 

distance_min=zeros(C+1,1);  %%各次迭代的最小的种群的距离

while C>=0

fprintf('迭代第%d次\n',C);

%%选择操作

nn=0;

for i=1:size(popm,1)

    len_1(i,1)=myLength(D,popm(i,:));

    jc=rand*0.3;

    for j=1:size(popm,1)

        if fitness(j,1)>=jc

        nn=nn+1;

        popm_sel(nn,:)=popm(j,:);

        break;

        end

    end

end

%%每次选择都保存最优的种群

popm_sel=popm_sel(1:nn,:);

[len_m len_index]=min(len_1);

popm_sel=[popm_sel;popm(len_index,:)];

 

%%交叉操作

nnper=randperm(nn);

A=popm_sel(nnper(1),:);

B=popm_sel(nnper(2),:);

for i=1:nn*Pc

[A,B]=cross(A,B);

popm_sel(nnper(1),:)=A;

popm_sel(nnper(2),:)=B;

end

%%变异操作

for i=1:nn

    pick=rand;

    while pick==0

         pick=rand;

    end

    if pick<=Pmutation

       popm_sel(i,:)=Mutation(popm_sel(i,:));

    end

end

%%求适应度函数

NN=size(popm_sel,1);

len=zeros(NN,1);

for i=1:NN

    len(i,1)=myLength(D,popm_sel(i,:));

end

maxlen=max(len);

minlen=min(len);

distance_min(C+1,1)=minlen;

fitness=fit(len,m,maxlen,minlen);

rr=find(len==minlen);

fprintf('minlen=%d\n',minlen);

R=popm_sel(rr(1,1),:);

for i=1:N

fprintf('%d ',R(i));

end

fprintf('\n');

popm=[];

popm=popm_sel;

C=C-1;

%pause(1);

end

figure(3)

plot_route(pos,R);

axis([-3 3 -3 3]);