粒子群算法

粒子群算法及Matlab实现

一.粒子群算法原理介绍

定义

粒子群优化算法(Particle Swarm optimization,PSO)又翻译为粒子群算法、微粒群算法、或微粒群优化算法。是通过模拟鸟群觅食行为而发展起来的一种基于群体协作的随机搜索算法。粒子群优化算法的基本思想是通过群体中个体之间的协作和信息共享来寻找最优解．
举一个简单的例子，假设有一个鸟群在一个山区觅食，只有这片地区的最高点才有食物。假设鸟群在一开始随机分布，所有个体之间的信息是可以互相交流。鸟群之间的个体相互交流，发现鸟A的位置是所有鸟之中最高的，于是所有的鸟都会向鸟A的方向聚集过去，直到在这个过程中鸟B突然发现自己经过的一个地方好像比鸟A的位置还要高，于是所有个体又会向鸟B的方向聚集过去，直到到达最高点。
将山峰高度表示为$z=f(x,y)$，$x,y$为位置坐标，于是这个问题就抽象成了一个求目标函数最大值的问题，其中$z=f(x,y)$为目标函数，$x,y$为决策变量。
在这个过程中，我们需要知道信息有每个个体的位置信息（即决策变量的值），速度（决定了个体向哪个方向移动）以及最优值，其中包括每个个体的局部最优值（即飞行过程中每个鸟自己发现的最高点）和所有个体的全局最优值（即飞行过程中整个鸟群所有个体发现的最高点）。

算法流程

首先随机初始化一个粒子群$B=(B_1,B_2,...B_n)$，其中每个粒子的位置信息表示为$X_i=[x_{i1},x_{i2},...x_{id}]$，速度为$V_i=[v_{i1},v_{i2},...v_{id}]$，然后粒子通过不断迭代来寻找最优解，第t+1次迭代中鸟群的位置和速度可以更新为：
$$\begin{array}{rl}& X_i^{t+1}=X_i^t+V_i^{t+1}\\ & V_i^{t+1}=\omega V_i^t+c_1{r}_1(P_{ibest}^t-X_i^t)+c_2{r}_2(G_{best}^t-X_i^t)\end{array}$$
$$\begin{gathered} \omega 为惯性系数，是一个常数，\omega 较大时，全局寻优能力强，局部寻优能力弱，\omega 较小时则相反 \\ {c}_i(i=1,2)都是常数，分别称为个体学习因子和社会学习因子，一般取c_1=c_2\in [0,4] \\ {r}_i(i=1,2)是位于(0，1)区间的随机数 \end{gathered}$$
$P_{ibest}为局部最优值，G_{best}为全局最优值$
上式中$X_i^{t+1}$由三个部分组成，第一个部分$\omega V_i^t$为惯性项，表示个体沿原先方向运动具有的惯性，第二项$c_1{r}_1(P_{ibest}^t-X_i^t)$叫做个体自信度，表示粒子受到自身历史最好位置的吸引力，即粒子向自身最好历史位置移动的程度，第三项$c_2{r}_2(G_{best}^t-X_i^t)$为全局自信度，表示粒子受到全局最优位置的吸引程度，即向全局最优位置移动的意愿。
如此不断迭代，直到到达最大迭代次数或找到最优解时停止。

Matlab代码实现

%此程序目的是求解目标函数最小值
clc
clear

c1=2;%每个粒子的个体学习因子，加速度常数
c2=2;%每个粒子的社会学习因子，加速度常数
w=0.8;%惯性因子

maxgen=300;%设置最大迭代次数
sizepop=50;%设置种群粒子数目
nVar=2;%决策变量个数
vmax=5;
vmin=-5;%设置飞行速度范围，防止每次迭代时步长过大错过最优解甚至出现不收敛的情况
spacemax=100;
spacemin=-100;%设置粒子分布范围，即最大搜索范围

