PSO（粒子群）算法学习笔记

PSO（粒子群）算法学习笔记

PSO算法的概述

• PSO算法是一种全局优化的算法，它模拟的是鸟群或者鱼群的一种彼此共享信息去搜寻食物的过程。
• PSO算法与遗传算法(GA)类似，但是其少了GA中“交叉”，“变异”的操作，因此总体来看PSO算法操作更加简单，更加容易实现
• PSO算法可以理解为使用一群粒子来模拟鸟群个体的变化，每一个粒子在每一轮都会进行下一位置的搜索，下一轮粒子飞行的方向由当前时刻的方向以及自身搜索到的最优解方向和目前搜索到的全局最优解方向决定，粒子飞行的位置与当前位置与下一轮搜索的方向共同决定。

PSO算法最核心的两个公式如下：

V k n e x t = w V k n o w + C 1 R 1 ( P b e s t k − S k n o w ) + C 2 R 2 ( g b e s t − S k n o w ) V^{next}_k=wV_k^{now}+C_1R_1(P_{best}^{k}-S_k^{now})+C_2R_2(g_{best}-S^{now}_k) Vknext​=wVknow​+C1​R1​(Pbestk​−Sknow​)+C2​R2​(gbest​−Sknow​)

S k n e x t = S k n o w + V k n e x t S_k^{next}=S_k^{now}+V_k^{next} Sknext​=Sknow​+Vknext​

• 公式变量说明：
  V k n e x t V_k^{next} Vknext​代表了下一轮次搜索的方向，
  V k n o w V_k^{now} Vknow​代表了当先搜索的方向信息。
  S k n o w S^{now}_k Sknow​代表的是当前第K个粒子搜索的位置信息， S k n e x t S_k^{next} Sknext​代表的第K个下一次搜索到位置信息。 P b e s t k P_{best}^{k} Pbestk​代表的是当前第K个粒子搜索到的个体极值的位置信息
  g b e s t g_{best} gbest​代表了当前全局最优的位置信息。
  C 1 C_1 C1​是个体学习因子
  C 2 C_2 C2​是社会学习因子
  w w w是方向向量的惯性权重值。

• 公式原理说明：第K个对象拥有自身搜索的最优解信息 P b e s t k P^k_{best} Pbestk​,这是从其以往的经验中 得出来的信息；拥有目前为止的全局最优解 g b e s t g_{best} gbest​的信息，这是粒子之间共享的信息，为了从当前位置得到下一位置信息 S k n e x t S_k^{next} Sknext​，粒子必须知道心得方向信息 V k n e x t V_k^{next} Vknext​,此时PSO算法要求粒子必须在一方面维持当前方向 V k n o w V_k^{now} Vknow​的惯性，另一方面要考虑指向对象最优解的方向 P b e s t k − S k n o w P_{best}^{k}-S_k^{now} Pbestk​−Sknow​ ,以及指向全局最优解的方向 g b e s t − S k n o w g_{best}-S^{now}_k gbest​−Sknow​ 。这三者综合可以得到公式（1），而公式（2）描述的是粒子位置信息的改变方式，即每一个粒子在原先位置的基础上，加上新的方向信息，得到新的位置信息

PSO算法的步骤

1. 首先利用随机数初始化所有粒子的初始位置S与初始方向向量V，并初始化相关参数C1,C2,W等
2. 设置成本函数，并以成本函数最小（大）作为全局变量的搜索目标
3. 进行一轮搜索，记录每一个粒子历史搜索的最优值，最优位置
4. 结束一轮搜索后，比较当前轮次最小（大）成本与全局最小（大）成本，满足条件则替换
5. 根据公式（1）（2）迭代得到下一轮次的搜索方向与搜索位置
6. 重复执行（3），直达达到迭代上限或者得到满意的解

一个栗子

这里采用的是Branin函数（该函数常被用于二维的测试函数），Branin函数定义如下：

Branin函数的全局最优解为0.397887358，有三个全局最优解
这里我们使用PSO算法求解，代码如下：

clf;
PARAMETER=2;
v_max=5*ones(1,PARAMETER);   % 速度限制 上限
v_min=-5*ones(1,PARAMETER);  % 速度限制  下限
up_scope=[15 10];      % x y 坐标范围限制  上限
down_scope=[-15,-10];  % x y 坐标的范围限制  下限
C1=0.5;% 个体学习因子
C2=0.5;% 社会学习因子
w=0.6; %方向向量的惯性权重值
MAX_ITER=800;  % 最大迭代次数
PARTICLE_NUMBER=30;
%% 初始化位置 方向
for i=1:PARTICLE_NUMBER
    S(i,:)=(-down_scope+up_scope).*(rand(PARAMETER,1))'+down_scope;  % 位置
    V(i,:)=(v_max-v_min).*(rand(PARAMETER,1))'+v_min;   % 方向
end
%对象最优解
cost_personal=(1e10.*ones(PARTICLE_NUMBER,1));
pos_personal=S;
%全局最优解
cost_global=100;
pos_global=[];
% 最终结果
TRUE=0.397887358;
err=1000;
iter=0;
while (iter<=MAX_ITER && err>1e-7)
    iter=iter+1; % 迭代加一
    %  一轮中各个个体对应的代价最小值  位置值
    for i=1:PARTICLE_NUMBER
        cost(i,1)=TestBR(S(i,:));
        if cost(i,1)<cost_personal(i,1)
            cost_personal(i,1)=cost(i,1);   % 个体代价更新
            pos_personal(i,:)=S(i,:);  % 位置信息更新
        end
    end
    [val,position]=min(cost_personal);
    if (val<cost_global)
        cost_global=val;
        pos_global=pos_personal(position,:);
    end
    V=w*V+C1*rand*(pos_personal(:,:)-S(:,:))+C2*rand*(ones(30,1)*(pos_global)-S(:,:));
    S=S+V;
    err=abs(cost_global-TRUE);
    fee(iter)=val;

    figure(1);
    plot(pos_personal(:,1),pos_personal(:,2),".r",'Markersize',15);
    axis([down_scope(1) up_scope(1) down_scope(2) up_scope(2)]);
    set(get(gca, 'XLabel'), 'String', 'x坐标');
    set(get(gca, 'YLabel'), 'String', 'Y坐标');
    set(get(gca, 'Title'), 'String', 'PSO算法寻优');
    pause(0.1);
end
disp("min");
disp(cost_global);
disp ("min location");
disp(pos_global);
figure(2);
plot(1:length(fee),fee,'b');
set(get(gca, 'XLabel'), 'String', '迭代次数');
set(get(gca, 'YLabel'), 'String', '当前轮次搜索最优值');
hold on ;

function y=TestBR(x)  %的时刻
    x1=x(1);   x2=x(2);
    y=(x2-5.1/(4*pi^2)*x1^2 +5/pi*x1-6)^2+10*(1-1/(8*pi))*cos(x1)+10;
end

刚开始时刻，30个点随机分布：

最终结果基本汇聚在一起：


迭代次数：