粒子群算法（PSO算法）

粒子群算法概述

1.粒子群优化算法（Particle Swarm Optimization，简称PSO）。粒子群优化算法是在1995年由Kennedy博士和Eberhart博士一起提出的，它源于对鸟群捕食行为的研究。

2.基本核心是利用群体中的个体对信息的共享从而使得整个群体的运动在问题求解空间中产生从无序倒有序的演化过程，从而获得问题的最优解。

3.粒子群算法模拟鸟群的捕食过程，将待优化的问题看做是捕食的鸟群，解空间看做是鸟群的飞行空间，空间的每只鸟的位置即是粒子群算法在解空间的一个粒子，也就是待优化问题的一个解。

粒子群算法满足以下假设：

1.粒子被假定没有体积没有质量，本身的属性只有速度和位置。

2.每个粒子在解空间中运动，它通过速度改变其方向和位置。

3.通常粒子将追踪当前的最优粒子以经过最少代数的搜索达到最优。

PSO算法模型

1.PSO算法首先在可行解空间中初始化一群粒子，每个粒子都代表极值最优化问题的一个潜在最优解，用位置、速度和适应度值三项指标表示该粒子的特征。

2.粒子在解空间中运动，通过跟踪个体极值Pbest和群体极值Gbest来更新个体位置。

3.粒子每更新一次位置和速度，就计算一次适应度值，并且通过比较新粒子的适应度值和个体极值、群体极值的适应度更新个体极值Pbest和群体极值Gbest的位置。

假设在一个n维的搜索空间，有m个粒子组成一个粒子群 ${X}=\left ( X{_{1}},X_{2},...,X_{m} \right )$ ，其中：

注：

掌握PSO算法需要从以下两个角度进行理解：

1.粒子在一定程度上离开原先的搜索轨迹，向新的方向进行搜索，体现了一种向未知区域开拓的能力，类似于全局搜索。

2.粒子在一定程度上继续在原先的搜索轨迹上更细一步的搜索，主要针对搜索过程中所搜索到的区域进行更进一步的搜索，类似于局部搜索。

基本粒子群算法

在找到这两个最优解时，粒子根据下面的式子来更新自己的速度和位置：

注：

Vij 指第i个粒子在第j个分量的速度，k代表第k次循环。r1和r2是介于（0,1）的随机数。c1和c2是学习因子。

关于学习因子c1和c2：

1.如果c1=c2=0，则说明粒子将以当前的飞行速度飞到边界。此时，粒子仅能搜索有限的区域，所以很难找到最优解。

2.如果c1=0，则为“社会”模型，粒子只有群体经验而缺乏认知能力。此时，算法收敛速度快，但容易陷入局部最优。

3.如果c2=0，则为“认知”模型，粒子没有社会共享信息，粒子之间没有信息的交互。此时，算法找到最优解的概率很小。

易知该公式由三部分组成：

1.第一部分为“惯性(inertia)”或“动量(momentum)”部分，反映了粒子的运动“习惯(habit)”，代表粒子有维持自己先前速度的趋势。

2.第二部分为“认知(cognition)”，反映了粒子对自身历史经验的记忆(memory)或回忆(remembrance)，代表粒子有向自身历史最佳位置逼近的趋势。

3.第三部分为“社会(social)”部分，反映了粒子间协同合作与知识共享的群体历史经验，代表粒子有向群体或者邻域历史最佳位置逼近的趋势。

标准粒子群算法

Shi和Eberhart于1998年在IEEE国际计算学术会议上发表了题为“A Modified Particle Swarm Optimizer”的论文，首次在速度更新公式中加入了惯性权重，如下式所示：

其中w>0叫做惯性因子。w越大，全局寻优能力越强，局部寻优能力越弱；w越小，全局寻优能力越弱，局部寻优能力越强。此时的粒子群算法被称作是标准粒子群算法。

经验：

实验表明，动态的w能获得更好的寻优结果。在搜索过程中可以对w进行动态调整：可以初始给w赋子较大正值，随着搜索的进行，可以线性地使逐渐减小，这样可以保证在算法开始时，各粒子能够以较大的速度步长在全局范围内探测到较好的区域：而在搜索后期，较小的值则保证粒子能够在极值点周围做精细的搜索，从而使算法有较大的概率向全局最优解位置收敛。对W进行调整，可以权衡全局搜索和局部搜索能力。

目前，采用较多的动态惯性权重值是线性递减权值策略，其表达式如下：

其中， $T_{max}$ 表示最大迭代次数； $w_{max}$ 和 $w_{min}$ 分别表示最大和最小惯性权重。在大多数应用中，通常有 $w_{max}=0.9$ , $w_{min}=0.4$ 。

