Matlab遗传算法与TSP问题的结合

Matlab遗传算法与TSP问题的结合

准备

在编写自己的代码之前应该了解的是算法的基本结构，具体可参考我之前的文章大纲或者这里简单函数优化亦或者是更为进阶的非线性规划问题
文章虽然简单，但是也能了解基本该设置的参数与具体的算法结构。

遗传算法+TSP问题

TSP问题是经典的NP问题，遗传算法在改问题上的解决并不完善。

实际上目前为止还没找到一个多项式时间的有效算法

但是对于简单的TSP问题，采用遗传算法是完全足够的。在本章中，我会详细的阐述轮盘选择，单点交叉算法，基本位变异运算。
为解决下列最短路问题${}^{[2]}$，即从A到G的最短路径问题。

TSP问题

在两者结合之前，应该先了解什么叫做TSP问题。展开TSP为：$TravelingSalesmanProblem$，即旅行商问题。假设有一个旅行商人要拜访$n$个城市，他必须选择所要走的路径，路径的限制是每个城市只能拜访一次，而且最后要回到原来出发的城市。路径的选择目标是要求得的路径路程为所有路径之中的最小值。（摘自百度百科）
或者说有一个搜索自然子集 X = { 1 , 2 , 3 , ⋅ ⋅ ⋅ , n } X=\left\{ 1,2,3,\cdot \cdot \cdot ,n \right\} X={1,2,3,⋅⋅⋅,n}
存在一个排列Y：
Y ( X ) = { V 1 , V 2 , ⋅ ⋅ ⋅ , V n } Y\left( X \right) =\left\{ V_1,V_2,\cdot \cdot \cdot ,V_n \right\} Y(X)={V1​,V2​,⋅⋅⋅,Vn​}
使得目标函数，走的路线和最小
$$Min=\sum _{i=1}^{n-1}d\left(V_i,V_{i+1}\right)+d\left(V_n,V_1\right)$$

邻接矩阵

邻接矩阵（Adjacency Matrix）是表示顶点之间相邻关系的矩阵。如若已知各点的"坐标"且为无向图的话可以使用以下函数获取邻接矩阵。(这里不仅指笛卡尔坐标系下的坐标，更可以指一个顶点包含的信息)

function [d]=distance_Adjacency(data)			%data为各顶点的坐标信息
    [m,n]=size(data);
    d=zeros(m,m);
    for i=1:m
        for j=i+1:m
            dm=abs(data(i,:)-data(j,:));		%计算距离，欧氏距离。
            d(i,j)=sqrt(sum(dm.^2,2));
        end
    end
    d=d+d';										%无向图，应为对称矩阵
end

但是很多情况遇见的不是无向图，有时更没有坐标值，只有两个顶点间的距离。
这里提供了笔者自己编写的函数$Get\_Adjacency$

• $Get\_Adjacency$

function Adjacency_matrix = Get_Adjacency(mode,s)
%%输入mode=1时需要输入xlsx的地址
%%输入mode=0时输入节点距离生成邻接矩阵
    if nargin==1
        if  mode==0
            disp('输入参数');
            [Adjacency_matrix,~]=TSP_myAdjacency();
        else
            err('参数错误');
        end
    elseif nargin==2
        if mode==1
            Adjacency_matrix=xlsread(s);
        else
            err('参数错误');    
        end
    end
end

输入mode=1时需要输入xlsx的地址
输入mode=0时输入节点距离生成邻接矩阵
在进行遗产算法之前先运行$Get\_Adjacency\left(0\right)$生成邻接矩阵。

• 无向图
  运行$test\_wx.m$
  
  上图仅为操作示例，保存在示例1.xlsx
• 有向图
  运行$test\_yx.m$
  
  需要设定断路值，在这里因为示例中的距离在10以内，因此设置了100，输入小于100值认为该路为通路。具体距离设置需要根据实际情况定夺。
  上图仅为操作示例，保存在示例2.xlsx

遗传算法

选择了谢菲尔德大学的Matlab遗传算法工具箱实现算法。
工具箱的导入极为简单，可以自行在网上下载导入。
由于本文中关于交叉与变异的操作方法与适应度无关，仅选择操作与适应度有关，因此选择操作决定了遗传算法整体的选取值趋向（最大值还是最小值）

种群初始化和计算初始种群适应度

• 初始参数设置和读取邻接矩阵

%%清屏
clear
clc
close all
tic
%% 初始化
Adjacency_matrix=Get_Adjacency(1,'power.xlsx');	%邻接矩阵
NIND=80;                                        %种群规模
MAXGEN=200;                                     %迭代次数
px=0.7;                                         %交叉概率
pm=0.01;                                        %变异概率
gen=1;
Fitmean=[];
Y=[];
[PM,PN]=size(Adjacency_matrix);

• 生成二进制初始种群代表路线经过此地
• 计算初始适应度

