遗传算法求解TSP问题

一、概念

TSP问题
TSP问题：已给一个n个点的完全图，每条边都有一个长度，求总长度最短的经过每个顶点正好一次的封闭回路。

路径的选择目标是要求得的路径路程为所有路径之中的最小值。

遗传算法
遗传算法（Genetic Algorithm）是模拟达尔文生物进化论的自然选择和遗传学机理的生物进化过程的计算模型，是一种通过模拟自然进化过程搜索最优解的方法。
遗传算法的这些性质，已被人们广泛地应用于组合优化、机器学习、信号处理、自适应控制和人工生命等领域。
其做法是把问题参数编码为染色体，再利用迭代的方式进行选择、交叉、以及变异等运算来交换种群中的染色体信息，最终生成符合优化目标的染色体。

遗传算法的主要特点是直接对结构对象进行操作，不存在求导和函数连续性的限定，具有更好的全局寻优能力，采用概率化的寻优方法，能自动获取和指导优化的搜索空间，自适应的调整搜索方向，不需要确定的规则。

二、算法

遗传算法的基本流程

1. 初始化 设置最大迭代次数ITER和种群大小M,个体染色体长度N(这里为城市的个数)，随机产生初始群体。
2. 适应度函数 计算群体中各个个体的适应度。
3. 选择运算 选择操作即从旧种群中以一定概率选择个体到新种群中，个体被选中的概率跟适应度值有关，个体适应度值越大，被选中的概率就越大。选择的目的是把优化的个体之间遗传到下一代或通过配对交叉产生新的个体再遗传到下一代。
4. 交叉运算 将交叉算子作用于群体，采用部分映射杂交，将父代样本部分结构加以替换重组而产生新的个体。
5. 变异运算 将变异算子作用于群体 变异操作采取随机选取两个点，将其对换位置即可。
6. 终止判断 迭代次数是否到达设定的最大值，到达最大值即可结束。

流程图

主要程序
程序一：主程序
计算路径的遗传算法main.m

%main
clear;
clc;
%%%%%%%%%%%%%%%输入参数%%%%%%%%
N=25;               %%城市的个数
M=100;               %%种群的个数
ITER=2000;               %%迭代次数
%C_old=C;
m=2;                %%适应值归一化淘汰加速指数
Pc=0.8;             %%交叉概率
Pmutation=0.05;       %%变异概率
%%生成城市的坐标
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);%随机排列，比如[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:N
    fprintf('%d ',R(i));%把R顺序打印出来
end
fprintf('\n');


fitness=fitness/sum(fitness);
distance_min=zeros(ITER+1,1);  %%各次迭代的最小的种群的路径总长
nn=M;
iter=0;
while iter<=ITER
    fprintf('迭代第%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.5
        A=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:M
        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(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:N
        fprintf('%d ',R(i));
    end
    fprintf('\n');
    popm=[];
    popm=popm_sel;
    iter=iter+1;
    %pause(1);
    end
%end of while
    
figure(3)
plot_route(pos,R);
axis([-3 3 -3 3]);
figure(4)
plot(distance_min);

程序二：交叉函数

%交叉操作函数  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
%%W为需要交叉的位数
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

程序三：变异函数

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

程序四：变异函数

%对调函数 exchange.m
function [x,y]=exchange(x,y)
temp=x;
x=y;
y=temp;
end

程序五：计算路径长度的函数

%染色体的路程代价函数  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

程序六：计算归一化适应值计算函数

%适应度函数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

程序七：绘制路径图函数

%连点画图函数 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

三、运行结果

1.在种群数量100，交叉概率0.8，迭代次数为2000时

N=25;               %%城市的个数
M=100;               %%种群的个数
ITER=2000;               %%迭代次数
m=2;                %%适应值归一化淘汰加速指数
Pc=0.8;             %%交叉概率
Pmutation=0.05;       %%变异概率

结果：在接近三百次迭代时趋于平缓。
2.在种群数量500，交叉概率0.8，迭代次数为2000时

结果：在五百次迭代时趋于平缓。
3.在种群数量1000，交叉概率0.8，迭代次数为2000时

结果：在四百次左右迭代时趋于平缓。
4.在种群数量100，交叉概率0.9，迭代次数为2000时


结果：在两百次四十次左右迭代时趋于平缓。
5.在种群数量100，交叉概率0.7，迭代次数为2000时


结果：在四百五十次迭代时趋于平缓。
6.在种群数量100，交叉概率0.8，加速指数1时

结果：在一百次左右迭代时趋于平缓。
7.在种群数量100，交叉概率0.8，加速指针3时


结果：在六百五十次左右迭代时趋于平缓。

四、总结

图4.1：基于七次实验得到的结果

M为种群数量
N为城市个数
ITER迭代次数2000
m为加速指数，,最好取为1,2,3,4 ,不宜太大
1.前三次实验时，路径最短的是第三次种群数量增大，取得的路径值越短，且迭代次数并没有大幅度增加，说明种群数量的增大有利于最短路径的生成。
2.在1，4，5次实验中改变交叉概率的大小，当交叉概率为0.9时，最短路径最短。
3.在1，6，7次实验中改变加速指数的大小，当加速指数为时1，最短路径最短。
4.遗传算法只能得到问题的近似最优解，而且对不同的问题，该算法的性能也不一样。