matlab遗传算法解tsp旅行商问题实战

最新推荐文章于 2024-07-21 10:17:10 发布

nanjc

最新推荐文章于 2024-07-21 10:17:10 发布

阅读量1.3k

点赞数 6

文章标签： matlab 开发语言

本文链接：https://blog.csdn.net/qq_49617724/article/details/129728785

版权

1.引言

个人在学习遗传算法解tsp旅行商问题时发现代码中有部分疏漏之处；有些大佬没标注的地方对我来说也有些难以理解，故重新整理一下并成功复现，也加了一些说明。

本文只相当于一个小笔记，想要仔细理解的请参考原文。我主要参考了以下两篇文章：链接1 ，链接2。

在对比链接1与链接2的模拟效果后，发现链接1采用轮盘赌法,进化逆转的结果明显较好，所以此处只复刻了链接1中算法（当然运算速度也更慢）。

2.代码复现

首先是原始数据，采用链接2中随机生成的原始城市数据。新建脚本运行就行。

city_location = zeros(30,2);
city_location(:,1) = 100 * rand(30,1);
city_location(:,2) = 100 * rand(30,1);
%更改初始城市数时把30全改掉就行
save city_location

结果如下：

然后直接上链接1的主程序（有点大）（其中修改了会一直划线且图像混杂的问题）：

clear;
clc;
close all;

load("city_location.mat"); % 加载数据
 
location = load("city_location.mat"); % 开始对数据进行可视化
x = location.city_location(:,1);
y = location.city_location(:,2);

plot(x,y,'o');
xlabel('城市横坐标');
ylabel('城市纵坐标');
grid on;

NIND = 100;         %种群大小
MAXGEN = 500;       %最大迭代次数；由于初始城市数较少，我这里500次迭代已经完全够用了
Pc = 0.9;           %交叉概率，相当于基因遗传的时候染色体交叉
Pm = 0.05;          %染色体变异
GGAP = 0.9;         %这个是代沟，通过遗传方式得到的子代数为父代数*GGAP
D = Distance(city_location);    %通过这个函数可以计算i,j两点之间的距离
N = size(D,1);      %计算有多少个坐标点
%%初始化种群
Chrom = InitPop(NIND,N);    %Chrome代表的种群
%%画出随机解得路线图
DrawPath(Chrom(1,:),city_location)
pause(0.0001)
%输出随机解的路线和总距离
disp('初始种群中的一个随机值')
OutputPath(Chrom(1,:));%其中一个个体
Rlength = PathLength(D,Chrom(1,:));
disp(['总距离;',num2str(Rlength)]);
disp('~~~~~~~~~~~~~~~~~~~~~~~~~~~~~')
%优化
gen = 0;

figure();
hold on;
box on;
xlim([0,MAXGEN])
trace = zeros(1,MAXGEN);
title('优化过程')
xlabel('迭代次数')
ylabel('当前最优解')

ObjV = PathLength(D,Chrom);%计算当前路线长度，即上面随机产生的那些个体路径
preObjV = min(ObjV);%当前最优解
while gen<MAXGEN
    %%计算适应度
    ObjV = PathLength(D,Chrom);     %计算路线长度
    pause(0.0001);
    preObjV = min(ObjV);
    trace(1,gen+1)=preObjV;
    FitnV = Fitness(ObjV);
    %选择
    SelCh = Select(Chrom,FitnV,GGAP);
    %交叉操作
    SelCH = Recombin(SelCh,Pc);
    %变异
    SelCh = Mutate(SelCh,D,Pm);
    %逆转操作
    SelCh = Reverse(SelCh,D);
    %重插入子代的新种群
    Chrom = Reins(Chrom,SelCh,ObjV);
    %更新迭代次数
    gen = gen + 1;
end
%画出最优解的路线图
plot([1:MAXGEN],trace(1,:));

ObjV = PathLength(D,Chrom);     %计算路线长度
[minObjV,minInd] = min(ObjV);
DrawPath(Chrom(minInd(1),:),city_location)


%%输出最优解的路线和距离
disp('最优解:')
p = OutputPath(Chrom(minInd(1),:));
disp(['旅行商走过的总距离：',num2str(ObjV(minInd(1)))]);
disp('~~~~~~~~~~~~~~~~~~~~~~~~~~~~~')

一共会出来四个图，分别是初始城市位置图，初始路线图，最终优化路线图，以及最后迭代效果图。还会有outputpath提供的最终点位顺序文字以及最终总路径长度。

函数：每一次新建函数，运行就行。

1.计算两两城市的距离并返回距离矩阵。类似