当然，除了线性递减权值策略，还有线性微分递减策略：

压缩粒子群算法

当 $w=\lambda$ 时，标准粒子算法就会退化为压缩粒子群算法：

此时的粒子群算法被称作是压缩粒子群算法。其中 $\lambda$ 为压缩因子：

离散粒子群算法

基本粒子群算法是在连续域中搜索函数极值的有利工具。为了处理离散空间的问题，Kennedy和Eberhart又提出了一种离散二进制版的粒子群算法。在此算法中，将离散空间问题映射到连续粒子空间，并适当修改粒子群算法来求解。在计算上仍保留经典粒子群算法中的速度和位置更新规则。

在离散粒子群算法中，将离散问题空间映射到连续粒子运动空间，并做适当的修改。仍然保留经典粒子群算法中的速度和位置更新规则。粒子在状态空间的取值只限于0，1两个值，而速度的每一个位代表的是粒子位置所对应的位取值为0/1的可能性。因此在离散粒子群算法中，粒子速度的更新公式依然保持不变，但是个体最优位置和全局最优位置每一位的取值只能为0/1。

更新方式：

，

其中，r ~ U(0 , 1) , $V_{ij}$ 表示位置 $X_{ij}$ 取值为1的可能性，其更新公式与基本粒子群算法一致。

粒子群算法流程

Step 1.种群初始化：群体规模 $N$ ，每个粒子的位置 $X$ 和速度 $V_{i}$ 。

Step 2.计算每个粒子的适应值 $fit_{i}$ ;。

Step 3.对每个粒子，用它的适应度值 $fit_{i}$ 和个体极值 $Pbest_{i}$ 比较。如果 $fit_{i}<Pbest_{i}$ ，则用 $fit_{i}$ 替换 $Pbest_{i}$ 。

Step 4.对每个粒子，用它的适应度值 $fit_{i}$ 和全局极值 $Gbest_{i}$ 比较。如果 $fit_{i}<Gbest_{i}$ ，则用 $fit_{i}$ 替换 $Gbest_{i}$ 。

Step 5.迭代更新粒子的速度 $V_{i}$ 和位置 $X$ 。

Step6.进行边界条件处理。

Step 7.判断算法终止的条件是否满足。若满足，则结束算法并输出优化结果；否则，返回Step2。

算法的参数说明

注：

边界条件处理：

当某一维或若干维的位置或速度超过设定值时，采用边界条件处理策略可将粒子的位置限制在可行搜索空间内，这样能避免种群的膨胀与发散，也能避免粒子大范围地盲目搜索，从而提高了搜素效率。具体的方法有很多种，比如通过设置最大位置限制 $x_{max}$ 和最大速度限制 $v_{max}$ ，当超过最大位置或最大速度时，在范围内随机产生一个数值代替，或者将其设置为最大值，即边界吸收。

例题

适应度函数fit.m

function y = fit(x)
% 函数用于计算粒子的适应度
% x:输入粒子
% y:粒子适应度值
y = x(1).^2 + x(2).^2 - 10*cos(2*pi*x(1)) - 10*cos(2*pi*x(2)) + 20;

主脚本main.m

