1.PSO概述

1.1基本思想

$\bullet$ 群体迭代，粒子在解空间(n维)追随最优的粒子进行搜索。

$\bullet$ 粒子群算法的思想源于对鸟群捕食行为的研究。

$\bullet$模拟鸟集群飞行觅食的行为，鸟之间通过集体的协作使群体达到最优目的。

PSO基础：信息的社会共享

---

1.2特点

$\bullet$ 简单易行

$\bullet$ 收敛速度快

$\bullet$ 设置参数少

---

1.3算法介绍

主要思路：

$\bullet$ 每个寻优问题解都被想象成一只鸟，称为“粒子” 。所有粒子都在一个D维空间进行搜索。

$\bullet$ 根据适应度函数确定每个粒子的适应值以判断当前位置好坏。

$\bullet$ 每个粒子带有记忆功能，可以记录搜索到的最佳位置。

$\bullet$ 每个粒子都有一个速度决定飞行的距离和方向。这个速度根据它本身的飞行经验以及同伴的飞行经验进行动态调整。

变量：在D维空间中，有N个粒子

变量名称	含义
xi=(xi1,xi2,…xiD)	粒子i的位置，表示D维空间的空间坐标
vi=(vi1,vi2,…viD)	粒子i的速度，分解在D维空间里的D个分速度
pbest1=(pi1,pi2,…piD)	粒子i个体经过的最好位置
gbest=(g1,g2,…gD)	种群经历过的最好位置

通常在每一维都限定了位置的变化范围以及速度的变化范围。

---

1.4算法流程图

---

1.5构成要素

$\bullet$ 群体大小m：

m是一个整型参数，当m很小时，陷入局部最优的可能性很大；当m很大时，PSO的优化能力好，但是当种群数目增长至一定水平时，再增长不再有显著的作用。

$\bullet$ 权重因子：参与速度更新公式如下（其中$v_{id}^k$表示第k次迭代粒子i在第d维空间中的速度，$r_1、r_2$表示两个随机函数，取值范围为[0,1]，增加随机搜索范围）
$$v_{id}^k=\omega v_{id}^{k-1}+c_1{r}_1(pbes{t}_{id}-x_{id}^{k-1})+c_2{r}_2(gbes{t}_{id}-x_{id}^{k-1})$$

$\bullet$ 惯性因子：ω
非负数，调节对解空间的范围搜索。当ω=0时，失去对粒子本身的速度的记忆。

---

$\bullet$ 学习因子：c₁、c₂
（1）$c_1{r}_1(pbes{t}_{id}-x_{id}^{k-1})$表示自我认知部分，其中$c_1$决定了自身经验对粒子速度的影响程度，它保证了粒子向自己的历史最优位置靠近。当$c_1=0$时，为无私型粒子群算法，“只有社会，没有自我”，容易丧失种群多样性，易陷入局部最优。
（2）$c_2{r}_2(gbes{t}_{id}-x_{id}^{k-1})$表示社会经验部分，其中$c_2$决定了群体经验对粒子个体速度的影响程度，它保证了粒子能够向群体中的其他粒子学习，使粒子在飞行时向领域内所有粒子找到的历史最优位置靠近。当$c_2=0$时，为自我认知型粒子算法，“只有自我，没有社会”，完全没有信息的社会共享，导致算法收敛速度缓慢。
（3）当$c_1、c_2$都不为0时，称为完全型粒子群算法，完全型粒子群算法更容易保持收敛速度和搜索效果的均衡，是较好的选择。

$\bullet$ 最大速度$V_m$：

用来维护算法的探索能力和开发能力的平衡。当$V_m$ 较大时，探索能力增强，但粒子容易飞过最优解；当$V_m$ 较小时，开发能力增强，但粒子容易陷入局部最优；$V_m$ 一般设为每维变量变化范围的10%~20%。

$\bullet$ 领域的拓扑结构：

包括两种结构，一种是将种群内$所有个体$都作为粒子的领域，另一种是只将群体中的$部分个体$作为粒子的领域。领域拓扑结构决定了群体历史最优位置。由此，将粒子群算法分为：全局粒子群算法和局部粒子群算法。
$\bullet$ $全局粒子群算法$算法收敛速度快，但是容易陷入局部最优。
$\bullet$ $局部粒子群算法$算法收敛速度慢，但是狠点陷入局部最优。

