编程语言
首页 > 编程语言> > 实战1:基于遗传算法解决旅行商问题的MATLAB编程(3)TSP算法编写问题

实战1:基于遗传算法解决旅行商问题的MATLAB编程(3)TSP算法编写问题

作者:互联网

目录

   0.先上程序

   1.绪论+研究背景

   2.研究方法

   3.使用遗传算法编写TSP问题

         3.0 初始化定义+城市坐标分布编码和显示

         3.1 距离函数

         3.2 适应度函数

         3.3 选择算子

         3.4 交叉算子

         3.5 变异算子

         3.6 迭代

         3.7 作图

              3.7.1 每代最小值散点图

              3.7.2 总适应度折线图

              3.7.3 最优路径图

    4.参考


 

     本节内容是为了对之前程序篇内容的解释,如下章节

实战1:基于遗传算法解决旅行商问题的MATLAB编程(1)程序篇

实战1:基于遗传算法解决旅行商问题的MATLAB编程(2)绪论、研究背景以及研究方法

 

3.使用遗传算法编写TSP问题

 

     3.0 初始化定义+城市坐标分布编码和显示

(1)自己给定城市的坐标,这里举例用20个城市。

%% Define cities coordination
city20=[29 54;71 76;74 78;81 71;25 38;58 35;4 50;
13 40;18 40;24 98;71 44;64 60;69 58;83 69;58 69;54 63;51 67;37 84;
41 94;3 85];

(2)定义初始化变量,如种群个数,编码对象个数,迭代次数,实际迭代次数,变异概率,交叉概率等等

N=20;                            %城市数
city_coordinate=city20;          %城市位置,前面已经给点城市的数量,则可以确定城市的坐标。(8,15,20,30)
s=100;                           %样本数(种群)
c=30;                            %替换个数
pc=0.8;                          %交叉概率crossover
pm=0.1;                          %变异概率mutation
times=3000;                      %最大迭代次数
time=0;                          %实际迭代次数

(3)初始化变量,不同种群位置的矩阵定义、最短距离,最短路径编码排布

%%初始化变量,设定总群,总适应度以及最短距离以及最小距离。
pop=zeros(s,N+1);                %初始种群+适应度,之所以加N+1,因为最后还要回到原点。
pop_fit_aver=[];%总适应度
min_dis=[];                      %最短距离
pop_min=[];%最短距离的基因

(4)初始化种群,使用randperm为每一个种群产生随机的编码(随机的城市排布)

%初始化,其中s为种群数
for i=1:s   
    pop(i,1:N)=randperm(N);    %%randperm函数产生1-N内任意排列的N维数组。
end

(5)标出城市的位置用圆点,给城市位置编码并作图

plot(city_coordinate(:,1),city_coordinate(:,2),'bo');%画平均适应度折线图  //一个列向量作为自变量,另一个列向亮作为应变量。
%%变量处理:城市间
for i=1:N
    test_t=num2str(i);                            %%将变量类型转换为浮点型数
    text(city_coordinate(i,1),city_coordinate(i,2),test_t);%标号      行不变,列改变
end
grid on;
title('城市坐标分布图');
xlabel('x');
ylabel('y');

(6)效果图,可以清晰知道城市的位置,城市的编码序号排布。

                                           

   3.1 距离函数

function [city_distance] = CityDistance(city_coordinate,N)%城市距离矩阵
    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
end

代码解释:       

         由于我们知道城市坐标的位置,这样可以使用勾股定理求出任意两个城市之间的距离

city_distance之所以赋值为N BY N,是因为自己跟自己的距离在接下来城市距离计算中也存在。

通过双for循环,配合城市坐标,采用勾股定理求出距离,并填入相应两个城市之中。

         回到主程序(GeneticAlgorithm_TSP.m)调用该函数。

city_distance=CityDistance(city_coordinate,N);%城市间距离

3.2 适应度函数

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

代码解释:

individual_distance 是用来接收不同种群的路径的距离(编码城市不同会产生不同的旅行距离)。

                 

其中:    individual_fit=1./individual_distance;    (适应度函数变换,因为我们要求最小值,取反,求最大值即可)

不懂的朋友可以参考:https://www.cnblogs.com/yanshw/p/14749809.html 这篇文章,讲的很详细。

sum =所有列举种群中距离之和。

回到主程序,调用适应度库函数。

%适应度函数编写
[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,:)];

此时 min_dis  是最短种群编码的路径距离。

此时pop =  (S x(N+1)),具体内容如下:

                                               

POP_min 具体内容如下

                                             

 

3.3 选择算子(selection.m)

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

代码解释:

其中c在之间定义,替换的个数=30。

sortrows函数在help中的解释:(升序)

             

           

           

size 函数

size(A)=返回矩阵尺寸的行向量

size(A,1dim)=返回矩阵的行数

size(A,2dim)=返回矩阵的列数

size(A,3dim) =返回矩阵的第三维度dimension

返回主函数,下面解释选择算子在主函数中的作用。

