一、实践目的

运用粒子群优化算法解决TSP问题。

二、实践内容

1、问题定义
旅行商问题，即TSP问题（Traveling Salesman Problem）又为旅行推销员问题、货郎担问题，是数学领域中著名问题之一。假设有一个旅行商人要拜访n个城市，他必须选择所要走的路径，路径的限制是每个城市只能拜访一次，而且最后要回到原来出发的城市。路径的选择目标是要求得的路径路程为所有路径之中的最小值。
2、实验数据

3、实验要求
（1）输出最终优化路径；输出搜索过程中路径变化曲线（初始，第100次，第200次，最终）。
（2）讨论遗传算法中各类参数对搜索性能的影响。

4、粒子群优化算法基本介绍
（1）起源
粒子群优化算法(Particle Swarm optimization,PSO)又翻译为粒子群算法、微粒群算法、或微粒群优化算法。是通过模拟鸟群觅食行为而发展起来的一种基于群体协作的随机搜索算法。通常认为它是群集智能 (Swarm intelligence, SI) 的一种。它可以被纳入多主体优化系统(Multiagent Optimization System, MAOS)。粒子群优化算法是由Eberhart博士和kennedy博士发明。
（2）基本思想
粒子群算法的思想源于对鸟群觅食行为的研究，鸟群通过集体的信息共享使群体找到最优的目的地。如下图1，设想这样一个场景：鸟群在森林中随机搜索食物，它们想要找到食物量最多的位置。但是所有的鸟都不知道食物具体在哪个位置，只能感受到食物大概在哪个方向。每只鸟沿着自己判定的方向进行搜索，并在搜索的过程中记录自己曾经找到过食物且量最多的位置，同时所有的鸟都共享自己每一次发现食物的位置以及食物的量，这样鸟群就知道当前在哪个位置食物的量最多。在搜索的过程中每只鸟都会根据自己记忆中食物量最多的位置和当前鸟群记录的食物量最多的位置调整自己接下来搜索的方向。鸟群经过一段时间的搜索后就可以找到森林中哪个位置的食物量最多（全局最优解）。
3、实验要求
（1）输出最终优化路径；输出搜索过程中路径变化曲线（初始，第100次，第200次，最终）。
（2）讨论遗传算法中各类参数对搜索性能的影响。

4、粒子群优化算法基本介绍
（1）起源
粒子群优化算法(Particle Swarm optimization,PSO)又翻译为粒子群算法、微粒群算法、或微粒群优化算法。是通过模拟鸟群觅食行为而发展起来的一种基于群体协作的随机搜索算法。通常认为它是群集智能 (Swarm intelligence, SI) 的一种。它可以被纳入多主体优化系统(Multiagent Optimization System, MAOS)。粒子群优化算法是由Eberhart博士和kennedy博士发明。
（2）基本思想
粒子群算法的思想源于对鸟群觅食行为的研究，鸟群通过集体的信息共享使群体找到最优的目的地。如下图1，设想这样一个场景：鸟群在森林中随机搜索食物，它们想要找到食物量最多的位置。但是所有的鸟都不知道食物具体在哪个位置，只能感受到食物大概在哪个方向。每只鸟沿着自己判定的方向进行搜索，并在搜索的过程中记录自己曾经找到过食物且量最多的位置，同时所有的鸟都共享自己每一次发现食物的位置以及食物的量，这样鸟群就知道当前在哪个位置食物的量最多。在搜索的过程中每只鸟都会根据自己记忆中食物量最多的位置和当前鸟群记录的食物量最多的位置调整自己接下来搜索的方向。鸟群经过一段时间的搜索后就可以找到森林中哪个位置的食物量最多（全局最优解）。

（3）求解思路
粒子群优化算法（PSO），粒子群中的每一个粒子都代表一个问题的可能解，通过粒子个体的简单行为，群体内的信息交互实现问题求解的智能性。
在TSP问题中，我们将每一条访问城市的顺序编码为一个个体，每个种群有n个个体，即有n种访问顺序，同时，每个个体又有13个染色体，即[2,13]的随机排列（城市1作为起始和终止城市，不进行粒子操作），代表访问城市的顺序。 通过每一代的演化，对粒子群进行位置、速度更新操作，选择合适个体（最优的顺序）。
编码：位置（巡回顺序）：符号编码。速度编码：定义为交换子序列。单个交换子 si ={n，m}对个体的第n个和第m个元素进行交换，则速度表示为交换子序列ss = {s1，s2…sm}由m个交换子组成，同时交换子序列有顺序，从第一个交换子操作。
适应度：为了能够让适应度高的个体保存下来，定义适应度为：

其中 distance 为每个个体（路径）的距离。显然，如果距离最短，则适应度最高，更利于遗传给后代。
速度更新：

