【Matlab群体智能算法第一期】粒子群算法及其变体（一）

一、基本粒子群算法原理

        上世纪末 Eberhart 和 Kennedy 首先提出用粒子群算法模拟鸟群的捕食行为。假设在一片空间内，固定一块食物的位置，鸟群事先并不知晓其存放地点，他们随机搜索食物，因此它们需要通过不断更新自己的速度和位置来靠近食物邻近的那只鸟以获得这块固定的食物。

        在 PSO 算法模型中，都是通过搜索区域中的单只鸟来解决优化问题的。单只鸟在搜索区域中以一种速度搜索，将对比自身与鸟群其余个体的适应度评价后改变飞行的方向和距离。种群中的每只鸟都是解空间中的一个“粒子”，种群代表所有粒子的集合。捕食过程的最终目的就是找到搜索区间内的食物，需要通过实时共享自己当前的位置追随当前的最优粒子，也就是与最优值适应度最高的个体，这样就能快速高效地找到食物。PSO 算法在每次迭代过程中，粒子通过个体之间的相互合作和群体中的信息共享，追踪个体极值与全局极值来更新自己的位置。图1为粒子群算法中粒子解移动的示意图。

图1 粒子群算法中粒子解移动示意图

         在传统粒子群算法中，其主要分为两大部分，分别为粒子种群初始化及种群位置的更新。而在种群初始化方面，粒子群算法对粒子的速度及位置进行初始化操作，初始化公式如下：

         上述更新方法中，方式一主要出现与经典粒子群算法中，而方式二在当前主流群体智能算法及PSO算法改进版本中较为常见。

        与方式一相比，采用方式二进行初始化时，粒子在初始化阶段即可防止种群解越限，提升了解的质量（这一方式对约束优化问题更为关键）。

        粒子群算法在种群初始化后，对初始种群进行排序，筛选出当前种群下粒子群算法最优适应度值及所对应的粒子解的位置。其Matlab参考代码如下：

[fitnessgbest bestindex]=min(fitness);
gbest=x(bestindex,:);

        对粒子进行初始化过程后，粒子群算法进入更新迭代过程。迭代过程中判别算法是否迭代完成主要有两个条件，其一是算法迭代过程中是否满足精度要求，其二是算法是否达到最大迭代次数。在实际求解过程中，多数问题均以是否达到最大迭代次数最为判别标准。

        同时，粒子群算法更新迭代的机理较为简单，便于实现。其主要原则为保留本次迭代过程中目标函数最小值所对应的种群解作为当前迭代过程中的最优解，并保存粒子的位置及目标函数值。对应的matlab代码如下：

for it=1:itex
%种群更新公式
(见下文)

%PSO种群解选择机制
for i=1:N
    if fitness(i)<fitness_best
        fitness(i)=fitness_best;%保留当前粒子解对应的目标函数值
        gbest=x(i,:);%将当前粒子位置作为最优解
    end
end

end

        根据上述种群选择机制，粒子群算法进行优化过程中最优种群选择。在结束当前种群对比，获取对应最优解后，粒子群算法进行种群位置的更新，以搜索更优粒子解。

        粒子群算法更新公式如下：

         上述粒子群算法种群更新公式中，主要分为三部分，即：记忆项、自身认知项与群体认知项。在记忆项中，惯性权重w的值影响算法寻优性能。当w值较大时，粒子群算法的全局搜索能力增强，但会削弱算法的局部搜索能力。当w值较小时，算法的局部搜索能力增强，但全局搜索能力降低。因此，惯性权重w的选择会对算法的寻优能力产生较大的影响，也因此针对该部分的改进策略是PSO算法的改进方向之一。

        根据上述理论，粒子群算法求解流程如下：

（1）算法参数初始化，设置PSO算法运行的基本参数，如惯性权重w，学习因子、种群数量、最大迭代次数等。

（2）种群初始化，根据PSO算法种群初始化公式，获取算法初始时刻粒子种群的位置。

（3）计算目标函数值，并获取全局最优值；

（4）更新PSO算法粒子速度及位置

（5）计算更新后种群适应度值，更新粒子历史最优位置

（6）更新群体全局最优位置

（7）算法收敛性判断，若满足结束条件，则寻优结束，获取最优粒子位置及适应度值；若为满足结束条件，则返回步骤（4）。

根据上述理论及算法求解流程，标准PSO算法参考代码（matlab语言）如下：

%% 标准粒子群算法主程序
clc;
clear all;
close all
%% 算法基本参数设置
c1=2; %学习因子1
c2=2;%学习因子2
w=0.7;%惯性权重
MaxDT=500;%最大迭代次数
D=20;%搜索空间维数（未知数个数）
N=30;%初始化群体个体数目
Vmax=1;%速度上限
Vmin=-1;%速度下限
popmax=100;%粒子上边界
popmin=-100;%粒子下边界
%% 种群初始化
for i=1:N
        pop(i,:)=popmin+(popmax-popmin)*rand(1,D);  %随机初始化位置

        V(i,:)=rand(1,D); %随机初始化速度
        fitness(i)=fitness_pso(pop(i,:));

    end
%计算各个粒子的适应度值并初始化Pi和Pg
[fitnessgbest bestindex]=min(fitness);
gbest=pop(bestindex,:);
pbest=pop;
fitnesspbest=fitness;

%% 粒子群算法更新迭代部分
for i=1:MaxDT
    for j=1:N
        %种群更新
        V(j,:)=w*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)=fitness_pso(pop(j,:));
  
       %个体极值更新
       if fitness(j)<fitnesspbest(j)
            pbest(j,:)=pop(j,:);
            fitnesspbest(j)=fitness(j);
       end

       %全局极值更新
       if fitness(j)<fitnessgbest
           gbest=pop(j,:);
           fitnessgbest=fitness(j);  
       end
       
    end
   %记录粒子全局最优解
   Fgbest(i)=fitnessgbest;
   
end
%% 结果可视化
figure
plot(Fgbest)
title(['适应度曲线 ' '终止次数=' num2str(MaxDT)]);
xlabel('进化代数');
ylabel('适应度')

二、PSO算法案例应用

        本部分通过测试函数实现PSO算法的求解过程，本文所用测试函数如下：

 其中，函数种群的上下限值分别为-100,100，维度20、30、60、100。具体计算结果如下：

（1）粒子维度为20时：

（2）粒子维度为30时：

（3）粒子维度为60时：

（4）粒子维度为100时：

        根据上述求解结果，标准粒子群算法在该测试函数下，维度为20、30、60、100时，其目标函数值分别为：0.0103、0.0397、1.7937、11.5576。

        由此可见，标准粒子群算法在处理高维度求解问题时，其求解效率与寻优能力逐渐降低。如何提升高维度下算法的寻优能力，亦是后续本专栏的重点。

        为了便于学习交流，本文给出的标准PSO算法所用测试函数的matlab代码如下：

function f=fitness_pso(x)
%写法一
f=sum(x.^2);

%写法二（适合初学者）
% f2=0;
% for L=1:size(x,2)
%     f2=f2+x(1,L)^2;
% end

三、本专栏后续内容预告

（1）针对PSO算法的改进之处进行分析与matlab代码复现

（2）粒子群算法中高维约束问题的matlab代码实现

（3）基于PSO算法的微电网系统调度案例分析与matlab代码分享

（4）其余群体智能优化算法及其改进（matlab）