%生成初始种群
Chrom=TSP_population(NIND,PN);						%初始化种群
Chrom=Limit_population(Chrom);						%对种群进行限制（起点与终点被选中）
[N,L]=size(Chrom);
%显示参数
fprintf('Number of generations:%d\n',MAXGEN);
fprintf('Population size:%d\n',N);
fprintf('Crossover probability:%.3f\n',px);
fprintf('Mutation probability:%.3f\n',pm);
%计算适应度，画图
FitnV=Fit_fun_TSP(Chrom,Adjacency_matrix);
figure(1);
grid on;
hold on;
plot(FitnV,'k*');

• $TSP\_population$

function Chrom=TSP_population(NIND,PRECI)%%创建二进制种群
    Chrom=round(rand(NIND,PRECI));			%round取整
end

• $Limit\_population$

function Chrom=Limit_population(Chrom)
%%对种群个体进行限制
    Chrom(:,1)=1;						%选中起点
    Chrom(:,end)=1;						%选中终点
end

在本章中解决的仅仅为简单的TSP问题，没有过多的约束，具体情况需要具体分析。

计算适应度

• 创建独立的.m文件计算适应度

function F=Fit_fun_TSP(Chrom,Adjacency_matrix)%种群个体适应值，v为种群，Adjacency_matrix为系数矩阵(邻接矩阵)，
    [N,~]=size(Chrom);
    F=zeros(N,1);
    for i=1:N
        fp=find(Chrom(i,:)==1);
        [~,n]=size(fp);
        for j=1:n-1
            F(i)=F(i)+Adjacency_matrix(fp(j),fp(j+1));	%根据邻接矩阵求距离
        end
    end 
end

个体选择

• 采用轮盘方法从种群中选择个体

Fitmax=max(FitnV);
FitnV=Fitmax-FitnV;									%对适应度进行处理
SelCh=Roulette_selection(Chrom,FitnV);              %轮盘选择

因为轮盘选择趋向为最大值，因此对适应度进行处理，使得选择操作整体为最小值

• $Roulette\_selection$

function SelCh= Roulette_selection(Chrom,FitnV)%%轮盘选择
    [N,L]=size(Chrom);
    S=sum(FitnV);							%S为种群适应度的累积值
    for i=1:N
        SI=S*rand(1);						%数值随机，保障多样性
        for j=1:N
            if SI<=sum(FitnV(1:j))			%当前的个体的适应度越大，其和大于SI的可能性就越大
                SelCh(i,:)=Chrom(j,:);		%选取当前个体
                break;
            end
        end
    end
end

轮盘赌局选择法，适应度越大，在轮盘上占的面积越大，被选中的概率也就越大，但也有可能选取到适应度小的个体，一定程度保证了种群的多样性。

交叉

• 采用单点交叉算子对个体的染色体进行重组

SelCh=Single_point_crossover(SelCh,px);         	%单点交叉运算

• $Single\_point\_crossover$

function SelCh = Single_point_crossover(Chrom,pc)
    [N,L]=size(Chrom);
    C(:,1)=rand(N,1)<=pc;			%取随机数，选取小于交叉概率的个体					
    fc=find(C(:,1)==1);
    for i=1:2:size(fc)				%间隔一个方便选取两个染色体
        if i>=size(fc)
            break;
        end
        site=fix(1+L*rand(1));		%随机标定单点位置
        %选取个体，当前个体与下一个体
        temp=Chrom(fc(i,1),:);		
        Chrom(fc(i,1),site:end)=Chrom(fc(i+1,1),site:end);
        Chrom(fc(i+1,1),site:end)=temp(:,site:end);
    end
    SelCh=Chrom;
end

变异

• 对个体进行变异

SelCh=Bit_mutation(SelCh,pm);                       %基本位变异运算

• $Bit\_mutation$

function SelCh = Bit_mutation(Chrom,pm)
    [N,L]=size(Chrom);
    M=rand(N,L)<=pm;				%取随机数，选取整个种群小于变异概率的基因位（与考虑个体情况）
    SelCh=Chrom-2.*(Chrom.*M)+M;	%若当前为1，则变为0，若当前为0，则变为1
end

计算当前种群适应度记录最优解

• 计算当前所有个体的适应度。
• 找到并记录当代最优个体，更新最优解与平均解。

FitnV=Fit_fun_TSP(Chrom,Adjacency_matrix);
    %保留中间运算结果    
    [sol,indb]=min(FitnV);
    Chrom(1,:)=Chrom(indb,:);
    media=mean(FitnV);
    Y=[Y sol];                                          %每代最优
    Fitmean=[Fitmean media];                            %每代平均路径值
    gen=gen+1;

当平均解与最优解接近重合时说明当前种群中大多接近最优解，在一定程度代表当前种群已经较为稳定。

算法计算结果

路径，距离与个体结果如下：

种群起始位置与最终位置对比如下：

进化过程如下：

可见最优值稳定，平均值贴近最优值，说明算法的效果尚可，因为算法具有一定的随机性，因此需要运行多次取最佳值。

代码地址

完整代码在 GitHub网站上，如果有条件可以支持一下。现在关注我你就是老粉了

参考文章

本文内容参考下列文章，如有任何问题还请联系该账号

[1]史峰. MATLAB智能算法30个案例分析[M]. 北京航空航天大学出版社, 2011.
[2]汪晓银，李治，周保平主编.数学建模与数学实验 第3版[M].北京：科学出版社.2019.