运用遗传算法解决TSP问题

求解15个城市TSP的GA问题，两城市间距离自拟或随机生成

问题分析：

用图论的术语来说，假设有一个图g=(v,e)，其中v是顶点集，e是边集，设d=(dij)是由顶点i和顶点j之间的距离所组成的距离矩阵，旅行商问题就是求出一条通过所有顶点且每个顶点只通过一次的具有最短距离的回路。
这个问题可分为对称旅行商问题(dij=dji,任意i,j=1,2,3，…,n)和非对称旅行商问题(dij≠dji,任意i,j=1,2,3，…,n)。
若对于城市v={v1,v2,v3,…,vn}的一个访问顺序为t=(t1,t2,t3,…,ti,…,tn),其中ti∈v(i=1,2,3,…,n)，且记tn+1= t1，则旅行商问题的数学模型为：
min l=σd(t(i),t(i+1)) （i=1,…,n）
旅行商问题是一个典型的组合优化问题，并且是一个np难问题，其可能的路径数目与城市数目n是成指数型增长的，所以一般很难精确地求出其最优解，本文采用遗传算法求其近似解。
流程图

步骤：
运行主程序GA_TSP.m
自拟的城市间距离在CityPosition.mat中

程序：

GA_TSP.m

  %遗传算法求解TSP问题(为选择操作从新设计后程序)
    %输入：
    %D       距离矩阵
    %NIND    为种群个数
    %X       参数是中国15个城市的坐标(初始给定)
    %MAXGEN  为停止代数，遗传到第MAXGEN代时程序停止,MAXGEN的具体取值视问题的规模和耗费的时间而定
    %m       为适值淘汰加速指数,最好取为1,2,3,4,不宜太大
    %Pc      交叉概率
    %Pm      变异概率
    %输出：
    %R       为最短路径
    %Rlength 为路径长度
    clear
    clc
    close all
    %% 加载数据
    load CityPosition
    %X=data;
    D=Distanse(X);  %生成距离矩阵
    N=size(D,1);    %城市个数
    %% 遗传参数
    NIND=100;       %种群大小
    MAXGEN=200;     %最大遗传代数
    Pc=0.9;         %交叉概率
    Pm=0.05;        %变异概率
    GGAP=0.9;       %代沟
    %% 初始化种群
    Chrom=InitPop(NIND,N);
    %% 画出随机解的路径图
    DrawPath(Chrom(1,:),X)
    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])
    title('优化过程')
    xlabel('代数')
    ylabel('最优值')
    ObjV=PathLength(D,Chrom);  %计算路径长度
    preObjV=min(ObjV);
    while gen<MAXGEN
        %% 计算适应度
        ObjV=PathLength(D,Chrom);  %计算路径长度
        % fprintf('%d   %1.10f\n',gen,min(ObjV))
        line([gen-1,gen],[preObjV,min(ObjV)]);pause(0.0001)
        preObjV=min(ObjV);
        FitnV=Fitness(ObjV);
        %% 选择
        SelCh=Select(Chrom,FitnV,GGAP);
        %% 交叉操作
        SelCh=Recombin(SelCh,Pc);
        %% 变异
        SelCh=Mutate(SelCh,Pm);
        %% 逆转操作
        SelCh=Reverse(SelCh,D);
        %% 重插入子代的新种群
        Chrom=Reins(Chrom,SelCh,ObjV);
        %% 更新迭代次数
        gen=gen+1 ;
    end
    %% 画出最优解的路径图
    ObjV=PathLength(D,Chrom);  %计算路径长度
    [minObjV,minInd]=min(ObjV);
    DrawPath(Chrom(minInd(1),:),X)
    %% 输出最优解的路径和总距离
    disp('最优解:')
    p=OutputPath(Chrom(minInd(1),:));
    disp(['总距离：',num2str(ObjV(minInd(1)))]);
    disp('-------------------------------------------------------------'

)

Distant.m

%% 计算两两城市之间的距离
%输入 a  各城市的位置坐标
%输出 D  两两城市之间的距离
function D=Distanse(a)
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

DrawPath.m
%% 画路径函数
%输入
% Chrom  待画路径   
% X      各城市坐标位置
function DrawPath(Chrom,X)
R=[Chrom(1,:) Chrom(1,1)]; %一个随机解(个体)
figure;
hold on
plot(X(:,1),X(:,2),'o','color',[0.5,0.5,0.5])
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,:);
row=size(A,1);
for i=2:row
    [arrowx,arrowy] = dsxy2figxy(gca,A(i-1:i,1),A(i-1:i,2));%坐标转换
    annotation('textarrow',arrowx,arrowy,'HeadWidth',8,'color',[0,0,1]);
end
hold off
xlabel('横坐标')
ylabel('纵坐标')
title('轨迹图')
box on

dsxy2figxy.m

function varargout = dsxy2figxy(varargin)
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)

Fitness.m

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


