在这里插入图片描述

作者：非妃是公主
专栏：《智能优化算法》
博客地址：https://blog.csdn.net/myf_666
个性签：顺境不惰，逆境不馁，以心制境，万事可成。——曾国藩

文章目录

专栏推荐
序
一、概论
二、粒子群算法原理
1. 算法原理
2. 机理分析

三、算法流程
1. 算法伪代码
2. 算法流程图

四、参数设置
1. 粒子种群规模 $N$
2. 惯性权重 $\omega$
Ⅰ. 固定权重
Ⅱ. 时变权重

3. 加速常数 $c_1$ 和 $c_2$
4. 粒子最大速度 $v_{max}$
5. 停止准则
6. 邻域结构的设定
7. 边界条件处理

五、仿真实例
1. 问题
2. 分析
3. matlab代码实现
4. 效果展示

六、算法改进
1. 压缩因子粒子群算法
2. 离散粒子群算法

the end……

专栏推荐

专栏名称	专栏地址
软件工程	专栏——软件工程
计算机图形学	专栏——计算机图形学
操作系统	专栏——操作系统
软件测试	专栏——软件测试
机器学习	专栏——机器学习
数据库	专栏——数据库
算法	专栏——算法

序

粒子群算法也是一种智能优化算法，该算法作为一种群体智能算法，主要借鉴了鸟群的思想（那为什么不叫鸟群算法，非要叫粒子群算法呢？我不得而知，感兴趣的朋友可以在评论区告诉我，我会标明出处后进行补充）。

一、概论

粒子群优化算法对鸟群捕食行为进行研究：一群鸟在区域中随机搜索食物，所有鸟知道自己距离食物的距离。

在粒子群算法中，每一个潜在解对应一只小鸟，可能是不考虑这只小鸟的重量和大小等物理特性，就抽象它为一个点，所以叫做粒子。

粒子群算法首先在给定的空间中随机初始化粒子群，待优化问题的变量决定了解空间的维数。每个粒子有了初始位置和速度后，通过迭代，我们就可以确定它的下一次的位置。

在每一次迭代中，每个粒子通过跟踪两个“极值”来更新自己在解空间中的空间位置和飞行速度：一个极值就是单个粒子本身在迭代过程中找到的最优解粒子，这个粒子叫做个体极值。另一个极值是种群所有粒子在迭代过程中找到的最优解粒子，叫做全局极值。这个方法叫做全局粒子群算法，也是粒子群算法的 v1 版本，但是由于它考虑了全局种群，因此非常容易陷入到局部最优。

对应的改进方法是，只用其中一部分作为邻居粒子，这些邻居粒子中的极值就是局部极值，用它来代替上面的全局极值，这种方法就叫做局部粒子群算法，也是粒子群算法的 v2 版本。

二、粒子群算法原理

1. 算法原理

在一个D维的目标搜索空间中，N个粒子组成一个种群，其中第i个粒子表示为一个D维的向量：

$X_i=(x_{i1},x_{i2},...,x_{iD}),\quad i=1,2,...,N$

第i个粒子的“飞行”速度也是一个D维的向量，记为：

$V_i=(v_{i1},v_{i2},...,v_{iD}),\quad i=1,2,...N$

迄今为止，第i个粒子的个体最优位置：

$p_{best}=(p_{i1},p_{i2},...,p_{iD}),\quad i=1,2,...N$

迄今为止，整个粒子群的最优位置：

$g_{best}=(g_1,g_2,...,g_D)$

对于每一次迭代，粒子根据如下公式进行更新：

$v_{ij}(t+1)=v_{ij}(t)+c_1r_1(t)[p_{ij}(t)-x_{ij}(t)]+c_2r_2(t)[g_{ij}(t)-x_{ij}(t)]$

$x_{ij}(t+1)=x_{ij}(t)+v_{ij}(t+1)$

其中， $c_1$ 和 $c_2$ 为学习因子，就是2个常数，自己设定；

$r_1,r_2$ 为 $[0, 1]$ 范围内的均匀随机数。

