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

时间:2024-03-12 12:26:06

注:本篇随笔依据《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;    %程序运行计时结束

实机运行结果