三、实践结果

1、代码实现
使用matlab语言编写程序：

city =[16.47 96.10;
    16.47 94.44;
    20.09 92.54;
    22.39 93.37; 
    25.23 97.24;
22.00 96.05; 
20.47 97.02; 
17.20 96.29;
16.30 97.38; 
14.05 98.12; 
16.53 97.38; 
21.52 95.59; 
19.41 97.13; 
20.09 94.55;];
citynum = size(city,1);               %城市数量
dist_city = zeros(citynum,citynum);   %初始化城市距离
for i = 1:citynum                     %计算城市之间的距离
    for j = 1:citynum                 
        link = (city(j,1)-city(i,1)).^2+(city(j,2)-city(i,2)).^2;
        dist_city(i,j)=sqrt(link);
    end
end
padai =  2000;                   %演化代数
pasize = 100;                   %粒子数目
padim = citynum-1;              %维度
pos = initpos( pasize,padim );  %初始化位置
v = initv( pasize,padim );%初始化速度
pbest =zeros(pasize,1); %单个粒子最优适应度
pid =zeros(size(pos));  %单个粒子对应的位置
w = 0.9;                %权重
maxgbest = zeros(1);    %整个过程中最优适应度
maxpgd = zeros(1,size(pos,2)); %整个过程全局最优位置
maxfitvalueall = [];    %各代最优适应度
for item = 1:padai     %演化
    fitvalue = fit_cal(pos,dist_city); %适应度计算
    %%%%%%  更新个体最优和全局最优   %%%%%%
    [ pbest,pid,gbest,pgd ] = fit_cmp(pos,fitvalue,pbest,pid);
    v = updatev( v,w,pos,pid,pgd); %更新速度
    pos = updatepos( pos,v );      %更新位置
    w = w - (w-0.01)/padai;        %权重更新
    maxfitvalueall =[maxfitvalueall,maxgbest];
     if (gbest>maxgbest)           %保留最优粒子
         maxgbest = gbest;         %最优适应度
         maxpgd = pgd;             %最优粒子
     end
end
[fin_fit,max_fit_index] =max(fitvalue);
lu_bin = pos(max_fit_index,:); %适应度最大的路径
if(maxgbest>fin_fit)           %最优粒子
    lu_bin = pgd;
end
%%%%%%   计算最短路径        %%%%%%
min_distance = 0;
    for j = 1:length(lu_bin)-1
        min_distance =min_distance +dist_city(lu_bin(1,j),lu_bin(1,j+1));
    end
min_distance = min_distance+dist_city(1,lu_bin(1,1))+...
    dist_city(lu_bin(1,end),1);
disp(strcat('最短距离为：',num2str(min_distance))) %输出最短距离
%%%%%%  将首末位置加入绘图   %%%%%%
dit_lu = zeros(citynum+1,2);
dit_lu(1,:) =city(1,:); 
dit_lu(end,:) =city(1,:); 
for i =2: citynum
    dit_lu(i,:) = city(lu_bin(i-1),:);
end
figure(1)
hold on;grid on;         %保持图像 开启网格
plot(city(:,1),city(:,2),'r*','linewidth',5)  %绘制城市
plot(dit_lu(:,1),dit_lu(:,2),'linewidth',2)   %绘制最优路径
title('10城市TSP问题')                         %标题
for i =1:citynum                              %添加城市名
    text(city(i,1)+0.1,city(i,2),strcat('城市',num2str(i)))
end
figure(2)
plot(maxfitvalueall)                         %绘制
grid on; xlabel('迭代次数/n')                 %开启网格
ylabel('适应度'); title('演化过程最大适应度')  %标题 左边轴名称

function [ newv ] = fiel(v , r)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%fiel(): 筛选欠r个交换子 r筛选比例
% v newv 筛选前后的速度
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
index = length(v)/2;  %获取交换子个数
len = round(index*r); %计算筛选个数
newv = v(1:2*len);    %保留前len个
end
function [ jhz ] = find_jhz(pbest,x)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%find_jhz寻求基本交换子 由x->pbest
% jhz 计算得到的交换子
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
len = length(pbest);    %获取长度
jhz = ones((len-1)*2,1);%初始化交互子
for i =1:len-1          %循环计算
    index = find(x == pbest(i));
    if i~=index         %是否需要交换
       jhz(2*i-1) = i;
       jhz(2*i) = index;
    end
    tt = x(index);      %交换两个位置值
    x(index) =x(i);
    x(i) = tt;