$i = 1, 2, ..., NP;$ 表示种群数量， $j = 1, 2, ..., D$ 表示种群维度。

2. 机理分析

观察速度更新的公式，主要分为三个部分，进行机理分析，如下：

第一部分：惯性、动量部分。保持原有的运动趋势。

第二部分：认知部分。根据以往的经验，忘自己的个体最优值附近进行一定的偏移。

第三部分：社会部分。反映了粒子间的合作与知识共享，向全局最优值附近有一定的偏向。

三、算法流程

1. 算法伪代码

初始化种群。包括：NP， $x_i$ 和速度 $v_i$ 。
计算每个粒子适应度fit[i]
对每个粒子，用它的适应度值fit[i]和个体极值 $p_{best}(i)$ 比较，如果 $fit[i]<p_{best}$ ，那么就更新 $p_{best}$ 。
将fit[i]和全局极值 $g_{best}$ 比较。如果 $fit[i]<g_{best}$ ，则更新 $g_{best}$ 。
迭代更新粒子的速度 $v_i$ 得到 $v_{i+1}$ 。
迭代更新 $x_i$ 得到 $x_{i+1}$
对 $x_{i+1}$ 进行边界处理。
判断算法是否满足终止条件。如果满足，结束算法。不满足，返回步骤2.

2. 算法流程图

四、参数设置

1. 粒子种群规模 $N$

一般设置粒子数为20~50。

对于大部分问题10个粒子即可取得很好的结果；
对于复杂的问题，粒子数量可以取到100或200。

另外，粒子数目越大，搜索空间越大，全局搜索能力更强，运行时间也越长。

2. 惯性权重 $\omega$

Ⅰ. 固定权重

在进化过程中保持不变，一般取值为 $[0.8, 1.2]$ .

Ⅱ. 时变权重

粒子在不同阶段拥有不同的搜索和开发能力。

$\omega=\omega_{max}-\frac{(\omega_{max}-\omega_{min})\cdot t}{T_{max}}$

其中， $T_{max}$ 表示最大进化代数； $\omega_{min}$ 表示最小惯性权重； $\omega_{max}$ 表示最大惯性权重； $t$ 表示当前迭代次数。在大多数应用中 $\omega_{max}=0.9$ ， $\omega_{min}=0.4$

3. 加速常数 $c_1$ 和 $c_2$

一般设置 $c_1=c_2=1.5$ .

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

设定 $v_{max}$ 和调整惯性权重是等效的，所以 $v_{max}$ 一般用于对种群初始化进行设定，即将 $v_{max}$ 设定为每维变量的变化范围，二不再对最大速度进行细致的选择及调节。

5. 停止准则

最大迭代次数、计算精度或最优解的最大停滞部署 $\bigtriangleup t$ .

6. 邻域结构的设定

全局最优解收敛块，但回陷入局部最优。
局部版本粒子群算法，收敛慢，但不容易陷入局部最优。
通常的做法是，全局粒子群算法寻找到最优解的大致范围，然后采用局部粒子群算法在最优点附近进行精细搜索。

7. 边界条件处理

当朝贡国最大位置限制 $x_{max}$ 和最大速度限制 $v_{max}$ 时，在范围内随机产生一个数值代替，或者将其设置为最大值，即边界吸收。

五、仿真实例

1. 问题

求下列函数的最小值，其中x的维数n=10。

$f(x)=\sum_{i=1}^{n}x_i^2,\quad(-20\leq x_i \leq 20)$

2. 分析

函数只有一个极小点 $x = (0, 0, ...0)$ ，理论最小值 $f (0, 0, ...0) = 0$ 。

初始化种群。N=100，粒子维数10，最大迭代次数200

参数设置如下表所示：

参数名称	参数大小
种群大小 $N$	100
维度 $D$	10
最大迭代次数	200
学习因子 $c_1=c_2$	1.5
惯性权重 $\omega$	0.8
$X_{max}$	20
$X_{min}$	20
$V_{max}$	10