（0 1 3

1 0 4

3 4 0）

function D = Distance(a)
%计算两两城市之间的距离
%输入的数据为所有城市的x、y坐标位置矩阵city_location，输出为两两城市的距离构成的矩阵D
 
row = size(a,1);
 
D = zeros(row,row);
for i = 1:row
    for j = i+1:row
        D(i,j) = ((a(i,1)-a(j,1))^2+(a(i,2)-a(j,2))^2)^0.5;
        D(j,i) = D(i,j);
    end
end

2.计算初始种群。

function Chrom = InitPop(NIND,N)
%输入：
%NIND：种群大小
%N：个体染色体长度（城市个数）
%输出：初始种群
Chrom = zeros(NIND,N); % 定义存储种群的变量
for i = 1:NIND
    Chrom(i,:) = randperm(N);
end

3.在城市之间画出旅行商行走路线图

function DrawPath(Chrom,X)
%输入：
%待画路线
%城市的坐标位置
%输出：
%旅行商的路线
 
R =  [Chrom(1,:) Chrom(1,1)]; %第一个随机解（个体）【Chrom(1,:)取第一行数据】，一共有n个城市，但是这里R有n+1个值，因为后面又补了一个Chrom(1,1)，“是为了让路径最后再回到起点”
figure;
hold on
%plot(X(:,1),X(:,2),'bo')%X(:,1)，X(:,2)分别代表的X轴坐标和Y轴坐标
plot(X(:,1),X(:,2),'o','color',[0.5,0.5,0.5])%X(:,1)，X(:,2)分别代表的X轴坐标和Y轴坐标，
plot(X(Chrom(1,1),1),X(Chrom(1,1),2),'rv','MarkerSize',20)%标记起点

for i=1:size(X,1)
    text(X(i,1)+0.05,X(i,2)+0.05,num2str(i),'color',[1,0,0]);
end    %标注城市序号

A = X(R,:);                         %A是第一个随机解R的坐标矩阵
row = size(A,1);                    %row为坐标数+1
for i = 2:row
    [arrowx,arrowy] = dsxy2figxy(gca,A(i-1:i,1),A(i-1:i,2));    %dsxy2figxy坐标转换函数，annotation画图范围坐标值必须在（0，1）区间内，对a的坐标进行转换
    annotation('textarrow',arrowx,arrowy,'HeadWidth',8,'color',[0,0,1]);%将这两个点连接起来
                                                                        %annotation中颜色由[0,0,1]数组形式设置，'textarrow'对应文本式箭头
end
hold off
xlabel('横坐标x')
ylabel('纵坐标y')
title('旅行商轨迹图')
box on

4.以文字形式输出城市间行走的路径

function p=OutputPath(R)
%打印路线函数
%以1->2->3的形式在命令行打印路线
%以文字形式输出结果路径

R = [R,R(1)];
N = length(R);
p = num2str(R(1));
for i = 2:N
    p = [p,'->',num2str(R(i))];
end
disp(p)

5.计算每一个个体（行走路线）的路线总长度

function len = PathLength(D,Chrom)
%计算所有个体的路线长度
%输入
%D两两城市之间的距离
%Chrom个体的轨迹

[~,col] = size(D); %返回D的列数 
NIND = size(Chrom,1);%NIND等于初始种群
len = zeros(NIND,1);%初始化一个大小等于NIND的len来记录长度
for i = 1:NIND
    p = [Chrom(i,:) Chrom(i,1)];%构造p矩阵保存路线图 将第一行路线提出 再加上第一个构成回路
    i1 = p(1:end-1);%i1从第一个开始遍历到倒数第二个
    i2 = p(2:end);%i2从第二个开始遍历到倒数第一个
    len(i,1) = sum(D((i1-1)*col+i2));%计算出每种路线（种群的个体）的长度，这里相当于把D矩阵沿行展开
   %https://www.ilovematlab.com/thread-498868-1-1.html
   %这里len的求值，利用了在矩阵中，各个元素都有唯一标号的特性。
   %举例：例如下面这个矩阵D代表距离矩阵，D(i,j)，则代表i和j城市的距离
   % D =
   %     0         0.5052    0.6137    0.2908
   %     0.5052    0         0.4236    0.4305
   %     0.6137    0.4236    0         0.5503
   %     0.2908    0.4305    0.5503    0
   %那么对应数值的索引就应该是
   %   D =
   %     1         5         9         13
   %     2         6         10        14
   %     3         7         11        15
   %     4         8         12        16
   %D是4x4矩阵说明，只有4个城市，那么假设一个Chrom解是[3 2 1 4]
   %则p = [3 2 1 4 3],i1 = [3 2 1 4],i2 = [2 1 4 3]，i1-1 = [2 1 0 4]
   %现在计算城市间的距离，首先计算城市3->城市2的距离，对应D(2,3)
   %而D(2,3)对应的索引是10，而10=2*4+2
   %其中等号右边的第一个2，代表i1-1的第一个元素，4代表列数，第二个2代表i2的第一个元素
   %以此类推可以找出城市间的距离，最后sum就OK
