matlab遗传算法求解（TSP）附matlab代码

行左云合

已于 2023-10-16 18:21:37 修改

阅读量699

点赞数 4

文章标签： matlab 启发式算法

于 2023-10-16 18:20:00 首次发布

本文链接：https://blog.csdn.net/weixin_43855315/article/details/133863375

版权

问题描述

旅行商问题(traveling saleman problem, TSP)又被称为推销员问题、货郎担问题，该问题是最基本的路线问题。该问题寻求单一旅行者由起点出发，通过所有给定的需求点之后，最后再回到起点的最小路径成本，路径的限制是每个需求点只能拜访一次。最早的旅行商问题的数学模型是由Dantzig(1959)等学者提出的。

基本流程

算法原理不多赘述，主要接实例验证

（1）初始化：设置进化代数计数器t=0，设置最大进化代数T，随机生成M个个体作为初始群体P(0)。

（2）个体评价：计算群体中各个个体的适应度。

（3）选择：将选择算子作用于群体。选择的目的是把优化的个体直接遗传到下一代或通过配对交叉产生新的个体再遗传到下一代。

（4）交叉运算：将交叉算子作用于群体。遗传算法中起核心作用的就是交叉算子。

（5）变异：将变异算子作用于群体。即是对群体中的个体串的某些的基因值作变动。群体P(t)经过选择、交叉、变异之后得到下一代群体P(t+1)。

算法实现

以求遍历个点总距离最小为例。数据输入分为两种，一种是坐标输入（距离矩阵需要根据坐标计算），另一种是直接序号输入（各点之间距离已知，但坐标未知）。输入不同仅影响主函数部分，对变异、交叉等函数并不影响。参考了该文章。

主函数（坐标已知距离未知时）

%% Define cities coordination
city20(:,1) = 100 * rand(20,1); 
city20(:,2) = 100 * rand(20,1); 
N=20;                            %城市数
city_coordinate=city20;
s=100;                           %样本数（种群）
c=30;                            %替换个数
pc=0.8;                          %交叉概率crossover
pm=0.1;                          %变异概率mutation
times=3000;                      %最大迭代次数
time=0;                          %实际迭代次数
%%初始化变量，设定总群，总适应度以及最短距离以及最小距离。
pop=zeros(s,N+1);                %初始种群+适应度,之所以加N+1，因为最后还要回到原点。
pop_fit_aver=[];%总适应度
min_dis=[];                      %最短距离
pop_min=[];%最短距离的基因
%初始化，其中s为种群数
for i=1:s   
    pop(i,1:N)=randperm(N);    %%randperm函数产生1-N任意排列的N维无重复数组，代表城市访问顺序。
end
%%求距离矩阵
    city_distance=zeros(N,N); %%通过zeros构造
    for i=1:N
        for j=1:N
            city_distance(i,j)=((city_coordinate(i,1)-city_coordinate(j,1))^2+...
                (city_coordinate(i,2)-city_coordinate(j,2))^2)^0.5;     %%matlab是以矩阵的思维编写的，则每个距离都会存入city_diatance
        end
    end
%%

%适应度函数编写
[individual_fit,sum,min1,min_index]=GroupFit(city_distance,N,pop,s);   
         %%通过城市距离；城市数。pop（城市的随机分布）。s是城市分布的例子个数
sumP=sum;
pop_fit_aver=[pop_fit_aver;sum];        %%构造总群适应度函数
min_dis=[min_dis;min1];                 %%构造最小距离  min_dis=min1
pop(:,N+1)=individual_fit;              %%列矩阵的个体适应度函数
pop_min=[pop_min;pop(min_index,:)];
%% 选择父代
pop=selection(pop,N,s,c);
for k=1:times
    time=time+1;
    E_new_new=[];   %子代
    for j=1:s/2    
        a=rand(1);
        b=rand(1);
    %%  交叉crossover
            if a>pc          
                ;
            else
                crosspoint=rand(1,2);
                crosspoint=floor(crosspoint*N)+1;
                [pop(j,:),pop(j+s/2,:)]=Crossover(pop(j,:),pop(j+s/2,:),crosspoint,N);       
            end
                        
%% 变异mutation
            if b>pm         
                ;
            else
                pop(j,:)=Mutation(pop(j,:),N);
                pop(j+s/2,:)=Mutation(pop(j+s/2,:),N);
            end
            E_new_new=[E_new_new;pop(j,:);pop(j+s/2,:)]; %每次迭代保留优秀种群写E_new_new中
  
    end
    [individual_fit,sum,min1,min_index]=GroupFit(city_distance,N,E_new_new,s); 
                                                  %%适应度向量，适应度函数，最小距离，最小距离下标
                                                  %%参数：城市距离矩阵，城市数，每代最优种群的几何，种群数
    sumS=sum;
    pop_fit_aver=[pop_fit_aver;sum];              %%适应度函数
    min_dis=[min_dis;min1];                       %%每次迭代的最小值的向量集合
    E_new_new(:,N+1)=individual_fit;              %%目标函数（城市之间的距离）
    pop_min=[pop_min;E_new_new(min_index,:)];
    if(abs(sumS-sumP)<0.001)%退出条件             %%迭代的条件
        break;
    end
    pop=selection(E_new_new,N,s,c);               %%选择算子
