PSO算法与MATLAB实现

1. PSO算法基础理论

决策过程中,人们常常会选择综合两种重要信息做出决策:一则是自己的主观经验,二则是其他人的经验。

类比鸟类,群鸟觅食的时候,每只鸟的初始状态都是出于随机位置,而且飞翔状态也是随机的。但是随着时间的推移这些初始处于随机位置的鸟类通过群内相互学习、信息共享和个体不断积累觅食的经验,自发组织积聚成为一个群落。局部最优即每只鸟能够记住自己找到的最好的位置。此外,还能记住群鸟中所有个体能够找到的最好位置,称为全局最优

整个鸟群的觅食中心都是去向全局最优移动的,在生物学中称为“同步效应”。

1.1 模型化

在群鸟觅食模型中,每个个体都可以被看成是一个粒子,鸟群则被看成一个粒子群。设在一个$D$维的目标空间里,有$m$个粒子组成的群里,其中第$i$个粒子的位置表示为$X_i$,即第$i$个粒子在$D$维搜索空间的位置是$X_i$。换言之,每个粒子的位置就是一个潜在解(意味着可能是你要找到的那个全局位置或者局部位置)。将$X_i$代入目标函数既可以计算出其适应值,根据适应值的大小衡量其优劣,粒子个体经历过的最好位置记为$P_i$,整个群体所有粒子经历过的最好位置记为$P_g$,粒子$i$的速度记为$V_i$。

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

$$ v_i=w*v_i+c_1*r_1*(p_i-x_i)+c_2*r_2*(p_g-x_i)*x_i=x_i+a*v_i $$

  • $w$为非负数,称为惯性因子
  • $c_1$和$c_2$称为加速常数。$c_1$是根据个体自身的经验进行判断的常数;$c_2$则是根据群体的经验
  • $r_1$和$r_2$为[0,1]范围内变换的随机数
  • $a$称为约束因子,目的是控制速度的权重
  • $v_i$即粒子$i$的飞翔速度$V_i$被一个最大速度$V~max~$所限制。如果当前时刻粒子在某纬度的速度$V_i$更新后超过改纬度的最大飞翔速度$V~max~$,则当前时刻纬度被限制为$V~max~$

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

1.2 PSO算法的约束优化

原则上PSO适用于非约束的优化问题(也就是没有约束条件的优化问题),但是近年来,有许多人将其推广到约束化问题,主要在于怎么处理好约束条件。基于PSO算法约束优化工作主要分为两类:

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

1.3 PSO的优缺点

优点:

  1. 不依赖于问题信息,采用实数求解,算法通用性比较强
  2. 调整的参数较少,原理简单,容易实现
  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. 输出

2.2 PSO算法的参数选取

参数1:粒子数

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

参数2:惯性因子

惯性因子$w$对于粒子群算法的收敛性起到很大作用,$w$值越大,粒子飞翔幅度越大,容易错失局部寻优能力,而全局搜索能力越强(越容易跳出局部最优点);反之,则局部能力越强,全局弱。

这里有个误区:就是惯性因子越大,其当前飞翔速度越大是错的,这里只是较小程度地离开原先的寻优轨道。

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

参数3:加速常数$c_1$和$c_2$

一般情况下取$c_1=c_2=2.0$。目前对于加速常数的取值,有点争议。

其中:

  • $c_1$是代表自身经验
  • $c_2$是代表群体经验
参数4:最大飞翔速度$V~max~$

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

2.3 PSO算法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 秒。
Last modification:August 15th, 2019 at 12:31 am
如果觉得我的文章对你有用,请随意赞赏