%% I. 清空环境
clc
clear
close all%% II. 绘制目标函数曲线
figure
[x,y] = meshgrid(-5:0.1:5,-5:0.1:5);
z = x.^2 + y.^2 - 10*cos(2*pi*x) - 10*cos(2*pi*y) + 20;
mesh(x,y,z)
hold on%% III. 参数初始化
c1 = 1.49445;
c2 = 1.49445;kmax = 100;   % 进化次数  
popsize = 100; %种群规模% 控制粒子速度
Vmax = 1;
Vmin = -1;
% 个体位置变化的最大范围
popmax = 5;
popmin = -5;%% IV. 产生初始粒子和速度
for i = 1:popsize% 随机产生一个种群pop(i,:) = 5*rands(1,2);     %初始各个粒子位置V(i,:) = rands(1,2);         %初始化各个粒子速度% 计算适应度fitness(i) = fit(pop(i,:));   %各个粒子的适应度值
end%% V. 个体极值和群体极值
[best_fitness best_index] = max(fitness);
Gbest = pop(best_index,:);     %粒子群体最优位置
Pbest = pop;                   %粒子个体最优位置
fitnessPbest = fitness;        %粒子个体最优适应度值
fitnessGbest = best_fitness;   %粒子群体最优适应度值%% VI. 迭代寻优
for i = 1:kmaxfor j = 1:popsize% 速度更新V(j,:) = V(j,:) + c1*rand*(Pbest(j,:) - pop(j,:)) + c2*rand*(Gbest - 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;% 适应度值更新fitness(j) = fit(pop(j,:)); endfor j = 1:popsize  % 个体最优更新if fitness(j) > fitnessPbest(j)Pbest(j,:) = pop(j,:);fitnessPbest(j) = fitness(j);end% 群体最优更新if fitness(j) > fitnessGbestGbest = pop(j,:);fitnessGbest = fitness(j);endend Best(i) = fitnessGbest;%用于存储每次迭代产生的最优的适应度值        
end
%% VII.输出结果
[fitnessGbest, Gbest]
plot3(Gbest(1), Gbest(2), fitnessGbest,'bo','linewidth',1.5)figure
plot(Best)
title('最优个体适应度','fontsize',12);
xlabel('进化代数','fontsize',12);ylabel('适应度','fontsize',12);

运行结果：

ans =80.7065   -4.5223    4.5231

适应度函数fun.m

function y = fun(x)
%函数用于计算粒子适应度值
%x：输入粒子 
%y：粒子适应度值 
y = 0.5 + (sin(sqrt(x(1)^2+x(2)^2))^2 - 0.5) / (1 + 0.001*(x(1)^2+x(2)^2))^2;
% y = 3*cos(x(1)*x(2)) + x(1) + x(2)^2;

主脚本main.m

%% I. 清空环境
% clc
% clear%% II. 绘制目标函数曲线
figure
[x,y] = meshgrid(-10:0.1:10,-10:0.1:10);
z = 0.5 + (sin(sqrt(x.^2+y.^2)).^2 - 0.5)./ (1 + 0.001*(x.^2+y.^2)).^2;
mesh(x,y,z)
hold on
%% III. 参数初始化
c1 = 1.49445;
c2 = 1.49445;kmax = 100;   % 进化次数
popsize = 100; %种群规模% 控制粒子速度
Vmax = 1;
Vmin = -1;
% 个体位置变化的最大范围
popmax = 10;
popmin = -10;%惯性权重最大值和最小值
Wmax = 0.9;
Wmin = 0.4;% %压缩因子
% phi = c1+c2;
% lamda = 2 / abs(2 - phi - sqrt(phi^2 - 4*phi));   % 压缩因子
%% IV. 产生初始粒子和速度
for i = 1:popsize% 随机产生一个种群pop(i,:) = 10*rands(1,2);     %初始各个粒子位置V(i,:) = rands(1,2);         %初始化各个粒子速度% 计算适应度fitness(i) = fun(pop(i,:));   %各个粒子的适应度值
end%% V. 个体极值和群体极值
[best_fitness best_index] = min(fitness);
Gbest = pop(best_index,:);     %粒子群体最优位置
Pbest = pop;                   %粒子个体最优位置
fitnessPbest = fitness;        %粒子个体最优适应度值
fitnessGbest = best_fitness;   %粒子群体最优适应度值%% VI. 迭代寻优
for i = 1:kmax%计算动态惯性权重w = Wmax - (Wmax - Wmin)*i / kmax;for j = 1:popsize% 速度更新V(j,:) = w*V(j,:) + c1*rand*(Pbest(j,:) - pop(j,:)) + c2*rand*(Gbest - pop(j,:));%标准粒子群算法%         V(j,:) = lamda*V(j,:) + c1*rand*(Pbest(j,:) - pop(j,:)) + c2*rand*(Gbest - pop(j,:));%压缩粒子群算法%         V(j,:) = V(j,:) + c1*rand*(Pbest(j,:) - pop(j,:)) + c2*rand*(Gbest - 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;% 适应度值更新fitness(j) = fun(pop(j,:));endfor j = 1:popsize% 个体最优更新if fitness(j) < fitnessPbest(j)Pbest(j,:) = pop(j,:);fitnessPbest(j) = fitness(j);end% 群体最优更新if fitness(j) < fitnessGbestGbest = pop(j,:);fitnessGbest = fitness(j);endendBest(i) = fitnessGbest;%用于存储每次迭代产生的最优的适应度值
end
%% VII.输出结果
disp(['当x=',num2str(Gbest(1)),'和','y=',num2str(Gbest(2)),'时','所求函数的最小值是',num2str(fitnessGbest)])plot3(Gbest(1), Gbest(2), fitnessGbest,'bo','linewidth',15)figure
plot(Best)
title('最优个体适应度','fontsize',12);
xlabel('进化代数','fontsize',12);ylabel('适应度','fontsize',12);

运行结果：

>> main
当x=7.2595e-06和y=-1.7485e-05时所求函数的最小值是3.5877e-10