end
[a,min_index]=min(min_dis);                       %%min函数可以求得 该向量min_dis最小值和最小值所在的index
time1=1:time+1;        %%为什么要加1，因为没有迭代前，原始种群也有一个最小值

time1=1:time+1;        %%为什么要加1，因为没有迭代前，原始种群也有一个最小值
figure%画平均适应度折线图
plot(time1,min_dis,'k.');
grid on;
title('每代最小值散点图');
xlabel('迭代的次数');
ylabel('最短的距离');

%%
最优序列储存在E_new_new（最后一行）中，最小距离储存于min_dis（最后一行）。
%%

主函数（各点之间距离已知，但无坐标）

最优序列储存在E_new_new（最后一行）中，最小距离储存于min_dis（最后一行）。

clear
N=[];                            %城市数
city_distance=zeros(N,N); %%通过zeros构造预分配内存
city_distance=[]; %直接输入各个点之间的距离矩阵（二维矩阵，比如横表头是1-9个点，竖表头也为这个顺序，中间内容为各个点之间的相互距离，该矩阵应有一条对角线均为0，并且数据关于该对角线对称）
s=100;                           %样本数（种群）
c=30;                            %替换个数
pc=0.8;                          %交叉概率crossover
pm=0.1;                          %变异概率mutation
times=3000;                      %最大迭代次数
time=0;                          %实际迭代次数
%%初始化变量，设定总群，总适应度以及最短距离以及最小距离。
pop=zeros(s,N+1);                %初始种群+适应度,之所以加N+1，因为最后还要回到原点。
pop_fit_aver=[];%总适应度
min_dis=[];                      %最短距离
pop_min=[];%最短距离的基因
%初始化，其中s为种群数
for i=1:s   
    pop(i,1:N)=randperm(N);   %%randperm函数产生1-N任意排列的N维无重复数组，代表城市访问顺序
end
 
%适应度函数编写
[individual_fit,sum,min1,min_index]=GroupFit(city_distance,N,pop,s);   
         %%通过城市距离；城市数。pop（城市的随机分布）。s是城市分布的例子个数
sumP=sum;
pop_fit_aver=[pop_fit_aver;sum];        %%构造总群适应度函数
min_dis=[min_dis;min1];                 %%构造最小距离  min_dis=min1
pop(:,N+1)=individual_fit;              %%列矩阵的个体适应度函数
pop_min=[pop_min;pop(min_index,:)];
%% 选择父代
pop=selection(pop,N,s,c);
for k=1:times
    time=time+1;
    E_new_new=[];   %子代
    for j=1:s/2    
        a=rand(1);
        b=rand(1);
    %%  交叉crossover
            if a>pc          
                ;
            else
                crosspoint=rand(1,2);
                crosspoint=floor(crosspoint*N)+1;
                [pop(j,:),pop(j+s/2,:)]=Crossover(pop(j,:),pop(j+s/2,:),crosspoint,N);       
            end
                        
%% 变异mutation
            if b>pm         
                ;
            else
                pop(j,:)=Mutation(pop(j,:),N);
                pop(j+s/2,:)=Mutation(pop(j+s/2,:),N);
            end
            E_new_new=[E_new_new;pop(j,:);pop(j+s/2,:)];           %%每次迭代保留优秀种群写入E_new_new中
  
    end
    [individual_fit,sum,min1,min_index]=GroupFit(city_distance,N,E_new_new,s); 
                                                  %%适应度向量，适应度函数，最小距离，最小距离下标
                                                  %%参数：城市距离矩阵，城市数，每代最优种群的几何，种群数
    sumS=sum;
    pop_fit_aver=[pop_fit_aver;sum];              %%适应度函数
    min_dis=[min_dis;min1];                       %%每次迭代的最小值的向量集合
    E_new_new(:,N+1)=individual_fit;              %%目标函数（城市之间的距离）
    pop_min=[pop_min;E_new_new(min_index,:)];
    if(abs(sumS-sumP)<0.001)%退出条件             %%迭代的条件
        break;
    end
    pop=selection(E_new_new,N,s,c);               %%选择算子
end
[a,min_index]=min(min_dis);                       %%min函数可以求得 该向量min_dis最小值和最小值所在的index
time1=1:time+1;        %%为什么要加1，因为没有迭代前，原始种群也有一个最小值
% 
time1=1:time+1;        %%为什么要加1，因为没有迭代前，原始种群也有一个最小值
figure%画平均适应度折线图
plot(time1,min_dis,'k.');
grid on;
title('每代最小值散点图');
xlabel('迭代的次数');
ylabel('最短的距离');

%%
最优序列储存在E_new_new（最后一行）中，最小距离储存于min_dis（最后一行）。

