粒子群算法PSO

粒子群算法PSO

简介

粒子群优化算法(PSO)是一种进化计算技术(evolutionary computation)，1995 年由Eberhart 博士和kennedy 博士提出，源于对鸟群捕食的行为研究 。该算法最初是受到飞鸟集群活动的规律性启发，进而利用群体智能建立的一个简化模型。粒子群算法在对动物集群活动行为观察基础上，利用群体中的个体对信息的共享使整个群体的运动在问题求解空间中产生从无序到有序的演化过程，从而获得最优解。

基本思想

粒子群算法是模拟群体智能所建立起来的一种优化算法，粒子群算法可以用鸟类在一个空间内随机觅食为例，所有的鸟都不知道食物具体在哪里，但是他们知道大概距离多远，最简单有效的方法就是搜寻目前离食物最近的鸟的周围区域。

对粒子群优化算法操作的一个简单解释如下：每一个粒子都代表了当前优化任务的一个可能解决方案。在每次迭代过程中，每个粒子都会朝着自己的最优解的方向加速，也会朝着种群中任何粒子迄今为止发现的全局最佳位置的方向加速。这意味着，如果一个粒子发现了一个更好的解，所有其他粒子都会靠近它，在这个过程中不断地搜索最优解。可以总结出粒子群算法地三条简单规则：（1）飞离最近的个体，以避免碰撞；（2）飞向目标；（3）飞向群体的中心。

假设存在一个维度为S的搜索空间，由$m$ 个粒子组成粒子种群，其中第$i$ 个粒子用一个$S$ 维的向量表示，具体为$Xi=(x_{i1},x_{i2},\dots ,x_{iS})$。将$X_i$代入目标函数就可以算出其适应度值，适应度值的大小对应着粒子位置即可行解的好坏。另外，第$i$ 个粒子移动的速度是$S$ 维向量，记为$Vi=(v_{i1},v_{i2},\dots ,v_{iS})$。第$i$个粒子搜索到的最优位置为$P_i=(p_{i1},p_{i2},\dots ,p_{iS})$，整个粒子群搜索到的最优位置为$P_g=(p_{g1},p_{g2},\dots ,p_{gS})$。

粒子位置和速度的更新公式如下：
KaTeX parse error: No such environment: eqnarray at position 8: \begin{̲e̲q̲n̲a̲r̲r̲a̲y̲}̲ V_i&=&V_i+c_1*…
$i＝1，2，\dots ，m$，$m$是该群体中粒子的总数；$V_i$ 是粒子的速度；$rand()$是介于$(0,1)$之间的随机数； $X_i$ 是粒子的当前位置。 $c_1$和$c_2$是学习因子， 在每一维，粒子都有一个最大限制速度$V_{max}$，如果 某一维的速度超过设定的$V_{max}$ ，那么这一维的速度 就被限定为$V_{max}$ 。以上面两个公式为基础，形成了后来PSO 的标准形式。

一般在实际应用中，会对速度更新公式进行修正，即加入惯性权重$w$：
$$V_i=w\ast V_i+c_1\ast rand()\ast (P_i-X_i)+c_2\ast rand()\ast (P_g-X_i)\text{\hspace{1em}(3)}$$

算法流程

标准PSO算法的流程：

Step1:初始化一群微粒(群体规模为m)，包括随机位置和速度；

Step2:评价每个微粒的适应度；

Step3:对每个微粒，将其适应值与其经过的最好位置$P_i$作比较，如果较好，则将其作为该粒子当前的最好位置$P_i$;

Step4:对每个微粒，将其适应值与全局最好位置$P_g$作比较，如果较好，则将其作为当前全局最好位置$P_g$;

Step5:根据(2)、(3)式调整微粒速度和位置；

Step6：未达到结束条件则转Step2。

迭代终止条件根据具体问题一般选为最大迭代次数$k$或(和)微粒群迄今为止搜索到的最优位置满足预定最小适应阈值$ϵ$。

算法实现（MATLAB）

一般的PSO实现¹：

function [xm,fv,pbest]=stdPSO(fun,sizepop,c1,c2,w,maxg,D)
% Standard PSO
% Input
%fitness: 待优化的目标函数
%sizepop: 粒子数目
%c1: 学习因子1
%c2: 学习因子2
%w: 惯性权重
%maxg: 最大迭代次数
%D: 自变量的维度
%Output 
%xm: 目标函数取最小值时的自变量值
%fv: 目标函数的最小值
%% 参数初始化
% 初始速度和种群上下边界值
Vmax=1;
Vmin=-1;
popmax=5;
popmin=-5;
%% 产生初始粒子和速度
for i=1:sizepop
    %随机产生一个种群
    pop(i,:)=5*rands(1,D); %初始种群
    V(i,:)=rands(1,D); %初始化速度
    %计算适应度
    fitness(i)=fun(pop(i,:)); %染色体的适应度
end
%找最好的染色体
[bestfitness bestindex]=min(fitness);
zbest=pop(bestindex,:); %全局最佳
gbest=pop; %个体最佳
fitnessgbest=fitness; %个体最佳适应度值
fitnesszbest=bestfitness; %全局最佳适应度值

