注:本篇随笔依据《Matlab在数学建模上的应用》中第5章介绍来写,主要介绍粒子群算法思想及其Matlab实现
(博客以及Matlab小白,若有不当欢迎指出)
粒子群算法(PSO)简介
PSO属于智能算法,智能算法都属于软计算(动态自适应的求解方式)。
PSO依托群鸟觅食模型(Boid模型)寻找最优值。
粒子群算法的基本理论
群鸟觅食模型中,每只鸟的飞行基于自身经验和群体经验。
Boid模型遵守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; %程序运行计时结束
实机运行结果