1.背景
1995年,James Kennedy和Russell Eberhart受到鸟群觅食行为的规律性启发,提出了粒子群优化算法(Particle Swarm Optimization, PSO) 。
2.算法原理
2.1算法思想
粒子由速度和位置影响,位置代表解空间中一个解,速度代表下一次位置更新的方向和距离。
粒子群算法基于群体信息共享,假设每个粒子第
i
+
1
i+1
i+1次运动方向受第
i
i
i次个体惯性运动方向,第
i
i
i次个体最优方向和第
i
i
i次群体最优方向影响。
2.2算法过程
群体位置和速度初始化:
X
=
l
b
+
r
a
n
d
∗
(
u
b
−
l
b
)
V
=
V
m
i
n
+
r
a
n
d
∗
(
V
m
a
x
−
V
m
i
n
)
X=lb+rand*(ub-lb)\\ V=Vmin+rand*(Vmax-Vmin)
X=lb+rand∗(ub−lb)V=Vmin+rand∗(Vmax−Vmin)
其中,
u
b
,
l
b
ub,lb
ub,lb分别代表粒子上下位置边界,
V
m
a
x
,
V
m
i
n
Vmax,Vmin
Vmax,Vmin分别代表粒子速度边界。
速度更新:由于粒子受个体惯性运动方向,个体最优运动方向和群体最优运动方向影响,速度更新由三者进行矢量叠加。
V
i
+
1
=
w
V
i
+
C
1
(
P
b
e
s
t
i
−
X
i
)
+
C
2
(
G
b
e
s
t
i
−
X
i
)
V^{i+1}=wV^{i}+C_1(P^{i}_{best}-X^{i})+C_2(G^{i}_{best}-X^{i})
Vi+1=wVi+C1(Pbesti−Xi)+C2(Gbesti−Xi)
其中,
w
,
C
1
,
C
2
w,C_1,C_2
w,C1,C2分别代表权重系数和学习因子(实际意义对下一次运动方向的权重)。
位置更新:粒子第
i
+
1
i+1
i+1次位置受第
i
i
i次位置和第
i
+
1
i+1
i+1次速度影响。
X
i
+
1
=
X
i
+
V
i
+
1
X^{i+1}=X^{i}+V^{i+1}
Xi+1=Xi+Vi+1
3.代码实现
% 粒子群算法主函数
function [Best_pos, Best_fitness, Iter_curve, History_pos, History_best] = pso(pop, dim, ub, lb, fobj, vmax, vmin, maxIter)
%input
%pop 种群数量
%dim 问题维数
%ub 变量上边界
%lb 变量下边界
%fobj 适应度函数
%vmax 粒子速度上边界
%vmin 粒子速度下边界
%maxIter 迭代次数
%output
%Best_pos 粒子全局最优位置
%Best_fitness 全局最优位置对应的适应度函数值
%Iter_curve 每次迭代最优适应度值
%Histroy_pos 记录每代粒子群位置
%History_best 记录每代粒子群最优位置
%% 设置学习因子c1,c2
c1 = 2;
c2 = 2;
%% 初始化种群位置和速度
for i=1:dim
X(:,i) = lb(i)+rand(pop,1).*(ub(i) - lb(i));
end
for i=1:dim
V(:,i) = vmin(i)+rand(pop,1).*(vmax(i) - vmin(i));
end
%% 计算适应度
fitness = zeros(1, pop);
for i = 1:pop
fitness(i) = fobj(X(i,:));
end
%% 将初始种群记为历史最优值
pBest = X;
pBestFitness = fitness;
%% 记录初始全局最优值
% 寻找种群中适应度最小的位置
[~, index] = min(fitness);
gBestFitness = fitness(index); %全局最优fitness
gBest = X(index, :); %全局最优粒子位置
% 位置&适应度更新
Xnew = X;
fitnessNew = fitness;
Iter_curve = zeros(1, maxIter);
%% 开始迭代
for t = 1:maxIter
for i = 1:pop
% 对每个粒子速度更新
r1 = rand(1, dim);
r2 = rand(1, dim);
V(i, :) = V(i, :) + c1 * r1 .* (pBest(i, :) - X(i, :)) + c2 * r2 .* (gBest - X(i, :));
% 速度边界检查
V(i, :) = BoundaryCheck(V(i,:), vmax, vmin, dim);
% 位置更新&边界检查
Xnew(i, :) = X(i, :) + V(i, :);
Flag4ub=Xnew(i,:)>ub;
Flag4lb=Xnew(i,:)<lb;
Xnew(i,:)=(Xnew(i,:).*(~(Flag4ub+Flag4lb)))+ub.*Flag4ub+lb.*Flag4lb;
% 计算适应度
fitnessNew(i) = fobj(Xnew(i, :));
% 更新粒子i历史最优值
if fitnessNew(i) < pBestFitness(i)
pBest(i, :) = Xnew(i, :);
pBestFitness(i) = fitnessNew(i);
end
% 更新全局最优值
if fitnessNew(i) < gBestFitness
gBestFitness = fitnessNew(i);
gBest = Xnew(i, :);
end
end
% 将上次迭代的位置和适应度当做本次迭代初始值
X = Xnew;
fitness = fitnessNew;
% 记录最优位置和适应度
Best_pos = gBest;
Best_fitness = gBestFitness;
Iter_curve(t) = gBestFitness;
History_best{t} = Best_pos;
History_pos{t} = X;
end
end
优化问题:
m
i
n
f
(
x
1
,
x
2
)
=
x
1
2
+
x
2
2
min f(x_1,x_2)=x_1^2+x_2^2
minf(x1,x2)=x12+x22
clear,clc,close all
x1 = -10 : .5 : 10;
x2 = x1;
[X1, X2] = meshgrid(x1, x2);
f = X1.^2 + X2.^2;
s = mesh(X1, X2, f, 'FaceAlpha',0.5);
s.FaceColor = 'flat';
% 粒子群参数设定
pop = 50;
dim = 2;
ub = [10, 10];
lb = [-10, -10];
vmax = [2, 2];
vmin = [-2, -2];
maxIter = 100;
fobj = @(x) fun(x);
% pso求解
[Best_pos, Best_fitness, Iter_curve, History_pos, History_best] = pso(pop, dim, ub, lb, fobj, vmax, vmin, maxIter);