适应度函数

距离取倒数作为适应度，分母越小（距离越小，求最小距离），适应度越大。

function [individual_fit,num,min_distance,index] = GroupFit(city_distance,N,pop,s)%种群适应度
    individual_distance=zeros(s,1);     %%目的是求得最短路劲的值，给定INDIVIDUAL为sx1的矩阵
    for j=1:s
        sum_distance=0;
        for i=1:N-1
            sum_distance=sum_distance+city_distance(pop(j,i),pop(j,i+1));
        end
        sum_distance=sum_distance+city_distance(pop(j,N),pop(j,1));        %%最后回到起点
        individual_distance(j,1)=sum_distance;                             %%individual_distance是一个列矩阵
    end
    [min_distance,index]=min(individual_distance);    %%找到最短的距离
    individual_fit=1./individual_distance;           %%分子越小，分母越大
    num=0;
    for i=1:s
      num=num+individual_fit(i,1);
    end
end

选择算子

function [pop_ok]=selection(pop,N,s,c)%选择父代
    pop=sortrows(pop,N+1);
    for i=1:c
        pop(i,:)=pop(s+1-i,:);
    end
    randIndex=randperm(size(pop,1));
    pop=pop(randIndex,:);
    pop_ok=pop;
end

交叉算子

function [a,b]=Crossover(pop1,pop2,crosspoint,N)%交叉
    A=pop1;
    if(crosspoint(:,1)<crosspoint(:,2))    %%由于交叉的位置是随机的，有可能是从小到大，也有可能是从大到小
        pop1(crosspoint(:,1):crosspoint(:,2))=pop2(crosspoint(:,1):crosspoint(:,2));
        pop2(crosspoint(:,1):crosspoint(:,2))=A(1,crosspoint(:,1):crosspoint(:,2));
        while 1
            tbl = tabulate(pop1(1:N));   %%Create Frequency table
            if (tbl(:,3)<=(100/N))   %%如果个数超过一个就退出循环
                break;
            end
            [pop1,pop2]=SwapRepeat(tbl,pop1,pop2,crosspoint(:,1),crosspoint(:,2),N);%%执行存在重复的编码问题
        end
    else
        pop1(crosspoint(:,2):crosspoint(:,1))=pop2(crosspoint(:,2):crosspoint(:,1));
        pop2(crosspoint(:,2):crosspoint(:,1))=A(1,crosspoint(:,2):crosspoint(:,1));
        while 1
            tbl = tabulate(pop1(1:N));
            if (tbl(:,3)<=(100/N))
                break;
            end
            [pop1,pop2]=SwapRepeat(tbl,pop1,pop2,crosspoint(:,2),crosspoint(:,1),N);
        end
    end
    a=pop1;b=pop2;
end

去重函数（matlab自带去重函数 unique函数，所以个人认为该部分可以简化，如果懒得改那就直接用）

function [a,b]=SwapRepeat(tbl,pop1,pop2,c1,c2,N)%基因去重
    i=100/N;
    for k=1:(c1-1)              %%从左侧开始找起
        if tbl(pop1(k),3)>i     %%数重复出现
            kk=find(pop1(c1:c2)==pop1(k))+c1-1;  %%发现重复数在原本数组中的位置
            kkk=pop1(k);        %%备份自变量在k时的值
            pop1(k)=pop2(kk);   %%pop1和pop2
            pop2(kk)=kkk;       %%pop2
        end
    end
    for k=c2+1:N                %%从右侧开始找起
        if tbl(pop1(k),3)>i
            kk=find(pop1(c1:c2)==pop1(k))+c1-1;
            kkk=pop1(k);
            pop1(k)=pop2(kk);
            pop2(kk)=kkk;
        end
    end
    a=pop1;
    b=pop2;
end

变异算子

function [a]=Mutation(pop0,N)%变异
    crosspoint=rand(1,2);       %%随机产生一个1x2的，范围在0-1内的向量Vector
    crosspoint=floor(crosspoint*N)+1;
    if(crosspoint(:,1)<crosspoint(:,2))
        sub=pop0(crosspoint(:,1):crosspoint(:,2));
        sub=SwapGene(sub,crosspoint(:,1),crosspoint(:,2));
        pop0(crosspoint(:,1):crosspoint(:,2))=sub;
    else
        sub=pop0(crosspoint(:,2):crosspoint(:,1));
        sub=SwapGene(sub,crosspoint(:,2),crosspoint(:,1));
        pop0(crosspoint(:,2):crosspoint(:,1))=sub;
    end
    a=pop0;
end

其中算子 SwapGene.m

function [a]=SwapGene(sub,c1,c2)%交换
    kk=ceil((c2-c1)/2);
    kkk=(c2-c1)+2; %本来数量角度需要加1，但是为了能够表示最后一位数之后的数，就需要加2，反正代码可以跑。
    for k=1:kk
        kkkk=sub(k);
        sub(k)=sub(kkk-k);
        sub(kkk-k)=kkkk;
    end
    a=sub;
end