PSO 粒子群优化算法
PSO粒子群优化算法
PSO的基本原理
PSO算法把优化问题的解抽象成了粒子。如果把粒子想象成一只鸟，那么从一组初始解出发寻求最优解的过程就类似于鸟群不断调整方向来寻找食物这一过程，其主要包括以下四个部分：
1. 一群鸟在一片区域内寻找食物；
2. 所有的鸟都不知道食物在哪里，但是它们知道当前位置和以往飞行的位置哪个离食物更近，以及通过群体写作知道那只鸟的位置离食物更近；
3. 每只鸟根据自己的飞行经验和群体中距离食物最近的鸟的指引不断地调整自己的飞行方向；
4. 所有的鸟会向离食物最近的鸟靠拢，去搜寻该鸟周围的区域，直至最终找到食物。
用算法语言描述就是，粒子就是鸟，变量空间就是一片区域，问题的可能解对应着每只鸟的位置，问题的最优解就是食物的位置。
数学描述，如图所示：
算法流程：
PSO算法的主要参数的选取分析
- 初始种群规模M：算法从M个初始解出发开始寻优，M越大代表种群包含得到粒子个数越多，算法的寻优能力越强，但相应的计算量越大，种群包含的粒子个数越少，寻优能力越弱，当相应的计算量越小。为了平衡寻优能力与计算量对算法效率的影响，一般将种群规模设为20~50。
- 粒子的维数N：粒子的维数就是粒子的位置向量的维数，粒子的位置向量代表问题的可能解。需要求解的实际问题一旦确定，粒子的维数也就确定了，因此N是一个固定值。
- 粒子的最大速度Vmax：为了避免算法收敛的过快，在实际应用中还需要用户设定一个最大速度限制值Vmax，以将粒子的运动速度控制在[-Vmax,Vmax]。如果Vmax的值过大，粒子可能会飞过最优区域，如果粒子速度过小会陷入局部最优。
- 学习因子：PSO算法的速度更新公式分为三部分：Vi(t-1)是”惯性“部分，表示粒子对上一次迭代速度的继承；C1 * rand() * (Pi(t-1) - Xi(t-1))是”认知“部分，该部分有Pi(t-1)与Xi(t-1)的差值来产生速度分量，因此这部分代表粒子向自身学习的能力，即认知能力；C2 * rand() * (G(t-1) - Xi(t-1))是”社会“部分，该部分由G(t-1)与Xi(t-1)的差值来产生速度分量，而G(t-1)是粒子群体的局部最优解。
PSO算法实际应用
例题：求函数f(x,y) = 3cos(xy)+x+y*y的最小值，其中，x的取值范围[-4.4],y的取值范围[-4,4];
示例代码：
%初始化种群
clear;
clc;
f = @(x,y) 3.* cos(x .* y) + x + y.*y;
x0 = [-4,0.05,4];
y0 = x0;
[X,Y] = meshgrid(x0,y0);
Z = f(X,Y);
figure(1);mesh(X,Y,Z);
colormap(parula(5));
N = 20; %种群大小
D = 2; %粒子维数为2
t = 200; %最大迭代次数
c1 = 2; %自我学习因子
c2 = 2; %群体学习因子
W=0.8; %惯性权重
x_limit = [-4,4];%位置大小范围
v_limit = [-1,1];%速度大小范围
x = x_limit(1)+(x_limit(2)-x_limit(1)).*rand(N,D); %初始化种群位置
v = rand(N,D); %初始化个体速度
xm = x; %每个个体的历史最佳位置
ym = zeros(1,D); %种群历史最佳位置
fxm = ones(N,1)*inf; %每个个体的历史最佳适应度
fym = inf; %种群的历史最佳适应度
hold on;
scatter3( x(:,1),x(:,2) ,f( x(:,1),x(:,2) ),'r*' );
figure(2)
record=[];
% 群体更新
iter = 1;
while iter < t
fx = f(x(:,1),x(:,2)); %个体当前适应度
for i = 1:N
if fx(i)<fxm(i)
fxm(i) = fx(i);%更新个体历史最佳适应度
xm(i,:) = x(i,:);%更新个体历史最佳位置
end
end
if min(fxm)<fym
[fym,nmin] = min(fxm);%更新群体历史最佳适应度
ym = xm(nmin,:);%更新群体历史最佳位置
end
v = vW + c1 * rand (xm-x)+c2rand(repmat(ym,N,1)-x);%速度更新
%边界速度处理
v(v>v_limit(2)) = v_limit(2);
v(v<v_limit(1)) = v_limit(1);
x = x + v;%位置更新
%位置边界处理
x(x>x_limit(2)) = x_limit(2);
x(x<x_limit(1)) = x_limit(1);
record(iter) = fym;%最大值记录
subplot(1,2,1)
mesh(X,Y,Z)
hold on
scatter3( x(:,1),x(:,2) ,f( x(:,1),x(:,2) ) ,'r*');title(['状态位置变化','-迭代次数：',num2str(iter)])
subplot(1,2,2);plot(record);title('最优适应度进化过程')
pause(0.01)
iter = iter+1;
end
figure(4);mesh(X,Y,Z); hold on
scatter3( x(:,1),x(:,2) ,f( x(:,1),x(:,2) ) ,'r*');title('最终状态位置')
disp(['最优值：',num2str(fym)]);
disp(['变量取值：',num2str(ym)]);
浙公网安备 33010602011771号