|
???????????设计八：粒子群算法优化求解
旅行商（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/80367158
2、算法生物学背景
??粒子群算法属于群智能算法的一种，使用过模拟鸟群捕食行为设计的。假设区域里只有一块食物，鸟群的任务是找到这块食物。鸟群在整个搜寻过程中，通过相互传递各自的信息，让其他的鸟知道自己的位置，通过这样的协作，来判断自己找到的是不是最优解，同时，也将最优解的信息传递给整个鸟群。
??最终，整个鸟群都能聚集在食物源周围，即找到了问题的最优解，即问题收敛。
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<genmax
% 6.1更新迭代次数与惯性因子
gen=gen+1;
w = ws - (ws-we)*(gen/genmax)^2;
% 6.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
% 6.3更新位置
pop=position_plus_velocity(pop,v);
% 6.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
% 6.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
%% 7.结果显示
[Shortest_Length,index] = min(Gbest_fitness);
Shortest_Route = Gbest(index,:);
disp(['最短距离:' num2str(Shortest_Length)]);
disp(['最短路径:' num2str([Shortest_Route Shortest_Route(1)])]);
toc;
%% 8.绘制历史最优图
for b=1:2000
if (b == 1 || b == 5 || b==20 || b == 100 || b == 500||b==2000) %改b的值保存图片就是对应的迭代次数结果
% figure;
a=Gbest(b,:);
plot([citys(a,1);citys(a(1),1)],[citys(a,2);citys(a(1),2)],'o-');
xlabel('城市位置横坐标x');
ylabel('城市位置纵坐标y');
title(['迭代次数:',num2str(b),' 优化最短距离:',num2str(Gbest_fitness(b))]);
filename=num2str(b);
fileformat='png';
saveas(gcf,filename,fileformat);
end
end
plot([citys(Shortest_Route,1);citys(Shortest_Route(1),1)],...
[citys(Shortest_Route,2);citys(Shortest_Route(1),2)],'o-');
title(['迭代次数:',int2str(gen),' 优化最短距离:',num2str(Shortest_Length)]);
% title(['粒子群算法优化路径(最短距离:' num2str(Shortest_Length) ')'])
grid on
pause(0.005);
for i = 1:size(citys,1)
text(citys(i,1),citys(i,2),[' ' num2str(i)]);
end
text(citys(Shortest_Route(1),1),citys(Shortest_Route(1),2),' 起点');
text(citys(Shortest_Route(end),1),citys(Shortest_Route(end),2),' 终点');
xlabel('城市位置横坐标');
ylabel('城市位置纵坐标');
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
%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
|