InitPop.m
%% 初始化种群
%输入：
% NIND：种群大小
% N：   个体染色体长度（这里为城市的个数）  
%输出：
%初始种群
function Chrom=InitPop(NIND,N)
Chrom=zeros(NIND,N);%用于存储种群
for i=1:NIND
    Chrom(i,:)=randperm(N);%随机生成初始种群
end

Mutate.m

%% 变异操作
%输入：
%SelCh  被选择的个体
%Pm     变异概率
%输出：
% SelCh 变异后的个体
function SelCh=Mutate(SelCh,Pm)
[NSel,L]=size(SelCh);
for i=1:NSel
    if Pm>=rand
        R=randperm(L);
        SelCh(i,R(1:2))=SelCh(i,R(2:-1:1));
    end
end

OutputPath.m

%% 输出路径函数
%输入：R 路径
function p=OutputPath(R)
R=[R,R(1)];
N=length(R);
p=num2str(R(1));
for i=2:N
    p=[p,'—>',num2str(R(i))];
end
disp(p)

PathLength.m

%% 计算各个体的路径长度
% 输入：
% D     两两城市之间的距离
% Chrom 个体的轨迹
function len=PathLength(D,Chrom)
[row,col]=size(D);
NIND=size(Chrom,1);
len=zeros(NIND,1);
for i=1:NIND
    p=[Chrom(i,:) Chrom(i,1)];
    i1=p(1:end-1);
    i2=p(2:end);
    len(i,1)=sum(D((i1-1)*col+i2));
end

Recombin.m

%% 交叉操作
% 输入
%SelCh  被选择的个体
%Pc     交叉概率
%输出：
% SelCh 交叉后的个体
function SelCh=Recombin(SelCh,Pc)
NSel=size(SelCh,1);
for i=1:2:NSel-mod(NSel,2)
    if Pc>=rand %交叉概率Pc
        [SelCh(i,:),SelCh(i+1,:)]=intercross(SelCh(i,:),SelCh(i+1,:));
    end
end

%输入：
%a和b为两个待交叉的个体
%输出：
%a和b为交叉后得到的两个个体

function [a,b]=intercross(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));
        y=find(b==b(i));
        i1=x(x~=i);
        i2=y(y~=i);
        if ~isempty(i1)
            a(i1)=a1(i);
        end
        if ~isempty(i2)
            b(i2)=b1(i);
        end
    end
end



%
% %交叉算法采用部分匹配交叉%交叉算法采用部分匹配交叉
% function [a,b]=intercross(a,b)
% L=length(a);
% r1=ceil(rand*L);
% r2=ceil(rand*L);
% r1=4;r2=7;
% if r1~=r2
%     s=min([r1,r2]);
%     e=max([r1,r2]);
%     a1=a;b1=b;
%     a(s:e)=b1(s:e);
%     b(s:e)=a1(s:e);
%     for i=[setdiff(1:L,s:e)]
%         [tf, loc] = ismember(a(i),a(s:e));
%         if tf
%             a(i)=a1(loc+s-1);
%         end
%         [tf, loc]=ismember(b(i),b(s:e));
%         if tf
%             b(i)=b1(loc+s-1);
%         end
%     end
% end

Reins.m

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

Reverse.m
%% 进化逆转函数
%输入
%SelCh 被选择的个体
%D     个城市的距离矩阵
%输出
%SelCh  进化逆转后的个体
function SelCh=Reverse(SelCh,D)
[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,:);

Select.m

%% 选择操作
%输入
%Chrom 种群
%FitnV 适应度值
%GGAP：代沟
%输出
%SelCh  被选择的个体
function SelCh=Select(Chrom,FitnV,GGAP)
NIND=size(Chrom,1);
NSel=max(floor(NIND*GGAP+.5),2);
ChrIx=Sus(FitnV,NSel);
SelCh=Chrom(ChrIx,:);

Sus.m

% 输入:
%FitnV  个体的适应度值
%Nsel   被选择个体的数目
% 输出:
%NewChrIx  被选择个体的索引号
function NewChrIx = Sus(FitnV,Nsel)
[Nind,ans] = size(FitnV);
cumfit = cumsum(FitnV);
trials = cumfit(Nind) / Nsel * (rand + (0:Nsel-1)');
Mf = cumfit(:, ones(1, Nsel));
Mt = trials(:, ones(1, Nind))';
[NewChrIx, ans] = find(Mt < Mf & [ zeros(1, Nsel); Mf(1:Nind-1, :) ] <= Mt);
[ans, shuf] = sort(rand(Nsel, 1));
NewChrIx = NewChrIx(shuf);

test.m

clear;
clc
for i=1:100
    i
    [L(i),P{i}]=GA_TSP;
end
[a,index]=min(L);
disp('最优解:')
disp(P{index})
disp(['总距离：',num2str(a)]);
A=P{index};
B=L(index);
save P A B