粒子群算法讲解以及MATLAB实现_粒子群matlab代码

系列文章目录

遗传算法原理以及matlab实现https://blog.csdn.net/stm8151k4/article/details/126653614

---

前言

  粒子群优化算法(PSO：Particle swarm optimization) 是一种进化计算技术（evolutionary computation）。源于对鸟群捕食的行为研究。粒子群优化算法的基本思想：是通过群体中个体之间的协作和信息共享来寻找最优解。

  PSO的优势：在于简单容易实现并且没有许多参数的调节。目前已被广泛应用于函数优化、神经网络训练、模糊系统控制以及其他遗传算法的应用领域。

---

一、粒子群算法的基本原理

  粒子群算法通过设计一种无质量的粒子来模拟鸟群中的鸟，粒子仅具有两个属性：速度和位置，速度代表移动的快慢以及方向（与粒子下次运动有关，也就是该粒子对应变量的更新），位置代表当前所处的位置（也就是当前时刻该粒子所对应的变量值）。每个粒子在搜索空间中单独的搜寻最优解，并将其记为当前个体极值，并将个体极值与整个粒子群里的其他粒子共享，找到最优的那个个体极值作为整个粒子群的当前全局最优解，粒子群中的所有粒子根据自己找到的当前个体极值和整个粒子群共享的当前全局最优解来调整自己的速度和位置。从而通过不断迭代，使得粒子群不断获得更优的群体最优解。

二、算法流程

1.流程图如

  整个流程图可以分为两部分，左边是初始化部分，右边是迭代寻优部分。

2.初始化部分

  初始化部分中，首先是初始化参数：

c1 = 1.49445;
c2 = 1.49445;
N = 50; %　种群规模
D = 2;
T = 100; %　迭代次数
wa = 0.9;
wi = 0.4;
Xmax = 4 ;
Xmin = -4 ;
Ymax = 6 ;
Ymin = 1 ;
Vmax = 0.5 ; %　最大飞行速度
Vmin = -0.5 ; %　最小飞行速度

  其次是初始化粒子群：在初始化粒子群中，首先是随机获取粒子群中的每个粒子，然后是对每个粒子进行计算确定其适应度，更新各粒子的个体极值，然后进行比较确定粒子群的最优极值。

for i=1:N
    Xm(i,1)=4*rands(1);    %初始位置  N*2的矩阵，第i行表示第i个粒子的位置
    Xm(i,2)=5*rand(1)+1;    %初始位置  N*2的矩阵，第i行表示第i个粒子的位置
    V(i,:)=rands(1,2);       %初始速度 
    
     if ((Xm(i,1))^2+(Xm(i,2))^2)>20
        fitness(i)=0; 
     else
        fitness(i)=Griewan(Xm(i,:));      
     end
     xm=Xm(i,:);            % xm为初始位置
     f(i,:)=fitness(i);   % f i   为粒子i初始适应度  
     %scatter3( xm(:,1),xm(:,2) ,fitness(i),'r*' );  %打印粒子初始位置
     %plot(xm,'ro')  
end
 
%%个体极值与群体极值
[bestfitness,bestindex]=max(fitness);   %bestfitness=最佳适应度，bestindex为最佳适应度的粒子，为第几个
pbest = Xm;               %个体最优初始为他本身位置  pbest存所有粒子的个人最优位置           
gbest = Xm(bestindex,:)  %全局最优为找出来的最佳适应值的那个粒子，第bestindex个粒子的位置
fitnesspbest=fitness ;    %个体历史最佳适应度为初始值
fitnessgbest=bestfitness;%全局历史最佳最佳适应度为找出来的最佳度

  这里的适应度函数，也就是我们需要得到解决的极值问题，此处以如下适应度函数为例：

f = 21.5+x1.*sin(4.*pi.*x1)+x2.*sin(20.*pi.*x2); %目标值 
((x1)^2+(x2)^2)<20                               %条件

3.迭代寻优部分

  迭代寻优是对整个种群的粒子进行更新速度和位置，以达到逐渐接近最优值的目的。每一次迭代都需要对所有粒子进行更新，而更新的主要对象是粒子的速度和位置。

  速度更新公式如下：

        

  其中d表示第d次迭代，i表是第i个粒子。v表示速度，所谓速度也就是下一步的飞行方向，该方向是一个向量，每个元素对应一个自变量，比如针对这个适应度函数由两个变量，则速度是一个1*2维的向量。其中pbest指的是种群最优极值对应的位置（也就是变量值），gbest表示自身粒子之前所遇到的最优位置，该公式表示下一步的飞行速度和方向由当前的飞行速度，此时的全局极值以及自身极值有关，其中w为当前速度的惯性，c1,c2对应各分量的权值，也表示对全局极值和自身极值的信任程度。

  更新速度后，接下里是对位置进行更新，更新这次迭代下该粒子的位置。

                                                    

  其中运动时间为1，所有运动步长就是对应的速度。

1，然后每一次迭代会更新每个粒子的速度和位置

2，更新位置后就可以计算每个粒子新的适应度，根据适应度计算每个粒子自身的最优极值，以及粒子群的全局极值，以供下一次迭代是计算速度和位置。

---

 三、matlab代码实现

首先是代码主体：pso.m

