粒子群优化算法（1）

背景：

粒子群优化算法起源于对鸟群觅食活动的分析。鸟群在觅食的时候通常会毫无征兆的聚拢，分散，以及改变飞行的轨迹，但是在不同个体之间会十分默契的保持距离。所以粒子群优化算法模拟鸟类觅食的过程，将待求解问题的搜索空间看作是鸟类飞行的空间，将每只鸟抽象成一个没有质量和大小的粒子，用这个粒子来表示待求解问题的一个可行解。所以，寻找最优解的过程就相当于鸟类觅食的过程。

粒子群算法也是基于种群以及进化的概念，通过个体间的竞争与协作，实现复杂空间最优解的求解。但是与遗传算法不同的是，他不会对每个个体进行“交叉”，“变异”等操作，而实以一定的规则，更新每个粒子的速度以及位置，使得每一个粒子向自身历史最佳位置以及全局历史最佳位置进行移动，从而实现整个种群向着最优的方向进化。

算法的具体描述：

在粒子群优化算法中，粒子之间通过信息共享机制，获得其它粒子的发现与飞行经历。粒子群算法中的信息共享机制实际上是一种合作共生的行为，在搜索最优解的过程中，每个粒子能够对自己经过的最佳的历史位置进行记忆，同时，每个粒子的行为有会受到群体中其他例子的影响，所以在搜索最优解的过程中，粒子的行为既受其他粒子的影响，有受到自身经验的指导。

粒子群优化算法对于鸟群的模拟是按照如下的模式进行的：假设一群鸟在空中搜索食物，所有鸟知道自己当前距离食物有多远（这里的远近会用一个值来衡量，适应度值），那么每只鸟最简单的搜索策略就是寻找距离目前距离食物最近的鸟的周围空间。因此，在粒子群算法中，每个粒子都相当于一只鸟，每个粒子有一个适应度值，还有一个速度决定他们的飞行的距离与方向。所有的粒子追随当前最优的粒子在解空间中搜索。每搜索一次，最优的粒子会发生变化，其他的粒子又会追随新的最优粒子进行搜索，如此反复迭代。

在迭代开始的时候，每个粒子通过随机的方式初始化在空间中的速度和位置，然后在迭代过程中，粒子通过跟踪两个极值来自己在解空间中的位置和速度，一个极值是单个粒子自身在迭代的过程中的最优位置（就是最优适应度值所对应的空间解），这个称之为粒子的个体极值。另一个极值是种群中所有的粒子在迭代过程中所找到的最优位置，这个成为全局极值。如果粒子只是跟踪一个极值的话，则算法称为局部粒子群算法或者全局粒子群算法。

基本的粒子群算法：

假设搜索空间的维度是D,空间中粒子的个数是N，则粒子的位置可以表示为一个向量：

$X_{i}=(x_{i,1},x_{i,2},\cdots ,x_{i,D})$

粒子的飞行速度也可以表示为一个向量：

$V_{i}=(v_{i,1},v_{i,2},\cdots ,v_{i,D})$

第i个粒子自身的最优个体极值可以表示为：

$p_{best}=(p_{i,1},p_{i,2},\cdots ,p_{i,D})$

整个粒子群的全局最优极值可以表示为：

$g_{best}=(g_{1},g_{2},\cdots ,g_{D})$

在算法的迭代过程中，所有的粒子按照如下的方式更新自己的速度以及位置:

$v_{i}(t+1)=v_{i}(t)+c_{1}r_{1}(t)[p_{i,best}(t)-x_{i}(t)]+c_{2}r_{2}(t)[g_{best}(t)-x_{i}(t)]$

$X_{i}(t+1)=X_{i}(t)+v_{i}(t+1)$

其中：

$c_{1},c_{2}$ 表示学习因子

$r_{1},r_{2}$ 表示[0, 1]范围内的均匀分布的随机数

$v_{i}\in [v_{min}, v_{max}]$ 表示粒子的速度，其中 $v_{min}, v_{max}%u9700%u8981%u624B%u52A8%u8BBE%u5B9A$ 需要手动设定

从上面的表达式可以看到，粒子的速度主要由三部分组成，第一部分可以看作是粒子的惯性，反映了粒子具有维持自身速度不变的趋势。第二部分代表了粒子的认知属性，表示粒子具有向自身历史最优位置移动的属性。第三部分代表粒子的“社会属性”,反映粒子之间通过信息共享，来分享种群中的最优群体历史经验。代表了种群中的粒子有向全局最佳位置逼近的趋势。这部分应该算是粒子群优化算法的核心思想。