%% 选择父代
pop=selection(pop,N,s,c);
for i=1:times
    time=time+1;
    E_new_new=[];   %子代
    for j=1:s/2    
        a=rand(1);
        b=rand(1);

代码解释:

X = rand(n) 返回一个 n×n 的随机数矩阵。

times:迭代次数

time:  实际迭代次数


3.4 交叉算子(crossover.m)

             交叉算子:交叉分为单点交叉和多点交叉。与传统的遗传算法不同的是,进行交叉不能整体替换,即替换后

结果访问的城市不能有重复或者依然是1-N的排列方式。

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

代码解释:

(1)交换本质:如同C语言中的交换变量:   

              temp =A

              A =B

              B =Ttemp

(2)tubelate 函数生成线性表格


SwapRepeat.m (交换函数,为了是将某一个种群中重复的编码去掉)

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

代码解释:

上图:简单明了。

               


主函数调用算子:

%%  交叉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

代码解释:

pc:交叉概率

X = rand(sz1,...,szN) 返回由随机数组成的 sz1×...×szN 数组,其中 sz1,...,szN 指示每个维度的大小。例如:rand(3,4) 返回一个 3×4 的矩阵。(生成一个由介于 0 和 1 之间的均匀分布的随机数组成的 5×5 矩阵)


3.5 变异算子(Mutation.m)

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

代码解释:

(1)crosspoint 接收随机变异点

(2)判断两个点是大小,从小到大还是从大到小

使用sub = 截取pop0中相应变异点的值

(3)使用SwapGene库函数交换变异点

(4)替换

 图片解释:

                         

其中算子 SwapGene.m

function [a]=SwapGene(sub,c1,c2)%交换
    kk=ceil((c2-c1)/2);
    kkk=(c2-c1)+2;
    for k=1:kk
        kkkk=sub(k);
        sub(k)=sub(kkk-k);
        sub(kkk-k)=kkkk;
    end
    a=sub;
end

代码解释:之所以kkk = c2-c1  本来数量角度需要加1,但是为了能够表示最后一位数之后的数,就需要加2

                                               

上述代码中交换的方式:

                                               

传统的穷举交换方式:

                                 

变异算子在主函数的调用

            
%% 变异mutation
            if b>pm         
                ;
            else
                pop(j,:)=Mutation(pop(j,:),N);
                pop(j+s/2,:)=Mutation(pop(j+s/2,:),N);
            end

代码解释:

pm:变异概率0.1


3.6 迭代所需求的内容:

(1)适应度函数sum

(2)目标函数集合(每代最小值的集合)min_dis

(3)最终迭代完后保留下来的种群的内容写入E_new_new中 

(4)每代种群最优编码方式集合 pop_min

(5) E_new_new=[];   每次迭代都美重新更新

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

作图:

(1)每代最小值散点图

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

 

效果图

                                       

(2)总适应度折线图

figure%画总适应度折线图
plot(time1,pop_fit_aver);
grid on;
title('总适应度折线图');
xlabel('迭代次数');
ylabel('每代总适应度');

 效果图

                                           

(3)最优路径图

figure%画最优路径
DrawPath(city_coordinate,pop_min,min_index,N)    %%参数:城市坐标,每次迭代最小种群的编码方式,最小在种群index,城市数量
grid on;
title('最优路径图');
xlabel('x');
ylabel('y');

效果图

                                       

这里讲一下DrawPath.m

function DrawPath(city_coordinate,E_new_new,min_index,N)%画路径图  %% %%参数:城市坐标,每次迭代最小种群的编码方式,最小在种群index,城市数量
    k=E_new_new(min_index,1:N)     %%去掉适应度的值
    %plot(kkk(:,1),kkk(:,2),'b');%画平均适应度折线图
    plot(city_coordinate(:,1),city_coordinate(:,2),'bo');
    hold on;
    for i=1:N-1
        plot([city_coordinate(k(i),1),city_coordinate(k(i+1),1)],[city_coordinate(k(i),2),city_coordinate(k(i+1),2)],'r','LineWidth',2);
        test_t=num2str(i);
        text(city_coordinate(k(i),1),city_coordinate(k(i),2),test_t);
        hold on;
    end
    test_t=[num2str(N)];
    text(city_coordinate(k(N),1),city_coordinate(k(N),2),test_t);
    plot([city_coordinate(k(N),1),city_coordinate(k(1),1)],[city_coordinate(k(N),2),city_coordinate(k(1),2)],'r','LineWidth',2)
end

注意点:最后一个城市返回出发点的城市。

    plot([city_coordinate(k(N),1),city_coordinate(k(1),1)],[city_coordinate(k(N),2),city_coordinate(k(1),2)],'r','LineWidth',2)

 

Refrence :

1.https://blog.csdn.net/qcyfred/article/details/76731706

2.https://www.cnblogs.com/yanshw/p/14749809.html

 

date:20220601/1:37:31 PM

         六一儿童节的礼物,加油!

标签:city,min,%%,pop,MATLAB,coordinate,遗传算法,TSP,distance
来源: https://www.cnblogs.com/sophiaecho/p/16314561.html