所有博文已迁移至个人网站：http://ravenxrz.ink，请勿再留言评论，本文新链接：http://ravenxrz.ink/2019/05/03/particle-swarm-optimization-pso-and-matlab-implementation.html

1.PSO相关知识介绍

1.1PSO算法的基础理论

人们在决策过程中常常会综合两种重要的信息：第一种是他们自己的经验，第二种是其他人的经验。
同样的道理，群鸟在觅食过程中，每只鸟的初始状态都是出于随机位置，且飞翔的方向也是随机的。但是随着时间的推移这些初始处于随机位置的鸟类通过群内相互学习、信息共享和个体不断积累觅食的经验，自发组织积聚成一个群落。每只鸟能够记住自己找到的最好的位置，称之为局部最优。此外，还能够记住群鸟中所有个体索恩能够找到的最好位置，称为全局最优，整个鸟群的觅食中心都是去向全局最优移动的，这在生物学中称为“同步效应”。
在群鸟觅食模型中，每个个体都可以被看成是一个例子，鸟群则被看成一个粒子群。设在一个D维的目标空间里，有m个粒子组成的群里，其中第i个粒子的位置表示为Xi，即第i个粒子在D维搜索空间的位置是Xi。换言之，每个粒子的位置就是一个潜在解，将Xi带入目标函数就可以计算出其适应值，根据适应值的大小衡量其优劣，粒子个体经历过的最好位置记为Pi，整个群体所有粒子经历过的最好位置记为Pg。粒子i的速度记为Vi。

粒子群算法采用下列公式对粒子所在的位置不断更新（单位时间为1）：

vi = w*vi+c1*r1*(pi-xi)+c2*r2*(pg-xi)
xi = xi+a*vi

• w 为非负数，称为惯性因子
• c1和c2称为加速常数。c1是根据个体自身的经验进行判断的常数；c2根据群体的经验
• r1和r2为[0 1]范围内变换的随机数
• a 称为约束因子，目的是控制速度的权重
• vi 即粒子i的飞翔速度Vi被一个最大速度Vmax所限制。如果当前时刻粒子在某纬度的速度Vi更新后超过该纬度的最大飞翔速度Vmax，则当前时刻该纬度被限制为Vmax。

迭代终止条件根据具体问题设定，一般为达到预定迭代次数或者粒子群目前为止搜索到最有位置满足目标函数的最小容许误差。

1.2PSO算法的约束优化

原则上PSO适应于非约束的优化问题。近年来，一些学者将其推广到约束优化问题，其关键在于如何处理好约束。基于PSO算法约束优化工作主要分为两类。

1. 罚函数法。罚函数的目的是将约束优化问题转化为无约束优化问题
2. 将粒子群的搜索范围都限制在条件约束族中、

但是目前有关运用PSO算法来处理多约束优化问题的理论还不成熟。

1.3PSO的优缺点

优点

1. 不依赖与问题信息，采用实数求解，算法通用性强。
2. 需要调整的参数少，原理简单，容易实现，这是PSO的最大优点
3. 协同搜索，同时利用个体局部信息和群体全局信息进行指导搜索。
4. 收敛速度快
5. 更容易飞跃局部最优信息。

缺点

1. 算法局部搜索能力较差，搜索精度不高。
2. 算法不能保证搜索到全局最优解，主要有一下两个方面：
   1. 有时粒子群在俯冲过程中会错失全局最优解。
   2. 应用PSO算法处理高位复杂问题是，算法可能过早收敛。
3. 反酸搜索性能对参数具有一定的依赖性
4. PSO算法是一种概率算法
5. 欠缺的生物学背景

2.PSO算法程序设计

2.1程序设计流程

1. 初始化相关参数
2. 评价每个粒子的初始适应值
3. 将初始适应值作为当前每个粒子的最优值，并记录当前的位置作为局部最优位置。
4. 将最佳初始适应值作为当前全局最优值，并记录当前位置。
5. 依据上文提到的计算速度和位置公式进行计算。这里要注意最大速度限幅度处理。
6. 比较当前适应值与之前的适应值，如果更优则进行更新。
7. 找到当前粒子群的全局最优
8. 重复5-7步知道达到最小误差或者达到最优大迭代次数。
9. 输出

粒子算法.png

2.2PSO算法的参数选取

参数1 粒子数

粒子数一般取值20-40,实验表明，对于大多数30个粒子就够用了。粒子数量越多没搜索范围越大，越容易找到全局最优解，但是程序运行时间越长。

参数2惯性因子