end
end
function [ fitvalue ] = fit_cal(pos,dist_city)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% fit_cal() 适应度计算
% pos 粒子群 dist_city 城市距离
% fitvalue = 2/distance;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[sizex,sizey] = size(pos);   %获取粒子群维度
newgen = ones(sizex,sizey+2);%初始化
newgen(:,2:end-1)=pos;       %加入首尾城市
fitvalue = zeros(sizex,1);   %初始化适应度
for i = 1:sizex    %循环计算每一个种群的适应度
    dist = 0;
    for j = 1:size(newgen,2)-1 %迭代计算距离
        dist =dist +dist_city(newgen(i,j),newgen(i,j+1));
    end
    fitvalue(i) = dist;        %保存距离值
end          
fitvalue = 2*ones(size(fitvalue))./fitvalue;%计算适应度
end
function [ pbest,pid,gbest,pgd ] = fit_cmp(pos,fitvalue,pbest,pid)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% fit_cmp（） 更新个体和全局最优粒子和适应度
% pos当前位置 fitvalue 当前适应度
% pbest 单个粒子最优适应度 pid 单个粒子对应的位置
% gbest 全局最优适应度 pgd 全局最优位置
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[max1,index1] = max(fitvalue);%找到最大适应度
gbest = max1;                 %最优适应度
pgd = pos(index1,:);          %最优粒子
for i = 1:size(pos,1)         %循环个体
    if(fitvalue(i)>pbest(i))  %找到最优适应度
        pid(i,:) = pos(i,:);  %最优适应度
        pbest(i) = fitvalue(i);%个体最优位置
    end
end
end
function [ pos ] = initpos( pasize,padim )
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% initpos(）初始化粒子群
% pasize粒子群大小 padim维度
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
pos =zeros(pasize,padim);         %初始化
for i = 1:pasize
    t = 2:padim+1;                %生成初始访问
    ranorder = randperm(padim);   %乱序排列
    pos(i,:)=t(ranorder);         %生成粒子
end
end
function [ v ] = initv( pasize,padim ) 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% initv( pasize,padim ) 初始化速度
% pasize粒子群大小 padim维度
%每行两两元素为一个交换子
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
v =floor(1+(padim-1)*rand(pasize,(padim-1)*2));
end
function [ pos ] = updatepos( pos,v )
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% updatepos( pos,v )对粒子的位置更新
% pos：位置    v：速度
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[vx,vy] = size(v);     %粒子群维度和大小
for i = 1:vx           %循环每个粒子
    for j =1:vy/2      %对每个交换子运算
        tt = pos(i,v(i,2*j-1)); 
        pos(i,v(i,2*j-1)) = pos(i,v(i,2*j));
        pos(i,v(i,2*j))=tt;  %更新位置
    end
end
end
function [ newv ] = updatev( v,w,pos,pid,pgd)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% updatev()对粒子的速度更新
% pos：位置    v：速度  w：权重
%pid 个体最优  pgd 全局最优
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[posx,posy] = size(pos);    %回去维度
v1=[];v2=[];v3=[];          %定义三个分量
v22=zeros(posx,(posy-1)*2); %初始化
v33=zeros(posx,(posy-1)*2);
% r1= round(1);r2=round(1);
r1= 0.7;r2=0.7;             %学习因子
for i = 1:posx              %计算当前位置与最优位置的差
    v22(i,:) = find_jhz(pid(i,:),pos(i,:));
    v33(i,:) = find_jhz(pgd,pos(i,:));
end
for i = 1:posx              %循环计算三个分量
    v1 = [v1;fiel(v(i,:),w)]; 
    v2 = [v2;fiel(v22(i,:),r1)];
    v3 = [v3;fiel(v33(i,:),r2)];
end
    newv = [v1 v2 v3 ];     %合并三个分量 更新速度
end

2、运算结果
（1）粒子群大小n=200，演化代数padai=100，最短距离为42.0138，路线与演化最大适应度如下图2所示。


（2）粒子群大小n=200，演化代数padai=600，最短距离为43.4401，路线与演化最大适应度如下图3所示。


（3）粒子群大小n=100，演化代数padai=100，最短距离为47.1552，路线与演化最大适应度如下图4所示。


（4）粒子群大小n=100，演化代数padai=600，最短距离为44.6518，路线与演化最大适应度如下图5所示。

四、实践总结

从上述结果可知，针对14个城市的TSP问题，粒子群优化算法求解得出的最短路径值为42.0138。
粒子群算法是基于概率的随机自搜索算法，同遗传算法一样都是不稳定的，每次的搜索结果都不尽相同。需要不断修改粒子群大小与演化次数来获取最优值。 但是，每次运行的结果不会相差很大，最短距离基本都在42-45之间，上面仅给出了四次结果。粒子群算法在计算最优解会出现陷入局部最优的情况，但是它的运行效率高，在实际的应用中有很好的效果。