遗传算法解决城市TSP问题

一、遗传算法与TSP问题简介

1.遗传算法是模拟达尔文生物进化论的自然选择和遗传学机理的生物进化过程的计算模型，是一种通过模拟自然进化过程搜索最优解的方法。遗传算法是从代表问题可能潜在的解集的一个种群开始的，而一个种群则由经过基因编码的一定数目的个体组成。每个个体实际上是染色体带有特征的实体。染色体作为遗传物质的主要载体，即多个基因的集合，其内部表现（即基因型）是某种基因组合，它决定了个体的形状的外部表现。因此，在一开始需要实现从表现型到基因型的映射即编码工作。初代种群产生之后，按照适者生存和优胜劣汰的原理，逐代演化产生出越来越好的近似解，在每一代，根据问题域中个体的适应度大小选择体，并借助于自然遗传学的遗传算子进行组合交叉和变异，产生出代表新的解集的种群。这个过程将导致种群像自然进化一样的后生代种群比前代更加适应于环境，末代种群中的最优个体经过解码，可以作为问题近似最优解。
其基本流程如图所示：

2.旅行商问题，即TSP问题（Traveling Salesman Problem）又译为旅行推销员问题、货郎担问题，是数学领域中著名问题之一。假设有一个旅行商人要拜访n个城市，他必须选择所要走的路径，路径的限制是每个城市只能拜访一次，而且最后要回到原来出发的城市。路径的选择目标是要求得的路径路程为所有路径之中的最小值。

二、遗传算法的基本原理

遗传算法是一类借鉴生物界的进化规律(适者生存，优胜劣汰遗传机制)演化而来的随机化搜索方法，其主要特点是直接对结构对象进行操作，不存在求导和函数连续性的限定，具有内在的隐并行性和更好的全局寻优能力，采用概率化的寻优方法，能自动获取和指导优化的搜索空间，自适应地调整搜索方向，不需要确定的规则。

三、遗传算法的主要步骤

1.初始化种群：在可行域内随机产生若干个解，称其为初始种群。
2.对种群内的每个解进行适应度评估，评价每个个体对环境的适应程度强弱。
3.适者生存，不适者淘汰，运用随机数等方式，将不适合的个体淘汰，适合的保留。
4.通过交叉互换、变异等方式对中云进行扰动，实质上就是在可行域内进行搜索，搜索出最适合环境的解。
上述的环境是指目标函数，需要求解最优值的函数，每个个体就是每个可行解，适应度评估就是将可行解带入函数进行计算然后对比其大小，如果可行解代入函数的结果最优，那这个解就是最优解。优胜劣汰是对优秀解的保留，较差解的淘汰，交叉和变异就是对解进行扰动，以进一步提高解的优度。

四、遗传算法的特点

1.遗传算法只对个体的基因进行操作，所以无论实际问题多么复杂，其稳定性都不会受到太大的影响。
2.遗传算法的搜索过程属于并行计算，能够很好地搜索解空间。
3.稳定性、鲁棒性强，适用于非线性、高维复杂优化问题。

五、算法实现代码模块

1.连点画图函数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;
endend

2.染色体的路程代价函数 mylength.m

function len=myLength(D,p)%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.适应度函数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

4.交叉操作函数 cross.m

function [A,B]=cross(A,B)
L=length(A);
if L<10W=L;
elseif ((L/10)-floor(L/10))>=rand&&L>10W=ceil(L/10)+8;
elseW=floor(L/10)+8;
end
%%W为需要交叉的位数
p=unidrnd(L-W+1);%随机产生一个交叉位置
%fprintf('p=%d ',p);%交叉位置
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,p+i-1),B(1,p+i-1));[A(1,x),B(1,y)]=exchange(A(1,x),B(1,y));
endend

5.对调函数 exchange.m

function [x,y]=exchange(x,y)
temp=x;
x=y;
y=temp;end

6.变异函数 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

7.主函数

