粒子群算法

最新推荐文章于 2022-09-04 14:27:48 发布

小嵌同学

最新推荐文章于 2022-09-04 14:27:48 发布

阅读量431

点赞数

分类专栏：算法文章标签：算法机器学习

本文链接：https://blog.csdn.net/weixin_45842280/article/details/120281420

版权

算法专栏收录该内容

2 篇文章 0 订阅

订阅专栏

设计八：粒子群算法优化求解

旅行商（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