%初始化种群
for i=1:sizepop
    population(i,:)=(spacemax-spacemin)*rand(1,nVar)+spacemin;%种群位置初始化
    speed(i,:)=(vmax-vmin)*rand(1,nVar)+vmin;%速度初始化，每个方向（决策变量）都有一个速度
    goal(i)=func(population(i,:));%计算每个个体的目标函数值
end

[bestgoal,bestindex]=max(goal);
gbest=population(bestindex,:);%第一代全局最优位置
pbest=population;%第一代的局部最优位置就是种群当前位置
pbestgoal=goal;%局部最优值
gbestgoal=bestgoal;%全局最优值

for i=1:maxgen
    for j=1:sizepop
        %w=0.99^i*w;%w不断迭代
        speed(j,:)=w.*speed(j,:)+c1*rand*(pbest(j,:)-population(j,:))+c2*rand*(gbest-population(j,:));%速度更新
        speed(j,find(speed(j,:)>vmax))=vmax;
        speed(j,find(speed(j,:)<vmin))=vmin;

        population(j,:)=population(j,:)+speed(j,:);
        speed(j,find(population(j,:)>spacemax))=-speed(j,find(population(j,:)>spacemax));%将碰壁的粒子方向反转
        speed(j,find(population(j,:)<spacemin))=-speed(j,find(population(j,:)<spacemin));
        population(j,find(population(j,:)>spacemax))=spacemax;
        population(j,find(population(j,:)<spacemin))=spacemin;

        goal(j)=func(population(j,:));
    end

    for j=1:sizepop
        %更改局部最优值和局部最优位置
        if goal(j)<pbestgoal(j)
            pbestgoal(j)=goal(j);
            pbest(j,:)=population(j,:);
        end

        %更新全局最优值和全局最优位置
        if goal(j)<gbestgoal
            gbest=population(j,:);
            gbestgoal=goal(j);
        end
    end
    yy(i)=gbestgoal;
end

plot(yy);
title('迭代过程中全局最优值变化过程');
xlabel('迭代次数');
ylabel('最优值');
fprintf('最优值%f\n',gbestgoal);
fprintf('决策变量%f',gbest);

目标函数：

function f = func(x)
f=sum(x.^2)+x(1)^2/exp(x(2));
end

求解结果为：

最优值0.000000
决策变量0.000000决策变量-0.000000

我们使用matlab工具箱求解，可以得到相似的结果：

[x,fval]=fmincon('func',rand(2,1))

结果：
x =

   1.0e-07 *

   -0.1094
   -0.0794


fval =

   3.0235e-16

约束的处理

上述的算法只对变量范围进行了约束，但是在实际问题中问题还常常具有等式约束和不等式约束，数学模型如下所示：
目标函数：$\min (f(\vec{X}))$
决策变量: $\overrightarrow{X}=[x_1,x_2,...,x_n]$
约束条件：$\{\begin{array}{l}a\le x_i\le b\\ g(\vec{X})\le c\\ h(\vec{X})=d\end{array}$

对于具有约束的最优化问题，可以通过惩罚函数将其转化为无约束的目标函数，及对于不满足约束条件的个体，我们可以将其的目标函数值（适应度）加上一个很大的数，从而排除掉这个个体所在位置（对于最小化问题而言，当我们给某个个体目标值加上一个很大的数后，就阻止了其他个体向该方向的移动，相当于排除了这片区域）。
例如对于约束条件：$x_1+x_2>1$我们只需要将func.m函数修改成如下所示：

function f = func(x)
f=sum(x.^2)+x(1)^2/exp(x(2));

if x(1)+x(2)<1
    f=f+10000;
end

结果为：

最优值0.606476
决策变量0.378878决策变量0.621122

使用matlab工具箱：

[x,fval]=fmincon('func',rand(2,1),-[1,1],-1)

解得：
x =

    0.3789
    0.6211


fval =

    0.6065

粒子群算法变异策略

为了避免粒子群算法陷入局部最优，我们可以引入变异操作，即当我们进行了二十次或者三十次迭代之后，我们将部分个体的位置重新打散，使其均匀分布到整个搜索空间中，读者可以自行尝试。