数学建模智能算法--粒子群算法（PSO）

概念

通过模拟鸟群觅食过程中的迁徙和群聚行为而提出的一种基于群体智能的全局随机搜索算法。与其他智能算法（比如蚁群算法，遗传算法），它也是基于种群或者进化的概念，通过个体之间的协作与竞争，实现复杂空间最优解的搜索。不同之处是PSO模拟社会。将每个可能产生的解表述为群中的一个微粒，每个微粒都具有自己的位置向量和速度向量，以及一个由目标函数决定的适应度。所有微粒在搜索空间中以一定速度飞行，通过追随当前搜索到的最优值来寻找全局最优值。在算法运行的过程中，当前个体最优值以及全局最优值都在每一次迭代过程中不断更新。

它根据以下三条规则来进行操作：

1. 飞离最近的个体，以避免碰撞
2. 飞向目标
3. 飞向群体的中心

基本控制参数

粒子种群规模

一般设置为 20-50，粒子种群是随机初始化的，粒子种群数目越多，在整个全局范围里面的分布越分散，搜索的空间范围越大，更容易发现全局最优解。当然，粒子种群规模越大，运行时间也就越长

惯性权重

惯性权重值表示对粒子当前速度继承的多少。权重值较小的时候，不容易从当前区域跳到全局最优的区域，当跳到局域最优的能力强。权重值较大的时候，局域寻优能力弱，全局寻优能力强。

权重的选择有两种，固定权重和时变权重

固定权重是指选择一个常数作为惯性权重值，进化过程中保持不变，一般选取为[0.8,1,2]

时变权重是指设定一个惯性权重的区间，在进化的过程中按照某种方法逐渐减小惯性权重。因为随着进化的过程，粒子的活跃范围逐渐变小，所以权重也很有必要变小以减弱粒子的全局寻优能力。

加速常数

加速常数有两个c1 和 c2，表征着个体经验和群体经验对粒子的影响，c1为0的话，说明个体经验对粒子没有作用，只有社会经验，收敛速度较快，但容易陷入局部最优。c2为0的话，说明个体之间没有信息的交互，找到最优解的概率较小。如果c1和c2均为零的话，我们也许可以认为这是一个蒙泰卡罗模拟，没有搜索。我们一般设置二者相等，这样二者就可以有相同的影响力。

粒子最大速度

粒子的速度在每一个维度上都需要有一个速度的限制范围，这就决定了粒子对空间搜索的力度。同样，参数设置的过大过小都是有弊端的。过大粒子会飞过优秀区域，过小无法对局部最优之外的区域进行搜索，无法跳出局部最优。

停止准则

即算法的停止条件，比如迭代次数，计算的精度等

边界条件

对于粒子的一个或者多个维度，我们需要对粒子的速度或者位置进行界定，以满足我们求解问题的条件，从而将粒子的位置限定在可行搜索范围内，减少不必要的搜索，提高搜索效率。

我们可以设置速度或者位置的范围，当例子超过范围的时候，我们可以在范围内产生一个随机数进行代替，或者将其设置为最大值。

运算流程

无变量之间约束的PSO函数寻优

这里以Griewank函数为例，求解最优值。

本函数在x = ( 0 , 0 ,…, 0 )处有全局极小点0。x可以是n维，可以在下列代码中通过 size 来指定参数的个数。这里设定为5，就是说这个函数为5维的函数

我们在这里指定每个参数的范围为[-5,5]，没有添加变量之间的约束，比如x1 + x2 > x3 这种约束。

function y=Griewank(x)
col = size(x,2);
y1 = 1 / 4000 * sum(x.^2 );
y2 = 1 ;
for  h = 1 :col
    y2 = y2 * cos(x(h) / sqrt(h));
end
y = y1 - y2 + 1 ;
end

pso 求解

%% 清空环境
clc % 清屏
clear all; % 删除工作区变量
close all; % 关掉图形窗口
%% 参数初始化
%两个加速参数表征个体经验和群体经验，一般设置为 相同
c1 = 1.4;
c2 = 1.4;
% 进化次数  
maxg=50;   
%种群规模
sizepop=200;   

% 边界条件的设定 
Vmax=1;
Vmin=-1;
popmax=5;
popmin=-5;

% 粒子中参数的个数
size = 5;

%% 产生初始粒子和速度
for i=1:sizepop
    pop(i,:)=5*rands(1,size);    %初始种群对于不同维度的参数范围，设定为不同，这里均默认为[-5,5]
    V(i,:)=rands(1,size);        %初始化速度
    fitness(i)=Griewank(pop(i,:));   %染色体的适应度
end

%找最好的个体
[bestfitness bestindex]=min(fitness);
zbest=pop(bestindex,:);   %全局最佳
gbest=pop;                %个体最佳
fitnessgbest=fitness;     %个体最佳适应度值
fitnesszbest=bestfitness; %全局最佳适应度值

%% 迭代寻优
for i=1:maxg
    for j=1:sizepop
        %速度更新
        V(j,:) = V(j,:) + c1*rand*(gbest(j,:) - pop(j,:)) + c2*rand*(zbest - pop(j,:));
        V(j,find(V(j,:)>Vmax))=Vmax;
        V(j,find(V(j,:)<Vmin))=Vmin;
        %种群更新
        pop(j,:)=pop(j,:)+0.5*V(j,:);
        pop(j,find(pop(j,:)>popmax))=popmax;
        pop(j,find(pop(j,:)<popmin))=popmin;
        %自适应变异
        if rand>0.8
            k=ceil(size*rand);
            pop(j,k)=rand;
        end
        %适应度值
        fitness(j)=Griewank(pop(j,:));
        %个体最优更新
        if fitness(j) < fitnessgbest(j)
            gbest(j,:) = pop(j,:);
            fitnessgbest(j) = fitness(j);
        end
        %群体最优更新
        if fitness(j) < fitnesszbest
            zbest = pop(j,:);
            fitnesszbest = fitness(j);
        end  
    end
    yy(i)=fitnesszbest;      
end

%% 执行结果
plot(yy,'Linewidth',2)
title(['适应度曲线  ' '终止代数＝' num2str(maxg)]);
grid on
xlabel('进化代数');ylabel('适应度');
% 结果输出
disp('最优结果参数：');
disp(zbest);
disp('最优结果：');
disp(Griewank(zbest));

执行结果如下：

最优结果参数：
    0.0093   -0.0127   -0.0331    0.0131   -0.0188

最优结果：
   3.2364e-04

我们可以通过增加迭代次数，加速参数，种群规模来找到都更加精确的值

%% 参数初始化
%两个加速参数表征个体经验和群体经验，一般设置为 相同
c1 = 1.49;
c2 = 1.49;
% 进化次数  
maxg=500;   
%种群规模
sizepop=2000;   

结果如下：

最优结果参数：
    0.0002    0.0013    0.0002    0.0024    0.0001
最优结果：
   1.1659e-06

另附有约束条件的粒子群算法