$\bullet$ 停止准则：一般有两种情况，一为达到最大迭代步数，二为出现可接受的满意解。

$\bullet$ 粒子空间的初始化：较好地选择粒子初始化空间，将大大缩短收敛时间。

在这些构成要素中，权重因子和领域定义影响较大。

---

2.MATLAB实例

求解如下高斯函数的最优值：
$$img=\frac{1}{2\pi {\sigma }^2}{e}^{-\frac{x^2+y^2}{2{\sigma }^2}}$$

一共由三个函数组成：

函数名称	作用
main.m	主函数，调用其它函数
update_par.m	更新每个粒子的信息
compute_fit.m	求粒子适应度

---

$\bullet$ main.m：

clear all;close all;clc;

%以二维空间为例
[x y]=meshgrid(-100:100,-100:100);
sigma=50;
img = (1/(2*pi*sigma^2))*exp(-(x.^2+y.^2)/(2*sigma^2)); %目标函数，高斯函数
mesh(img);
hold on;
n=10;   %粒子群粒子个数

%初始化粒子群，定义结构体
%结构体中八个元素，分别是粒子坐标，粒子速度，粒子适应度，粒子最佳适应度，粒子最佳坐标
par=struct([]);
for i=1:n
    par(i).x=-100+200*rand();   %[-100 100]对x位置随机初始化
    par(i).y=-100+200*rand();   %[-100 100]对y位置随机初始化
    par(i).vx=-1+2*rand();      %[-1 1]对vx速度随机初始化
    par(i).vy=-1+2*rand();      %[-1 1]对vy速度随机初始化
    par(i).fit=0;               %粒子适应度为0初始化
    par(i).bestfit=0;           %粒子最佳适应度为0初始化
    par(i).bestx=par(i).x;      %粒子x最佳位置初始化
    par(i).besty=par(i).y;      %粒子y最佳位置初始化
end
par_best=par(1);    %初始化粒子群中最佳粒子

%设置迭代次数最大为1000次
for k=1:1000    
    plot3(par_best.x+100,par_best.y+100,par_best.fit,'g*'); %画出最佳粒子的位置，+100为相对偏移
    for p=1:n
        [par(p) par_best]=update_par(par(p),par_best);  %更新每个粒子信息         
    end  
end

---

$\bullet$ update_par.m：

function [par par_best]=update_par(par,par_best)
    
    %Px=Px+Pv*t,这里t=1,Px为当前粒子的位置，Pv为当前粒子的速度
    par.x=par.x+par.vx;   
    par.y=par.x+par.vy;   
    
    par.fit=compute_fit(par);    %计算当前粒子适应度
    
    %Pv=Pv+(c1*rand*(Gx-Px))+(c2*rand*(PBx-Px))
    %这里c1,c2为加速因子
    %Gx为具有最佳适应度粒子的位置
    %PBx为当前粒子的最佳位置
    c1=1;
    c2=1;
    %速度更新公式，上面已经讲过
    par.vx=par.vx+c1*rand()*(par_best.x-par.x)+c2*rand()*(par.bestx-par.x);   
    par.vy=par.vy+c1*rand()*(par_best.y-par.y)+c2*rand()*(par.besty-par.y);
 
    if par.fit>par.bestfit      %如果当前粒子适应度要好于当前粒子最佳适应度
        par.bestfit=par.fit;    %则更新当前粒子最佳适应度
        par.bestx=par.x;        %更新当前粒子最佳位置
        par.besty=par.y;
        if par.bestfit>par_best.fit     %如果当前粒子最佳适应度好于种群最佳粒子适应度
            par_best.fit=par.bestfit;   %则更新最佳粒子适应度
            par_best.x=par.x;           %更新最佳粒子位置
            par_best.y=par.y;
        end
    end

end

---

$\bullet$ compute_fit.m：

function re=compute_fit(par)
    x=par.x;
    y=par.y;
    sigma=50;
    if x<-100 || x>100 || y<-100 || y>100
        re=0;        %超出范围适应度为0
    else            %否则适应度按目标函数求解
        re= (1/(2*pi*sigma^2))*exp(-(x.^2+y.^2)/(2*sigma^2)); 
    end
end

---

结果：