Matlab应用笔记--粒子群算法

注：本篇随笔依据《Matlab在数学建模上的应用》中第5章介绍来写，主要介绍粒子群算法思想及其Matlab实现
（博客以及Matlab小白，若有不当欢迎指出）

粒子群算法（PSO）简介

PSO属于智能算法，智能算法都属于软计算（动态自适应的求解方式）。
PSO依托群鸟觅食模型（Boid模型）寻找最优值。

粒子群算法的基本理论

群鸟觅食模型中，每只鸟的飞行基于自身经验和群体经验。
Boid模型遵守3个行为准则

1. 冲突避免：避免碰撞。
2. 速度匹配：配合群体中心的移动速度。
3. 群体中心：个体向群体中心移动。

一只鸟是一个粒子，m只鸟的鸟群为m个粒子组成的粒子群。
每个粒子都在D维空间中，第i个粒子的位置表示为$X_{i}=(x_{i}^{1},x_{i}^{1},x_{i}^{2},...,x_{i}^{D})$（相当于函数的一个解）
粒子个体经过的最优位置记为$P_{i}=(p_{i}^{1},p_{i}^{1},p_{i}^{2},...,p_{i}^{D})$
所有粒子经理过的最优位置记为$P_{g}=(p_{g}^{1},p_{g}^{1},p_{g}^{2},...,p_{g}^{D})$
粒子速度记为$V_{i}=(v_{i}^{1},v_{i}^{1},v_{i}^{2},...,v_{i}^{D})$
则粒子速度更新公式为：$v_{i}^{d}=\omega v_{i}^{d}+c_{1}r_{1}(p_{i}^{d}-x_{i}^{d})+c_{2}r_{2}(p_{i}^{d}-x_{i}^{d})$
粒子位置更新公式为：$x_{i}^{d}=x_{i}^{d}+\alpha v_{i}^{d}$

其中
$i=1,2,...,m$
$d=1,2,...,D$
$\omega$为惯性因子（非负数）
$c_{1}、c_{2}$为加速常数（非负数）
$r_{1}、r_{2}$为$[0,1]$内的随机数
$\alpha$为约束因子，用来控制速度的权重
此外还应该有$v_{i}^{d}\in [-v_{max}^{d},v_{max}^{d}]$，$v_{max}^{d}$是粒子在每一维度方向上的最大速度。

PSO算法的优缺点

优点：
1.通用性强。
2.调整的参数少，原理简单。
3.协同搜索，个体局部信息和群体局部信息指导搜索。
4.收敛速度快，对内存和CPU要求低。
5.更容易飞跃局部最优信息。

缺点：
1.局部搜索能力差，搜索精度不高。
2.算法不能绝对保证搜索到全局最优解。

PSO算法实现的流程图

PSO算法的参数影响介绍以及选取建议

各参数影响介绍
1.粒子数$m$：粒子数越多，搜索范围越大，越容易得到全局最优解，算法运行的时间也越长。
2.惯性因子$\omega$：$\omega$越大全局搜索能力越强，但搜索精度越低；反之，$\omega$越小全局搜索能力越弱，但搜索精度越高。
3.加速常数$c_{1}、c_{2}$：$c_{1}$和$c_{2}$是调整自身经验和社会经验在其运动中索契作用的权重。

1）$c_{1}=0$ $c_{2}\neq 0$收敛快，但容易陷入局部最优点。
2）$c_{1}\neq 0$ $c_{2}=0$等价于运行$m$个单粒子，得到最优解的概率非常小。
3）$c_{1}=0$ $c_{2}=0$没有任何经验借鉴，各粒子盲目乱飞，很难得到最优解。

4.最大飞翔速度$v_{max}：若不限制飞行速度（过大），则有可能飘向无穷远处；若飞行速度太小，则很难跳出局部最优点。

各参数取值建议
1.粒子数$m$：一般取值为$20\sim 40$。
2.惯性因子$\omega$：若$\omega$为变量应在迭代开始时设大，在迭代过程中逐步变小。若为定值一般取$0.6\sim 0.75$。
3.加速常数$c_{1}、c_{2}$：一般取$c_{1}=c_{2}=2.0$
4.最大飞翔速度$v_{max}$：一般取每维变化范围的$10％\sim 20％$。

例题及相关函数模版

function main()
clc;clear;close all;
 
tic;    %程序运行计时
 
%设置相关参数值
E0 = 0.001;             %允许误差
MaxNum = 100;           %粒子最大迭代次数
narvs = 1;              %目标函数的自变量个数
particlesize = 30;      %粒子群规模
c1 = 2;                 %个体学习因子，加速常数1
c2 = 2;                 %社会学习因子，加速常数2
w = 0.6;                %惯性因子
vmax = 0.8;             %粒子最大飞行速度
 
%初始化粒子的速度和位置以及目标函数
v = 2*rand(particlesize,narvs);                         %粒子的飞行速度
x = -5 + 10*rand(particlesize,narvs);                   %粒子所在的位置
fitness = @(x) 1/(1+(2.1*(1-x+2*x.^2).*exp(-x.^2/2)));  %适应度函数
%用匿名函数定义适应度函数（速度比单独编写一个m文件要快）
%或者用inline定义  inline('1/(1+(2.1*(1-x+2*x.^2).*exp(-x.^2/2)))','x');
 
%初始化粒子群的个体和全局最大值
f(particlesize,1) = 0;
for i=1:particlesize
    for j=1:narvs
        f(i,j) = fitness(x(i,j));
    end
end
personalbest_x = x;
personalbest_fval = f;
[globalbest_fval,i] = min(personalbest_fval);
globalbest_x = personalbest_x(i,:);
 
%迭代运算
k = 1;
while k <= MaxNum
    %更新粒子的速度和位置信息
    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
    %更新粒子群的个体和全局最大值
    for i = 1:particlesize
        for j = 1:narvs
            f(i) = fitness(x(i,j));
        end
        if f(i) < personalbest_fval
            personalbest_x(i,:) = x(i,:);
            personalbest_fval(i) = f(i);
        end
    end
    [globalbest_fval, i] = min(personalbest_fval);
    globalbest_x = personalbest_x(i,:);
    %检查是否符合停止迭代的条件
    if abs(globalbest_fval) < E0,break,end
    k = k+1;
end
 
%输出最大值结果
Value1 = 1/globalbest_fval-1;Value1 = num2str(Value1);
%strvat指令可以实现字符的组合输出
disp(strcat('the maximun 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_fval-1,'kp','linewidth',4);
legend('目标函数','搜索到的最大值');
xlabel('x');ylabel('y');
grid on;
 
toc;    %程序运行计时结束

实机运行结果