惯性因子w对于粒子群算法的收敛性起到很大作用，w值越大，粒子飞翔幅度越大，容易错失局部寻优能力，而全局搜索能力越强；反之，则局部能力越强，全局弱。这里需要注意一个误区：认为惯性因子越大则继承当前飞翔速度也就越大，所以应该是较小程度地离开原先的寻优轨道。错误理解的原因是，w越大，粒子飞翔的速度和位置更新的幅度也就越大，因此偏离原先原先寻优轨道的程度也就越大。（和中学的学习的加速度对速度理解类似？？）

通常的做法是在迭代开始时讲惯性因子设置得较大，然后在迭代过程中逐步减小。w的取值是[0 1]，如果w取定值，那么建议取0.6-0.75.

参数3加速常熟c1和c2

一般情况下c1= c2=2.0。目前对于加速常数c1和c2的确切取值，学术界有较大的分歧。

其中：

• c1 代表自身经验
• c2 代表群体经验

参数4最大飞翔速度Vmax

粒子群算法是通过调整每一次迭代时每一个粒子在每一维上移动一定的距离进行的，速度的改变是随机的，而不希望不收控制的粒子搜索轨道被扩展到问题空间越来越广阔的距离，并达到无穷。如果要粒子进行有效的搜索，就必须对参数Vmax进行限制。参数Vmax有利于防止搜索范围毫无意义的发散。关于Vmax的确定，需要有一定的先验。为了跳出局部最优，需要较大的寻优补偿，而在接近最优值时，采用更小的步长会更好。如果Vmax不变，通常设置为没维变化范围的10%-20%

2.3PSO算法MATLAB实现例题

maxf(x) = 2.1(1-x+2*x^2)*exp(-x^2/2)

matlab代码如下：

clc;
clear;
close all;
tic;                              %程序运行计时
E0=0.001;                        %允许误差
MaxNum=100;                    %粒子最大迭代次数
narvs=1;                         %目标函数的自变量个数
particlesize=30;                    %粒子群规模
c1=2;                            %每个粒子的个体学习因子，也称为加速常数
c2=2;                            %每个粒子的社会学习因子，也称为加速常数
w=0.6;                           %惯性因子
vmax=0.8;                        %粒子的最大飞翔速度
x=-5+10*rand(particlesize,narvs);     %粒子所在的位置
v=2*rand(particlesize,narvs);         %粒子的飞翔速度
fitness = @(x)(1/(1+(2.1*(1-x+2*x.^2).*exp(-x.^2/2))));%定义适应度函数
f = zeros(1,particlesize);      % 预分配
for i=1:particlesize
    for j=1:narvs
        f(i)=fitness(x(i,j));
    end
end
personalbest_x=x;
personalbest_faval=f;
[globalbest_faval,i]=min(personalbest_faval);
globalbest_x=personalbest_x(i,:);
k=1;
while k<=MaxNum
    for i=1:particlesize
        for j=1:narvs
            f(i)=fitness(x(i,j));
        end
        if f(i)<personalbest_faval(i) %判断当前位置是否是历史上最佳位置
            personalbest_faval(i)=f(i);
            personalbest_x(i,:)=x(i,:);
        end
    end
    [globalbest_faval,i]=min(personalbest_faval);
    globalbest_x=personalbest_x(i,:);
    for i=1:particlesize %更新粒子群里每个个体的最新位置
        v(i,:)=w*v(i,:)+c1*rand*(personalbest_x(i,:)-x(i,:))...
            +c2*rand*(globalbest_x-x(i,:));
        for j=1:narvs    %判断粒子的飞翔速度是否超过了最大飞翔速度
            if v(i,j)>vmax
                v(i,j)=vmax;
            elseif v(i,j)<-vmax
                v(i,j)=-vmax;
            end
        end
        x(i,:)=x(i,:)+v(i,:);
    end
    if abs(globalbest_faval)<E0,break,end
    k=k+1;
end
Value1=1/globalbest_faval-1; Value1=num2str(Value1);
% strcat指令可以实现字符的组合输出
disp(strcat('the maximum value','=',Value1));
%输出最大值所在的横坐标位置
Value2=globalbest_x; Value2=num2str(Value2);
disp(strcat('the corresponding coordinate','=',Value2));
x=-5:0.01:5;
y=2.1*(1-x+2*x.^2).*exp(-x.^2/2);
plot(x,y,'m-','linewidth',3);
hold on;
plot(globalbest_x,1/globalbest_faval-1,'kp','linewidth',4);
legend('目标函数','搜索到的最大值');
xlabel('x');
ylabel('y');
grid on;
toc;

程序输出结果：

the maximum value=5.1985
the corresponding coordinate=-1.1617
时间已过 0.250482 秒。