%main
clear;
clc;
%%%%%%%%%%%%%%%输入参数%%%%%%%%
N=25;               %%城市的个数
M=100;               %%种群的个数
ITER=2000;               %%迭代次数
%C_old=C;
m=2;                %%适应值归一化淘汰加速指数
Pc=0.8;             %%交叉概率
Pmutation=0.05;       %%变异概率
%%生成城市的坐标pos=randn(N,2);
%randn是一种产生标准正态分布的随机数或矩阵的函数，属于MATLAB函数
%返回一个N*2的随机项的矩阵。如果N不是个数量，将返回错误信息。%%生成城市之间距离矩阵
D=zeros(N,N);
%zeros(N,N)产生N×N的double类零矩阵
for i=1:Nfor j=i+1:Ndis=(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);
%zeros(M,N)产生M×N的double类零矩阵
for i=1:Mpopm(i,:)=randperm(N);%随机排列，比如[2 4 5 6 1 3]
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);%找到最小值的下标，赋值为rr
R=popm(rr(1,1),:);%提取该染色体，赋值为R
for i=1:Nfprintf('%d ',R(i));%把R顺序打印出来
end
fprintf('\n');fitness=fitness/sum(fitness);
distance_min=zeros(ITER+1,1);  %%各次迭代的最小的种群的路径总长
nn=M;
iter=0;
while iter<=ITERfprintf('迭代第%d次\n',iter);%%选择操作p=fitness./sum(fitness);q=cumsum(p);%累加for i=1:(M-1)len_1(i,1)=myLength(D,popm(i,:));r=rand;tmp=find(r<=q);popm_sel(i,:)=popm(tmp(1),:);end [fmax,indmax]=max(fitness);%求当代最佳个体popm_sel(M,:)=popm(indmax,:);%%交叉操作nnper=randperm(M);
%    A=popm_sel(nnper(1),:);%   B=popm_sel(nnper(2),:);%%for i=1:M*Pc*0.5A=popm_sel(nnper(i),:);B=popm_sel(nnper(i+1),:);[A,B]=cross(A,B);%      popm_sel(nnper(1),:)=A;%      popm_sel(nnper(2),:)=B; popm_sel(nnper(i),:)=A;popm_sel(nnper(i+1),:)=B;end%%变异操作for i=1:Mpick=rand;while pick==0pick=rand;endif pick<=Pmutationpopm_sel(i,:)=Mutation(popm_sel(i,:));endend%%求适应度函数NN=size(popm_sel,1);len=zeros(NN,1);for i=1:NNlen(i,1)=myLength(D,popm_sel(i,:));endmaxlen=max(len);minlen=min(len);distance_min(iter+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:Nfprintf('%d ',R(i));endfprintf('\n');popm=[];popm=popm_sel;iter=iter+1;%pause(1);end
%end of whilefigure(3)
plot_route(pos,R);
axis([-3 3 -3 3]);
figure(4)
plot(distance_min);

六、控制变量分析各参数对解的影响

1.种群规模对解的影响
不变的参数如下：
N=25; //城市的个数
ITER=2000; /迭代次数
m=2; //适应值归一化淘汰加速指数
Pc=0.8; //交叉概率
Pmutation=0.05; //变异概率
①M=20; //种群的个数

②M=50; //种群的个数

③M=100; //种群的个数

结论：种群规模越大算法结果越精确，适应度越好，但是运行时间越久。
2.迭代次数对解的影响
不变的参数如下：
N=25; //城市的个数
M=100; //种群的个数
m=2; //适应值归一化淘汰加速指数
Pc=0.8; //交叉概率
Pmutation=0.05; //变异概率
①ITER=1000; /迭代次数

②ITER=1500; /迭代次数

③ITER=2000; /迭代次数

迭代次数越大解的正确率越高。当迭代次数达到1000时算法运行结果趋于稳定(稳定得到最短路径)，这时算法的效率比较高。算法的迭代次数最优值为2000。
3.城市数量对解的影响
不变的参数如下：
M=100; //种群的个数
ITER=2000; /迭代次数
m=2; //适应值归一化淘汰加速指数
Pc=0.8; //交叉概率
Pmutation=0.05; //变异概率
①N=15; //城市的个数

②N=20; //城市的个数

③N=25; //城市的个数

结论：当城市数量较多时，大于50个城市，迭代多次，GA仍然不收敛，可能的问题是陷入了局部最优解
4.交叉概率对解的影响
不变的参数如下：
N=25; //城市的个数
M=100; //种群的个数
ITER=2000; /迭代次数
m=2; //适应值归一化淘汰加速指数
Pmutation=0.05; //变异概率
①Pc=0.6; //交叉概率

②Pc=0.7; //交叉概率

③Pc=0.8; //交叉概率

结论：交叉概率过低将得不到最优解，交叉概率越高平均适应度越好。
5.变异概率对解的影响
不变的参数如下：
N=25; //城市的个数
M=100; //种群的个数
ITER=2000; /迭代次数
m=2; //适应值归一化淘汰加速指数
Pc=0.8; //交叉概率
①Pmutation=0.01; //变异概率

②Pmutation=0.05; //变异概率

③Pmutation=0.1; //变异概率

结论：变异概率过高或者过低都影响得到最优解。

七、总结

由于遗传算法的整体搜索策略和优化计算是不依赖梯度信息，只需要影响搜索方向的目标函数和相应的适应度函数，所以它的应用比较广泛。利用遗传算来进行大规模问题的组合优化是一种比较有效的方法，因为在目前的计算上利用枚举法求解最优解有一定的难度。但是遗传算法也有不足之处，它对算法的精确度、可行度、计算复杂性等方面还没有有效的定量分析方法。通过本文的算法也可以清晰地认识到，遗传算法所求得的解不一定是最优解。