对于群智能进化算法，算法的核心思想无非就是两点，1.算法的探索能力，总结起来说，就是算法向未知区域开拓的能力，也即算法的全局搜索能力。2.算法的局部搜索能力，指的是粒子在原来搜索的轨迹上进行更细一步的搜索，是算法在原来位置位置的较小领域内进行搜索的能力。

所以在上面的粒子群所发的基础上，提出了带有惯性权重的改进粒子群算法：

$v_{i}(t+1)=w v_{i}(t)+c_{1}r_{1}(t)[p_{i,best}(t)-x_{i}(t)]+c_{2}r_{2}(t)[g_{best}(t)-x_{i}(t)]$

惯性权重 $w$ 表示在多大程度上表示原来的权重， $w$ 越大，表示全局的搜索能力越强，局部搜索能力越弱，反之，则表示全局搜索能力越弱，局部搜索能力越强。

经过实验表明， $w$ 取值为0.8~1.2之间的时候，粒子群算法由更快的收敛速度。

粒子群算法的改进：

上面介绍了带有惯性权重的粒子群优化算法，对于惯性权重的设计，也可以采用动态的方法进行。在迭代开始的时候，可以设定一个较大的惯性权重，使得算法具有较强的全局搜索能力。随着算法迭代次数的增加，可以采用线性或者非线性的方式减小惯性权重系数。保证粒子能够具有较强的局部搜索能力，在极值点周围进行细致的全局搜搜。现在采用较多的调节权重的方法是线性递减权值策略：

$w = w_{max}-\frac{(w_{max}-w_{min}).t}{T_{max}}$

一般取：

$w_{max}=0.9,w_{min}=0.4$

一种收敛速速更快的粒子群算法：

压缩因子粒子群算法是一种收敛速度较快的算法，而且能够的到质量较高的解：压缩因子粒子群算法的速度更新公式为：

$v_{i}(t+1)=\lambda v_{i}(t)+c_{1}r_{1}(t)[p_{i,best}(t)-X_{i}(t)]+c_{2}r_{2}(t)[p_{g}(t)-X_{i}(t)]$

其中， $\lambda$ 为压缩因子：

$\lambda =\frac{2}{|2-\varphi -\sqrt{\varphi ^{2}-4\varphi }|}$

$\varphi =c_{1}+c_{2}$

-------------------------------------------------------------------------------------------------------------------------------------------------------------

粒子群优化算法的一般流程：

在粒子群算法中，一些参数的选取对于算法的性能会产生较大的影响：
1. 种群规模NP

一般取100-200之间，不能过大，也不宜过小

2. 惯性权重

控制算法的全局搜索能力以及局部搜索能力。一般可以取值[0.8, 1.2]

3. 加速常数 $c_{1},c_{2}$ （也成为学习因子）

分别调节粒子向历史最优位置和全局最优位置移动的步长。分别决定了粒子的个体经验和粒子的群体经验对粒子运行轨迹的影响。通常取 $c_{1}=c_{2}=1.5$

4. 粒子的最大速度 $v_{max}$ ：

用来对例子的速度进行限制。粒子的速度限制范围不宜过大或者过小。研究指出，调节 $v_{max}$ 和调节惯性权重的作用是一样的，所以一般只有 $v_{max}$ 限制速度的大小，通过调节惯性权重系数来调节粒子群算法的全局搜索能力和局部搜索能力。

5. 停止准则：

在前面介绍遗传算法的时候，一般采取的停止准则是最大迭代次数，最大停止准则还可以是计算精度，或者最优解的停止步数。

6. 边界条件处理：

与遗传算法中的处理方法一样

程序算例：

假设需要求解多项式的最小值，多项式的表达式如下所示：

$f(x)=\sum_{i=1}^{n}x_{i}^{2}$

其中，n=20; x的范围是[-20, 20]

压缩因子粒子群算法：