clc;
clear;
%%绘制目标函数曲线
figure
x = -4:0.01:4;
y = 1:0.01:6;
[X,Y] = meshgrid(x,y);
[row,col] = size(X);
for l = 1:col
	for h = 1:row
        z(h,l) = Griewan([X(h,l),Y(h,l)]);
	end
end
surf(X,Y,z);
title('PSP最优个体适应度','fontsize',12);
shading interp
hold on
 
%%参数初始化
c1 = 1.49445;
c2 = 1.49445;
N = 50; %　种群规模
D = 2;
T = 100; %　迭代次数
wa = 0.9;
wi = 0.4;
Xmax = 4 ;
Xmin = -4 ;
Ymax = 6 ;
Ymin = 1 ;
Vmax = 0.5 ; %　最大飞行速度
Vmin = -0.5 ; %　最小飞行速度
 
% Xm(i,:)粒子i位置
t1=cputime;
%%产生初始粒子和速度
for i=1:N
    Xm(i,1)=4*rands(1);    %初始位置  N*2的矩阵，第i行表示第i个粒子的位置
    Xm(i,2)=5*rand(1)+1;    %初始位置  N*2的矩阵，第i行表示第i个粒子的位置
    V(i,:)=rands(1,2);       %初始速度 
    
     if ((Xm(i,1))^2+(Xm(i,2))^2)>20
        fitness(i)=0; 
     else
        fitness(i)=Griewan(Xm(i,:));      
     end
     xm=Xm(i,:);            % xm为初始位置
     f(i,:)=fitness(i);   % f i   为粒子i初始适应度  
     %scatter3( xm(:,1),xm(:,2) ,fitness(i),'r*' );  %打印粒子初始位置
     %plot(xm,'ro')  
end
 
%%个体极值与群体极值
[bestfitness,bestindex]=max(fitness);   %bestfitness=最佳适应度，bestindex为最佳适应度的粒子，为第几个
pbest = Xm;               %个体最优初始为他本身位置  pbest存所有粒子的个人最优位置           
gbest = Xm(bestindex,:)  %全局最优为找出来的最佳适应值的那个粒子，第bestindex个粒子的位置
fitnesspbest=fitness ;    %个体历史最佳适应度为初始值
fitnessgbest=bestfitness;%全局历史最佳最佳适应度为找出来的最佳度
 
 
 
 
%%迭代寻优
for i=1:T    %从第一代到第T代，当前代数为i
    w=wa-(wa-wi)*i/T;     %wa=0.9 wi=0.4 w=0.9-0.5*(i/T)   w(0.9-0.4)前期w大，适合保持自己速度充分寻找，后期w小便于向其他学习
    for j=1:N  %N个粒子，当前为第j个
        %速度更新
        V(j,:)=w*V(j,:)+c1*rand*(pbest(j,:)-Xm(j,:))+c2*rand*(gbest-Xm(j,:));  %速度矩阵第j行更新
        V(j,find(V(j,:)>Vmax))=Vmax; 
        V(j,find(V(j,:)<Vmin))=Vmin;
        %种群更新
                            %位置矩阵第j行更新
        if ((Xm(j,1)+V(j,1))^2+(Xm(j,2)+V(j,2))^2)<20
            Xm(j,:)=Xm(j,:)+V(j,:);
            Xm(j,find(Xm(j,1)>Xmax))=Xmax;
            Xm(j,find(Xm(j,1)<Xmin))=Xmin;
            Xm(j,find(Xm(j,2)>Ymax))=Ymax;
            Xm(j,find(Xm(j,2)<Ymin))=Ymin;
            %适应度值更新
            fitness(j)=Griewan(Xm(j,:));                %第j个粒子的适应度更新
        end
 
    end
%      dt=plot3(Xm(:,1),Xm(:,2),fitness(:),'bo','linewidth',2);   % 动态绘图
%      pause(0.1)  
%      delete(dt)
 
    for j=1:N
        %个体最优更新
        if fitness(j)>fitnesspbest(j)   %第j个粒子更新后适应度
            pbest(j,:)=Xm(j,:);
            fitnesspbest(j)=fitness(j);
        end
        
        if fitness(j)>fitnessgbest
            gbest=Xm(j,:);
            fitnessgbest=fitness(j);
        end    
    end
    fitness;
    yy(i)=fitnessgbest;
end  
t2=cputime;
t=t2-t1; 
da=['种群: ',num2str(N),'；迭代: ',num2str(T),'；目标值: ',num2str(fitnessgbest),'；运行时间: ',num2str(t)];disp(da)
%%输出结果
 
plot3(gbest(1),gbest(2),fitnessgbest,'ro','linewidth',2)   % 打印最佳粒子的位置
figure
plot(yy)
title('PSP最优个体适应度','fontsize',12);
xlabel('进化代数','fontsize',12);
ylabel('适应度','fontsize',12);

 其次是适应度函数：Griewan.m

function f=Griewan(x)
x1=x(1);
x2=x(2);
f = 21.5+x1.*sin(4.*pi.*x1)+x2.*sin(20.*pi.*x2); %目标值 
end

最后是判断变量是否有效：panduan.m

function xx=panduan(x)
x1=x(1);
x2=x(2);
 if ((x1)^2+(x2)^2)<20
     xx=1;
 else
     xx=0;
 end   
end

总结

代码可以直接运行，如果有问题可以直接留言或者私信，后期会继续跟新粒子群算法的一些实际应用，比如背包问题。