%% 迭代寻优
for i=1:maxg   %迭代次数
    for j=1:sizepop
        %速度更新
        V(j,:) = w*V(j,:) + c1*rand*(gbest(j,:) - pop(j,:)) + c2*rand*(zbest- pop(j,:));
        V(j,find(V(j,:)>Vmax))=Vmax;
        V(j,find(V(j,:)<Vmin))=Vmin;
        %种群更新
        pop(j,:)=pop(j,:)+V(j,:);
        pop(j,find(pop(j,:)>popmax))=popmax;
        pop(j,find(pop(j,:)<popmin))=popmin;
        %自适应变异
        if rand>0.8
            k=ceil(2*rand);
            pop(j,k)=rand;
        end
        %适应度值
        fitness(j)=fun(pop(j,:));
        %个体最优更新
        if fitness(j) < fitnessgbest(j)
            gbest(j,:) = pop(j,:);
            fitnessgbest(j) = fitness(j);
        end
        %群体最优更新
        if fitness(j) < fitnesszbest
            zbest = pop(j,:);
            fitnesszbest = fitness(j);
        end
    end
    yy(i)=fitnesszbest;
end
% 结果输出
xm=zbest; %最佳个体值
fv=fitnesszbest;
pbest=yy;

测试函数Rastrigin
$$f(x)=A\ast n+\sum _{i=1}^n[x_i^2-A\ast \cos (2\pi x_i)]$$

function F= Rastrigin(x)
%Rastrigin function
% f(x)=An+\sum_{i=1}^n[x_i^2-Acos(2pix_i)]
% when A=10 and x_i\in[-5.12,5.12]. It has a global minimum at x=0 where f(x)=0

[m,n]=size(x);
A=10;
F=zeros(m,1);
for i=1:m
    F(i)=A*n+sum(x(i,:).^2-A*cos(2*pi*x(i,:)));
end

绘出测试图形：

% %% 经典测试函数
function DrawRastrigin()
% % finding the mimimum value.
% clc %清屏
% clear all; %删除workplace 变量
% close all; %关掉显示图形窗口
% 绘制Rastrigin 函数图形
x = [-5 : 0.05 : 5 ];
y = x;
[X,Y] = meshgrid(x,y);
[row,col] = size(X);
for l = 1 :col
for h = 1 :row
z(h,l) = Rastrigin([X(h,l),Y(h,l)]);
end
end
mesh(X,Y,z);
shading interp
end

利用PSO进行求解：

clear
clc
sizepop=200;
c1=1.49445;
c2=1.49445;
w=0.5;
maxg=500;
D=2;
t1=clock;
[xm,fv,pbest]=stdPSO(@Rastrigin,sizepop,c1,c2,w,maxg,D);
t2=clock;
disp(['PSO时间：',num2str(etime(t2,t1)),'s']);
plot(pbest,'x-');
title(['时间为',num2str(etime(t2,t1)),'s','; 最佳值为：',num2str(fv)])

可看出，该函数的寻优是收敛的，最优个体值接近于理论值(0,0)；但是当维度增加时，标准PSO就不是这么有效了。比如维度增加到10维，最终结果为：

结果离最优值还有差距，这就需要对标准PSO进行改进了。

PSO算法的向量化实现

在进行MATLAB计算时，应减少循环的出现，尽量用向量或矩阵计算代替循环以此来优化程序运行速度。下面，我将用向量化的方式来重写标准PSO算法程序。

function [xm,fv,pbest]=stdPSO_v1(fun,sizepop,c1,c2,w,maxg,D)
% Standard PSO
% Input
%fitness: 待优化的目标函数
%sizepop: 粒子数目
%c1: 学习因子1
%c2: 学习因子2
%w: 惯性权重
%maxg: 最大迭代次数
%D: 自变量的维度
%Output 
%xm: 目标函数取最小值时的自变量值
%fv: 目标函数的最小值

%% 参数初始化
% 初始速度和种群上下边界值
Vmax=1;
Vmin=-1;
popmax=5;
popmin=-5;
%% 产生初始粒子和速度
% for i=1:sizepop
%     %随机产生一个种群
%     pop(i,:)=5*rands(1,D); %初始种群
%     V(i,:)=rands(1,D); %初始化速度
%     %计算适应度
%     fitness(i)=fun(pop(i,:)); %染色体的适应度
% end
pop=5*rands(sizepop,D);
V=rands(sizepop,D);
fitness=fun(pop);

%找最好的染色体
[bestfitness,bestindex]=min(fitness);
zbest=pop(bestindex,:); %全局最佳
gbest=pop; %个体最佳
fitnessgbest=fitness; %个体最佳适应度值
fitnesszbest=bestfitness; %全局最佳适应度值
pbest=zeros(1,maxg);

%% 迭代寻优
for i=1:maxg   %迭代次数
    %速度更新
    V=w*V+c1*rand(sizepop,1).*(gbest - pop)+c2*rand(sizepop,1).*(zbest- pop);
    V(V>Vmax)=Vmax;
    V(V<Vmin)=Vmin;
    %种群更新
    pop=pop+V;
    pop(pop>popmax)=popmax;
    pop(pop<popmin)=popmin;
    %自适应变异
    tmp=rand(sizepop,D)>0.8;
    tmpr=rand(sizepop,D);
    pop(tmp)=tmpr(tmp);
    %适应度值
    fitness=fun(pop);
    %个体最优更新
    tmp=fitness<fitnessgbest;
    if any(tmp) == true
        gbest(tmp,:)=pop(tmp,:);
        fitnessgbest(tmp)=fitness(tmp);
    end
    %群体最优更新
    tmp=fitness<fitnesszbest;
    if any(tmp) == true
        tmpz=pop(tmp,:);
        tmpf=fitness(tmp);
        [fitnesszbest,tmpi]=min(tmpf);
        zbest=tmpz(tmpi,:); %全局最佳
    end
    pbest(i)=fitnesszbest;
end
% 结果输出
xm=zbest; %最佳个体值
fv=fitnesszbest;

重新对上面的最优化问题求解，得到：

时间对比为0.902s比0.084s，向量化的程序运行时间只有原始程序的1/10。

reference

---

1. MATLAB 优化算法案例分析与应用 ↩︎