设计八：粒子群算法优化求解
旅行商（TSP）问题，是数学领域中著名问题之一，被证明具有NPC计算复杂性。该问题假设有一个旅行商人要拜访n个城市，他必须选择所要走的路径，路径的限制是每个城市只能拜访一次，而且最后要回到原来出发的城市。路径的选择目标是要求得的路径路程为所有路径之中的最小值。从TSP问题测试数据集中选择一组测试数据，选择一种改进的粒子群算法，编写Python程序，求解最优路径问题。设计能够进行人机交互的可视化界面，使用图形方式表达求解过程，统计并分析实验数据，给出相应的数据表格和图形。学生也可以选择其他复杂问题进行求解。
1、算法参考学习https://finthon.com/python-pso/ https://www.omegaxyz.com/2018/01/12/python_pso/ https://blog.csdn.net/wang454592297/article/details/803671582、算法生物学背景
粒子群算法属于群智能算法的一种，使用过模拟鸟群捕食行为设计的。假设区域里只有一块食物，鸟群的任务是找到这块食物。鸟群在整个搜寻过程中，通过相互传递各自的信息，让其他的鸟知道自己的位置，通过这样的协作，来判断自己找到的是不是最优解，同时，也将最优解的信息传递给整个鸟群。
最终，整个鸟群都能聚集在食物源周围，即找到了问题的最优解，即问题收敛。
3、PSO算法基本原理概述：
从鸟群的视角：假设一个区域，所有的鸟都不知道食物的位置，但是它们知道当前位置离食物还有多远。
从PSO算法的视角:每个解看作一只鸟，称为“粒子(particle)”，所有的粒子都有一个适应值，每个粒子都有一个速度决定它们的飞翔方向和距离，粒子们追随当前最优粒子在解空间中搜索。初始化一群随机粒子，通过迭代找到最优。每次迭代中，粒子通过跟踪“个体极值（pbest）”和“全局极值(gbest)”来更新自己的位置。
4、算法基本步骤:
1）初始化 初始化粒子群（粒子群共有n个粒子）：给每个粒子赋予随机的初始位置和速度。 2）计算适应值 根据适应度函数，计算每个粒子的适应值。 3）求个体最佳适应值 对每一个粒子，将其当前位置的适应值与其历史最佳位置（pbest）对应的适应值比较， 如果当前位置的适应值更高，则用当前位置更新历史最佳位置。 4）求群体最佳适应值 对每一个粒子，将其当前位置的适应值与其全局最佳位置（gbest）对应的适应值比较， 如果当前位置的适应值更高，则用当前位置更新全局最佳位置。 5）更新粒子位置和速度 根据公式更新每个粒子的速度与位置。 6）判断算法是否结束 若未满足结束条件，则返回步骤2，若满足结束条件则算法结束，全局最佳位置（gbest） 即全局最优解。
5、算法matlab代码
%TSPPSO.m
%% 1.清空环境变量
close all;
clear all;
clc;
tic;
%% 2.导入数据
% %%load citys_data.mat; %数据集的变量名为citys
citys=[1304 2312;3639 1315;4177 2244;3712 1399;3488 1535;3326 1556;...
2328 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];
%citys=[1304 2312;3639 1315;4177 2244;3712 1399;3488 1535];
%% 3.计算城市间相互距离
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
%% 4.初始化参数
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; %惯性因子最小值
NC=1;
%% 5.产生初始粒子
% 5.1随机产生粒子初始位置和速度
for i=1:m
pop(i,:)=randperm(n);
v(i,:)=randperm(n);
end
% 5.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
% 5.3计算个体极值和群体极值
Pbest_fitness=fitness;
Pbest=pop;
[Gbest_fitness(1),min_index]=min(fitness);
Gbest(1,:)=pop(min_index,:);
Length_ave(1)=mean(fitness);
%% 6.迭代寻优
% figure(1);
while gen
%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
%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
%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