end

6.轮盘赌法选择个体

function NewChrIx = Sus(FitnV,Nsel)
%输入：
%FitnV 个体的是适应度值
%Nsel 被选个体的数目
%输出：
%NewChrIx 被选择个体的索引号
[Nind,ans] = size(FitnV);%Nind为FitnV的行数也就是个体数 ans为列数1
cumfit = cumsum(FitnV);%对适应度累加 例如 1 2 3 累加后 1 3 6 用来计算累积概率
trials = cumfit(Nind)/Nsel * (rand + (0:Nsel-1)');%cumfit(Nind)表示的是矩阵cumfit的第Nind个元素 A.'是一般转置，A'是共轭转置 rand返回一个在区间 (0,1) 内均匀分布的随机数
%cumfit(Nind)/Nsel平均适应度 * （rand +(0:Nsel-1)'）从0开始到89的转置矩阵（行矩阵变列矩阵）加上每一项加上一个0-1的随机数
%生成随机数矩阵 用来筛选
Mf = cumfit(:,ones(1,Nsel));%将生成的累积概率 复制90份 生成一个100*90的矩阵
Mt = trials(:,ones(1,Nind))';
[NewChrIx,ans] = find(Mt<Mf & [zeros(1,Nsel);Mf(1:Nind-1,:)]<= Mt);%寻找满足条件的个体 返回返回数组 X 中每个元素的行和列下标
[ans,shuf] = sort(rand(Nsel,1));%生成Nsel*1的随机数矩阵  按升序对 A 的元素进行排序 返回选择的shuf 随机打乱个体顺序 保证后续遗传算子操作的随机性
NewChrIx = NewChrIx(shuf);%返回shuf索引的矩阵元素

7.适应度函数（其实根本不需要函数的。另外，如果要用链接2的代码，这里记得把fintness多了的n删了）

function FitnV = Fitness(len) %适应度函数
    %输入：
    %len 个体的长度（TSP的距离）
    %输出：
    %FitnV 个体的适应度值
FitnV = 1./len;

8.利用轮盘赌法，select选择优秀个体

function SelCh = Select(Chrom,FitnV,GGAP)
%输入：
%Chrom 种群
%FitnV 适应度值
%GGAP 选择概率
%输出：
%SelCh 被选择的个体
NIND = size(Chrom,1);%种群个体总数
NSel = max(floor(NIND * GGAP+.5),2);%确定下一代种群的个体数，如果不足二个就计为二个
ChrIx = Sus(FitnV,NSel);
SelCh = Chrom(ChrIx,:);

9.交叉，基因交换，产生新的路线（换两个城市的排序）

function [a,b] = intercross(a,b)
%输入：
%a和b为两个待交叉的个体
%输出：
%a和b为交叉后得到的两个个体
L = length(a);
%随机产生交叉区段
r1 = randsrc(1,1,[1:L]); % 这里先随机选出两个位置，位置不能超过
r2 = randsrc(1,1,[1:L]); % 
if r1~=r2
    a0 = a;
    b0 = b;
    s = min([r1,r2]);
    e = max([r1,r2]);
    for i = s:e
        a1 = a;
        b1 = b;
        %第一次互换
        a(i) = b0(i);
        b(i) = a0(i);
        %寻找相同的城市
        x = find(a==a(i)); % 如果有相同的城市，x会得到包含i的两个值，y同理
        y = find(b==b(i));
        %第二次互换产生新的解
        i1 = x(x~=i);      % 将位置不是i但相同的城市标记出来
        i2 = y(y~=i);
        if ~isempty(i1)
            a(i1)=a1(i);
        end
        if ~isempty(i2)
            b(i2)=b1(i);   % 注意，原文这里有误，应该是b(i2)
        end
    end
end
% 这里增加一段代码，r1=r2时，两个个体只在一点交叉
if r1 == r2
    a0 = a;
    b0 = b;
    a(r1) = b0(r1);
    b(r1) = a0(r1);
    x = find(a==a(r1));
    y = find(b==b(r1));
    i1 = x(x~=r1);
    i2 = y(y~=r1);
    if ~isempty(i1)
        a(i1) = a0(r1);
    end
    if ~isempty(i2)
        b(i2) = b0(r1);
    end
end

10.概率选择交叉的个体

function SelCh = Recombin(SelCh,Pc)
%交叉操作
%输入：
%SelCh 被选择的个体
%Pc    交叉概率
%输出：
%SelCh 交叉后的个体
 
NSel = size(SelCh,1);
for i = 1:2:NSel - mod(NSel,2) % 若不减去余数且NSel如果是奇数，最后一个i会等于NSel，i+1报错
    if Pc>=rand %交叉概率PC
        [SelCh(i,:),SelCh(i+1,:)] = intercross(SelCh(i,:),SelCh(i+1,:));
    end
end

11.变异（个体一定几率进行变异）

function SelCh = Mutate(SelCh,D,Pm)
%变异操作
%输入：
%SelCh  被选择的个体
%Pm  变异概率
%输出：
%SelCh  变异后的个体
index = SelCh;
col = size(SelCh,2); %返回SelCh的列数
%lenSelCh = [];
lenSelCh = PathLength(D,SelCh);
[NSel,L] = size(SelCh);
for i = 1:NSel
    if Pm >= rand
        R = randperm(L);
        index(i,R(1:2)) = index(i,R(2:-1:1)); % 将个体i中R(1)和R(2)这两个位置的城市互换
        p = [index(i,:) index(i,1)];
        i1 = p(1:end-1);
        i2 = p(2:end);
        DIndexi = sum(D((i1-1)*col+i2));              % 计算出变异后个体的路径距离
        if DIndexi < lenSelCh(i)                 % 如果变异后路径距离比变异前更小，则保留变异
            SelCh(i) = index(i);
        end
    end
end

12.进化逆转（我是没搞懂，不过没出错）

function SelCh = Reverse(SelCh,D)
%进化逆转函数
%输入：
%SelCh  被选择的个体
%D  各城市的距离矩阵
%输出：
%SelCh  进化逆转后的个体
 
[row,col] = size(SelCh);
ObjV = PathLength(D,SelCh);
SelCh1 = SelCh;
for i = 1:row
    r1 = randsrc(1,1,[1:col]);
    r2 = randsrc(1,1,[1:col]);
    mininverse = min([r1 r2]);
    maxinverse = max([r1 r2]);
    SelCh1(i,mininverse:maxinverse) = SelCh1(i,maxinverse:-1:mininverse);
end
ObjV1 = PathLength(D,SelCh1);%计算路线长度
index = ObjV1<ObjV;
SelCh(index,:)=SelCh1(index,:);

13.更新种群（根据适应度迭代出新种群）

function Chrom = Reins(Chrom,SelCh,ObjV)
%重插入子代的种群
%输入：
%Chrom      父代的种群
%SelCh      子代的种群
%ObjV       父代适应度
%输出：
%Chrom      组合父代与子代后得到的新种群
NIND = size(Chrom,1);
NSel = size(SelCh,1);
[TobjV,index] = sort(ObjV);
Chrom =  [Chrom(index(1:NIND-NSel),:);SelCh];

14.最后漏了一个dsxy2figxy函数，把坐标值转换到0,1之间满足annotation函数需要的。没太看懂这个怎么做的，不过能用就行

function varargout = dsxy2figxy(varargin) %%将两点坐标值转化至0,1之间
if length(varargin{1}) == 1 && ishandle(varargin{1}) ...
                            && strcmp(get(varargin{1},'type'),'axes')   
    hAx = varargin{1};
    varargin = varargin(2:end);
else
    hAx = gca;
end
if length(varargin) == 1
    pos = varargin{1};
else
    [x,y] = deal(varargin{:});
end
axun = get(hAx,'Units');
set(hAx,'Units','normalized'); 
axpos = get(hAx,'Position');
axlim = axis(hAx);
axwidth = diff(axlim(1:2));
axheight = diff(axlim(3:4));
if exist('x','var')
    varargout{1} = (x - axlim(1)) * axpos(3) / axwidth + axpos(1);
    varargout{2} = (y - axlim(3)) * axpos(4) / axheight + axpos(2);
else
    pos(1) = (pos(1) - axlim(1)) / axwidth * axpos(3) + axpos(1);
    pos(2) = (pos(2) - axlim(3)) / axheight * axpos(4) + axpos(2);
    pos(3) = pos(3) * axpos(3) / axwidth;
    pos(4) = pos(4) * axpos(4 )/ axheight;
    varargout{1} = pos;
end
set(hAx,'Units',axun)