遗传算法GA求解TSP问题

假设有一个旅行商人要拜访全国31个省会城市，它需要选择所要走的路径，路径的限制是每个城市只能拜访一次，而且最后要回到原来出发的城市。对路径选择的要求是：所选路径的路成为所有路径之中的最小值。
全国31个省会城市的坐标为

[1304 2312;3639 1315;4177 2244;3712 1399;3488 1535;3326 1556;3238 1229;4196 1044;

4312  790;4386  570;3007 1970;2562 1756;2788 1491;2381 1676;1332  695;3715 1678;

3918 2179;4061 2370;3780 2212;3676 2578;4029 2838;4263 2931;3429 1908;3507 2376;

3394 2643;3439 3201;2935 3240;3140 3550;2545 2357;2778 2826;2370 2975]

clear all;
close all;
clc;
C=[1304 2312;3639 1315;4177 2244;3712 1399;3488 1535;3326 1556;...
    3238 1229;4196 1044;4312  790;4386  570;3007 1970;2562 1756;...
    2788 1491;2381 1676;1332  695;3715 1678;3918 2179;4061 2370;...
    3780 2212;3676 2578;4029 2838;4263 2931;3429 1908;3507 2376;...
    3394 2643;3439 3201;2935 3240;3140 3550;2545 2357;2778 2826;...
    2370 2975]; 
N=size(C,1);%1表示按行计算，2表示按列计算
D=zeros(N);%城市距离矩阵
%%%%%求任意两城市距离
for i=1:N
    for j=1:N
        D(i,j)=((C(i,1)-C(j,1))^2+(C(i,2)-C(j,2))^2)^0.5;
    end
end
NP=200;%种群规模
G=1000;%遗传次数
f=zeros(NP,N);%用于存储种群
F=[];%种群更新中间存储（选择，交叉，遗传
for i=1:NP
    f(i,:)=randperm(N);%初始化种群
end
R=f(1,:);           %R存储历代最优路径，初始化时取种群第一个个体做最优
len=zeros(NP,1);    %路径长度
fitness=zeros(NP,1);%种群归一化后的适应度值
gen=0;
while gen
    %%%%%%%%%%计算路径长度
    for i=1:NP
        len(i,1)=D(f(i,N),f(i,1));
        for j=1:N-1
            len(i,1)=len(i,1)+D(f(i,j),f(i,j+1));
        end
    end
    maxlen=max(len);
    minlen=min(len);
    %%%%%%%%%%更新最短路径
    rr=find(len==minlen);
    R=f(rr(1,1),:);
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%归一化适应值
    for i=1:length(len)
        fitness(i,1)=((1-(len(i,1)-minlen)/(maxlen-minlen+0.001)));
    end
    %%%%%选择操作
    nn=0;
    for i=1:NP
        if fitness(i,1)>=rand%遍历每个个体，适应度大于某个随机值时选择
            nn=nn+1;
            F(nn,:)=f(i,:);
        end
    end
    [aa,bb]=size(F);
    while aa
        nnper=randperm(nn);%任选两个进行交叉和变异操作
        A=F(nnper(1),:);
        B=F(nnper(2),:);
        %%%%交叉操作
        W=ceil(N/10);%交叉点个数
        p=unidrnd(N-W+1);   %交叉点起始位置，把p到p+w位置上的元素交换
        for i=1:W
            x=find(A==B(p+i-1));
            y=find(B==A(p+i-1));
            temp=B(p+i-1);%A与B交换p+i-1的位置
            B(p+i-1)=A(p+i-1);
            A(p+i-1)=temp;
            temp=A(x);%A与B交换x和y的位置
            A(x)=B(y);
            B(y)=temp;
        end
        %%%%%%%%%%%%%%%%%%%变异操作
        p1=floor(1+N*rand());
        p2=floor(1+N*rand());
        while p1==p2
            p1=floor(1+N*rand());
            p2=floor(1+N*rand());
        end
        temp=A(p1);%将A 的p1和p2位置交换
        A(p1)=A(p2);
        A(p2)=temp;
        temp=B(p1);%将B的p1和p2的位置交换
        B(p1)=B(p2);
        B(p2)=temp;
        F=[F;A;B];
       [aa,bb]=size(F);%a的值每次+2
    end
    %%%%%%%%后事处理
    if aa>NP
        F=F(1:NP,:);%保持种群规模为NP
    end
    f=F;
    f(1,:)=R;%将种群第一个个体赋值为上一代最优个体
    clear F;%清除种群更新中间存储
    gen=gen+1;
    Rlength(gen)=minlen;%记录历代最优路径长度
end
figure
for i=1:N-1
    plot([C(R(i),1),C(R(i+1),1)],[C(R(i),2),C(R(i+1),2)],'bo-');
    hold on;
end
plot([C(R(N),1),C(R(1),1)],[C(R(N),2),C(R(1),2)]);
title(['优化最短距离:',num2str(minlen)]);
figure
plot(Rlength);
xlabel('迭代次数');
ylabel('目标函数值');
title('适应度进化曲线')