%遗传算法解决TSP问题（旅行商）
C=[1304 2312;3639 1315;4177 2244;3712 1399;3488 1535;3325 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]; %31个省会坐标
N=size(C,1); %TSP问题的规模，即城市数目
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,:); %存储最优种群
len=zeros(NP,1); %存储路径长度
fitness=zeros(NP,1); %存储归一化适应值
gen=0;
%遗传算法循环
while gen<G
%计算路径长度
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<NP
nnper=randperm(nn);
A=F(nnper(1),:);
B=F(nnper(2),:); %以一对染色体为基础方便交叉操作
%交叉操作
W=ceil(N/10); %交叉点个数
P=unidrnd(N-W+1); %保证交叉点数为W，交叉范围为（P，P+W）
for i=1:W
x=find(A==B(P+i-1));
y=find(B==A(P+i-1)); %找到一堆染色体交叉点相同的位置（防止两次出现同一个城市）
temp=A(P+i-1);
A(P+i-1)=B(P+i-1);
B(P+i-1)=temp; %A B染色体在P+i-1基因数进行交叉
temp=A(x);
A(x)=B(y);
B(y)=temp; %防止出现两次同一个城市，其实质是将A染色体中的A（P+i-1）与和B（P+i-1）数相同的A（x）进行交换位置
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)=A(P2);
A(P2)=temp; %A染色体基因变异
temp=B(P1);
B(P1)=B(P2);
B(P2)=temp; %B染色体基因变异
F=[F;A;B];
[aa,bb]=size(F);
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(1)
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(2)
plot(Rlength)
xlabel('迭代次数');
ylabel('目标函数值');
title('适应度变化曲线');