利用上述推导描述的算法框架进行实现。

3. matlab代码实现

核心代码如下，我已经逐行注释，可以自行阅读代码流程。

%%%%%%%%%%%%%%%%%粒子群算法求函数极值%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%初始化%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear all;              %清除所有变量
close all;              %清图
clc;                    %清屏
N=100;                  %群体粒子个数
D=10;                   %粒子维数
T=200;                  %最大迭代次数
c1=1.5;                 %学习因子1
c2=1.5;                 %学习因子2
w=0.8;                  %惯性权重
Xmax=20;                %位置最大值
Xmin=-20;               %位置最小值
Vmax=10;                %速度最大值
Vmin=-10;               %速度最小值
%%%%%%%%%%%%%%%%初始化种群个体（限定位置和速度）%%%%%%%%%%%%%%%%
x=rand(N,D) * (Xmax-Xmin)+Xmin; % 初始化粒子位置
v=rand(N,D) * (Vmax-Vmin)+Vmin; % 初始化粒子速度
%%%%%%%%%%%%%%%%%%初始化个体最优位置和最优值%%%%%%%%%%%%%%%%%%%
p=x;                            % 初始化个体最优位置
pbest=ones(N,1);                % 初始化个体最优值
for i=1:N
    pbest(i)=func1(x(i,:));     % 计算个体适应度（最优）值
end
%%%%%%%%%%%%%%%%%%%初始化全局最优位置和最优值%%%%%%%%%%%%%%%%%%
g=ones(1,D);                    % 初始化全局最优位置
gbest=inf;                      % 将全局最好点设置为无穷大
for i=1:N                       % 更新全局最优解
    if(pbest(i)<gbest)          % 如果找到
        g=p(i,:);               % 全局最优位置
        gbest=pbest(i);         % 全局最优值
    end
end
gb=ones(1,T);                   % 追踪 T 次全局最优解
%%%%%%%%%%%按照公式依次迭代直到满足精度或者迭代次数%%%%%%%%%%%%%
for i=1:T                       % 开始迭代
    for j=1:N                   % 遍历所有个体
        %%%%%%%%%%%%%%更新个体最优位置和最优值%%%%%%%%%%%%%%%%%
        if (func1(x(j,:))<pbest(j))
            p(j,:)=x(j,:);              % 更新个体最优位置
            pbest(j)=func1(x(j,:));     % 更新个体最优值
        end
        %%%%%%%%%%%%%%%%更新全局最优位置和最优值%%%%%%%%%%%%%%%
        if(pbest(j)<gbest)
            g=p(j,:);                   % 更新全局最优位置
            gbest=pbest(j);             % 更新全局最优值
        end
        %%%%%%%%%%%%%%%%%跟新位置和速度值%%%%%%%%%%%%%%%%%%%%%
        v(j,:)=w*v(j,:)+c1*rand*(p(j,:)-x(j,:))...
            +c2*rand*(g-x(j,:));        % 更新速度
        x(j,:)=x(j,:)+v(j,:);           % 更新位置
        %%%%%%%%%%%%%%%%%%%%边界条件处理%%%%%%%%%%%%%%%%%%%%%%
        for ii=1:D                               % 遍历所有维度
            if (v(j,ii)>Vmax)  ||  (v(j,ii)< Vmin)
                v(j,ii)=rand * (Vmax-Vmin)+Vmin; % 处理速度边界
            end
            if (x(j,ii)>Xmax)  ||  (x(j,ii)< Xmin)
                x(j,ii)=rand * (Xmax-Xmin)+Xmin; % 处理位置边界
            end
        end
    end
    %%%%%%%%%%%%%%%%%%%%记录历代全局最优值%%%%%%%%%%%%%%%%%%%%%
    gb(i)=gbest;
end
disp('最优个体：')
disp(g);                         % 最优个体         
disp('最优值：')
disp(gb(end));                   % 最优值
figure
plot(gb)
xlabel('迭代次数');
ylabel('适应度值');
title('适应度进化曲线')