clear all;
clc
N=100;           % 个体数
D=10;            % 每个个体的维数
gen=50;          % 进化代数
c1=1.2;          %cognition
c2=1.5;          % social
w=0.8;           % interial
X_max=20;
X_min=-20;
V_max=10;
V_min=-10;
%% 压缩粒子群算法
phi=c1+c2;
lambda=2/abs(2-phi-sqrt(phi^2-4*phi));
%%%      初始化种群个体     %%%
x=X_min+rand(N,D)*(X_max-X_min);
v=V_min+rand(N,D)*(V_max-V_min);
%%  初始化局部最有位置和最优值
partial_best=x;   %  局部最优值% 全局最优位置和最优值
pbest=[N,1];
for i=1:Npbest(i)=func(x(i,:));
end[minum,index1]=min(pbest);
global_best=x(index1,:);            % 全局最优位置
generation_best=zeros(gen,1);       % 每一次迭代中的最优值
% 迭代部分
for i=1:genfor j=1:Nif (func(x(j,:)))<pbest(j)      % 个体历史最优partial_best(j,:)=x(j,:);   % 最优位置pbest(j)=func(x(j,:));      % 最优速度endend%% 全局最优解和全局最优值[minum,index1]=min(pbest);        global_best=partial_best(index1,:);generation_best(i)=minum;       % 每一代中的最小值%% 更新速度和位置信息for j=1:Nv(j,:)=v(j,:)*lambda+c1*rand()*(partial_best(j,:)-x(j,:))+c2*rand()*(global_best-x(j,:));   x(j,:)=x(j,:)+v(j,:);end%% 边界条件处理（两种方法1,2）for i=1:Nfor j=1:Dif (v(i,j)<V_min)||(v(i,j)>V_max)v(i,j)=V_min+rand()*(V_max-V_min);end           if (x(i,j)<X_min)||(x(i,j)>X_max)x(i,j)=X_min+rand()*(X_max-X_min);endendendendfigure(1);
plot(generation_best);
title('寻优曲线');
xlabel('进化代数');
ylabel('适应度');%   适应度函数
function f=func(buf)f=0;for i=1:10f=f+buf(i)*buf(i);end
end

运行结果：

测试函数：
计算函数的最小值，函数的表达式如下所示：

$f(x)=3cos(x)+x+y^{2}$

其中x,y的范围是[-4, 4]

函数图像如下所示：

代码：

% 绘制函数的图像
X=-4:0.1:4;
Y=-4:0.1:4;
[x,y]=meshgrid(X,Y);
z=3.*cos(x.*y)+x+y.^2;
figure(1);
surf(x,y,z);
title('函数图像');
xlabel('X');
ylabel('Y');
zlabel('Z');
figure(2);
%%
% f=3*cos(xy)+x+y*y 的最小值
% 采用动态的惯性权重系数%%
% clearvars -EXCEPT count;
clc;
N=100;      % particle number
D=2;        % munber of variable
gen=50;     % iternation number;
c1=1.5;
c2=1.5;
w_min=0.4;  % 速度最大最小值
w_max=0.8;
X_min=-4;
X_max=4;
V_min=-1;
V_max=1;
%% 初始化种群个体  位置和速度信息
x=X_min+rand(N,D)*(X_max-X_min);
v=V_min+rand(N,D)*(V_min-V_max);
partial_best=x;                % 局部最优
global_best=ones(1,D);         % 全局最优
pbest=ones(N,1);
for i=1:Npbest(i)=func(x(i,:));
end
% 寻找全局最优点
[minum,index]=min(pbest);
global_best=x(index,:);
SIGN=inf;
for i=1:gen% 计算每个粒子的历史最优位置for j=1:Nif func(x(j,:))<pbest(j)partial_best(j,:)=x(j,:);pbest(j)=func(partial_best(j,:));endend%%%   动态惯性权值  %%%w=w_max-((w_max-w_min)/gen)*i;  %% 寻找全局最优值[minum,index]=min(pbest);global_best=partial_best(index,:);gen_best(i)=minum;%% 更新粒子群的速度和位置信息for m=1:Nv(m,:)=w*v(m,:) + c1*rand()*(partial_best(m,:)-x(m,:)) + c2*rand()*(global_best-x(m,:));x(m,:)=x(m,:) + v(m,:);end%% 边界条件处理for m=1:Nfor n=1:Dif (x(m,n)<X_min)||(x(m,n)>X_max)x(m,n)=X_min+rand()*(X_max-X_min);endif (v(m,n)<V_min)||(v(m,n)>V_max)v(m,n)=V_min+rand()*(V_max-V_min);endendendendfigure(2);
plot(gen_best);
title('进化曲线');
xlabel('迭代次数');
ylabel('适应度');function f=func(buf)f=3*cos(buf(1)*buf(2))+buf(1)+buf(2)^2;
end

运行结果：