【建模算法】基于粒子群算法（PSO）求解TSP问题（matlab求解）
TSP (traveling salesman problem，旅行商问题)是典型的NP完全问题，即其最坏情况下的时间复杂度随着问题规模的增大按指数方式增长，到目前为止还未找到一个多项式时间的有效算法。本文探讨了使用matlab软件，基于粒子群算法求解TSP问题。
一、问题描述
本案例以31个城市为例，假定31个城市的位置坐标如表1所列。寻找出一条最短的遍历31个城市的路径。
| 城市编号 | X坐标 | Y坐标 | 城市编号 | X坐标 | Y坐标 |
|---|---|---|---|---|---|
| 1 | 1.304 | 2.312 | 17 | 3.918 | 2.179 |
| 2 | 3.639 | 1.315 | 18 | 4.061 | 2.37 |
| 3 | 4.177 | 2.244 | 19 | 3.78 | 2.212 |
| 4 | 3.712 | 1.399 | 20 | 3.676 | 2.578 |
| 5 | 3.488 | 1.535 | 21 | 4.029 | 2.838 |
| 6 | 3.326 | 1.556 | 22 | 4.263 | 2.931 |
| 7 | 3.238 | 1.229 | 23 | 3.429 | 1.908 |
| 8 | 4.196 | 1.044 | 24 | 3.507 | 2.376 |
| 9 | 4.312 | 0.79 | 25 | 3.394 | 2.643 |
| 10 | 4.386 | 0.57 | 26 | 3.439 | 3.201 |
| 11 | 3.007 | 1.97 | 27 | 2.935 | 3.24 |
| 12 | 2.562 | 1.756 | 28 | 3.14 | 3.55 |
| 13 | 2.788 | 1.491 | 29 | 2.545 | 2.357 |
| 14 | 2.381 | 1.676 | 30 | 2.778 | 2.826 |
| 15 | 1.332 | 0.695 | 31 | 2.37 | 2.975 |
| 16 | 3.715 | 1.678 |
二、粒子群算法（PSO）介绍
1.粒子群算法的由来及思想
粒子群算法最早是由两名美国的科学家基于群鸟觅食，寻找最佳觅食区域的过程所提出来的，作为一种智能算法，PSO模拟的就是最佳决策的过程，鸟群觅食类似于人类的决策过程，想想在你做出选择之前，是不是会受到自己的经验（局部最优）以及周围人的经历（全局最优）的影响？同样的道理，群鸟在觅食的过程当中，每一只鸟的初始位置都处于随机状态，当然也不知道最佳的觅食点在何处，并且每只鸟的飞行方向也是随机的。可以认为，在觅食的初期，鸟群的运动轨迹都是杂乱无章的，随着时间的推移，处于随机位置的鸟类通过群内的相互学习、共享觅食信息，每一只鸟在每一次的觅食过程中结合自己的经验和同伴传送的信息估计目前所处的位置能够找到食物有多大的价值。基于这样的搜索方式，粒子群算法（Particle Swarm Optimization,PSO）应运而生。
2.基本概念
每只鸟在某处位置能够找到食物的可能可以通过适应值来刻画，每只鸟能够记住自己的觅食位置，并且找到其中的最佳位置（局部最优），鸟群中所有个体的最佳位置就可以看做整个鸟群的最佳觅食点（全局最优）。可以预见的是，整个鸟群的觅食活动总体上一定是往这个全局最优的觅食区域运动的，通过鸟群觅食位置的不断移动即不断地迭代，速度的不断更新，鸟群往该最优位置步步逼近。
迭代终止条件根据具体问题而定，一般达到预定最大迭代次数或者粒子群目前为止搜寻到的最优位置满足目标函数的最小容许误差。
3.PSO算法的主要实现步骤
- 初始化粒子群。包括粒子的初始位置及速度，惯性因子等参数值
粒子数M一般选取20~40个，不过对于一些特殊的难题需要更多的粒子，粒子数量越多，搜索范围就越广，越容易找到全局最优解，但是也会带来更大的计算消耗。
-
评价各个粒子的初始适应值。
-
将初始的适应值作为各个粒子的局部最优解，保存各粒子的最优位置。并找到其中的最优值，作为全局最优解的初值，并记录其位置
-
更新粒子速度及位置
-
计算更新后粒子的适应值，更新每个粒子的局部最优值以及整个粒子群的全局最优值。
-
重复4~5直至满足迭代结束条件
流程图如下：
三、粒子群算法求解TSP问题
在TSP问题中，粒子的位置可以使用路径来表示，速度如何表示却是一个难题。基本粒子群算法的速度和位置更新公式不再适用，我们重新定义了速度和位置的更新公式。
基本粒子群算法速度位置更新公式：
本文重新定义了速度，以及三个运算规则：位置-位置，常数*速度，位置+速度
1.速度定义为交换序列，位置1+速度=位置2定义为位置1通过速度这个交换序列变成位置2，速度中每一个索引对应的数字意为：位置1中索引为该数字的城市与位置1中索引为速度索引的城市对调，例如：
位置1=[5 4 3 2 1]
位置2=[1 2 3 4 5]
位置2[1]=位置1[5]
因此位置1[1]与位置1[5]互换位置，速度[1]=5
交换位置之后：
位置1=[1 4 3 2 5]
位置2=[1 2 3 4 5]
位置2[2]=位置1[4]
因此位置1[2]与位置1[4]互换位置，速度[2]=4
交换位置之后：
位置1=[1 2 3 4 5]
位置2=[1 2 3 4 5]
位置2[3]=位置1[3]
因此位置1[3]与位置1[3]互换位置，速度[3]=3
由此可以推导，速度=[5 4 3 4 5]
位置2-位置1即为速度
2.常数*位置定义为：速度中的每一个值以该常数的概率保留，避免过早陷入局部最优解。
四、运行结果
优化结果
优化后的路线图：
迭代过程图：
最优解：
25->30->31->27->28->26->20->24->19->17->3->18->22->21->23->6->5->16->4->8->9->10->2->7->13->12->14->15->1->29->11->25
总距离：16.711
当然，PSO是一个基于概率的随机自搜索算法，每次的搜索结果都不尽相同，但是若算法的收敛性表现的较好，也可以认为对算法的设计是合理的。
五、matlab源代码
主函数my_main.m如下：
%% 粒子群算法求解TSP问题
clear,clc;
%% 1.读入数据
data=importdata('p_xy.xlsx');
citys=data.data.Sheet1(:,2:3); %各城市坐标数据
%% 2.计算距离矩阵
n=size(citys,1);
D=zeros(n,n);
for i=1:n
for j=i+1:n
D(i,j)=sqrt(sum((citys(i,:)-citys(j,:)).^2));
D(j,i)=D(i,j);
end
end
%% 3.初始化参数
c1=0.1; %个体学习因子
c2=0.075; %社会学习因子
w=1; %惯性因子
m=500; %粒子数量
pop=zeros(m,n); %粒子位置
v=zeros(m,n); %粒子速度
gen=1; %迭代计数器
genmax=2000; %迭代次数
fitness=zeros(m,1); %适应度函数值
Pbest=zeros(m,n); %个体极值路径
Pbest_fitness=zeros(m,1); %个体极值
Gbest=zeros(genmax,n); %群体极值路径
Gbest_fitness=zeros(genmax,1); %群体极值
Length_ave=zeros(genmax,1); %各代路径的平均长度
ws=1; %惯性因子最大值
we=0.8; %惯性因子最小值
%% 4.产生初始粒子
% 4.1随机产生粒子初始位置和速度
for i=1:m
pop(i,:)=randperm(n);
v(i,:)=randperm(n);
end
% 4.2计算粒子适应度函数值
for i=1:m
for j=1:n-1
fitness(i)=fitness(i) + D(pop(i,j),pop(i,j+1));
end
fitness(i)=fitness(i) + D(pop(i,end),pop(i,1));
end
% 4.3计算个体极值和群体极值
Pbest_fitness=fitness;
Pbest=pop;
[Gbest_fitness(1),min_index]=min(fitness);
Gbest(1,:)=pop(min_index,:);
Length_ave(1)=mean(fitness);
%% 5.迭代寻优
while gen<genmax
% 5.1更新迭代次数与惯性因子
gen=gen+1;
w = ws - (ws-we)*(gen/genmax)^2;
% 5.2更新速度
%个体极值修正部分
change1=position_minus_position(Pbest,pop);
change1=constant_times_velocity(c1,change1);
%群体极值修正部分
change2=position_minus_position(repmat(Gbest(gen-1,:),m,1),pop);
change2=constant_times_velocity(c2,change2);
%原速度部分
v=constant_times_velocity(w,v);
%修正速度
for i=1:m
for j=1:n
if change1(i,j)~=0
v(i,j)=change1(i,j);
end
if change2(i,j)~=0
v(i,j)=change2(i,j);
end
end
end
% 5.3更新位置
pop=position_plus_velocity(pop,v);
% 5.4适应度函数值更新
fitness=zeros(m,1);
for i=1:m
for j=1:n-1
fitness(i)=fitness(i) + D(pop(i,j),pop(i,j+1));
end
fitness(i)=fitness(i) + D(pop(i,end),pop(i,1));
end
% 5.5个体极值与群体极值更新
for i=1:m
if fitness(i)<Pbest_fitness(i)
Pbest_fitness(i)=fitness(i);
Pbest(i,:)=pop(i,:);
end
end
[minvalue,min_index]=min(fitness);
if minvalue<Gbest_fitness(gen-1)
Gbest_fitness(gen)=minvalue;
Gbest(gen,:)=pop(min_index,:);
else
Gbest_fitness(gen)=Gbest_fitness(gen-1);
Gbest(gen,:)=Gbest(gen-1,:);
end
Length_ave(gen)=mean(fitness);
end
%% 6.结果显示
[Shortest_Length,index] = min(Gbest_fitness);
Shortest_Route = Gbest(index,:);
disp('最优解：')
p=OutputPath(Shortest_Route);
disp(['总距离：',num2str(Shortest_Length)]);
disp('-------------------------------------------------------------')
%% 7.绘图
%画出最优解的路线图
DrawPath(Shortest_Route, citys);
figure;
plot(1:genmax,Gbest_fitness,'b',1:genmax,Length_ave,'r:')
legend('最短距离','平均距离')
xlabel('迭代次数')
ylabel('总距离')
title('各代最短距离与平均距离对比')
函数position_minus_position.m文件：
function change=position_minus_position(best,pop)
%记录将pop变成best的交换序列
for i=1:size(best,1)
for j=1:size(best,2)
change(i,j)=find(pop(i,:)==best(i,j));
temp=pop(i,j);
pop(i,j)=pop(i,change(i,j));
pop(i,change(i,j))=temp;
end
end
end
函数constant_times_velocity.m文件：
function change = constant_times_velocity(constant,change)
% 以一定概率保留交换序列
for i=1:size(change,1)
for j=1:size(change,2)
if rand>constant
change(i,j)=0;
end
end
end
end
函数position_plus_velocity.m文件：
function pop = position_plus_velocity(pop,v)
%利用速度记录的交换序列进行位置修正
for i=1:size(pop,1)
for j=1:size(pop,2)
if v(i,j)~=0
temp=pop(i,j);
pop(i,j)=pop(i,v(i,j));
pop(i,v(i,j))=temp;
end
end
end
end
函数position_plus_velocity.m文件：
function pop = position_plus_velocity(pop,v)
%利用速度记录的交换序列进行位置修正
for i=1:size(pop,1)
for j=1:size(pop,2)
if v(i,j)~=0
temp=pop(i,j);
pop(i,j)=pop(i,v(i,j));
pop(i,v(i,j))=temp;
end
end
end
end
完整源代码下载：https://download.csdn.net/download/baidu/85284316
1090
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
