【智能算法】粒子群算法（PSO）原理及实现

原创

小O的算法实验室

修改于 2024-03-05 13:52:04

6780

修改于 2024-03-05 13:52:04

文章被收录于专栏：智能算法

1.背景

1995年，James Kennedy和Russell Eberhart受到鸟群觅食行为的规律性启发，提出了粒子群优化算法(Particle Swarm Optimization, PSO) 。

2.算法原理

2.1算法思想

粒子由速度和位置影响，位置代表解空间中一个解，速度代表下一次位置更新的方向和距离。

粒子群算法基于群体信息共享，假设每个粒子第i+1 次运动方向受第i 次个体惯性运动方向，第i 次个体最优方向和第i 次群体最优方向影响。

2.2算法过程

群体位置和速度初始化:

X=lb+rand*(ub-lb)\\V=Vmin+rand*(Vmax-Vmin)

其中，ub,lb 分别代表粒子上下位置边界，Vmax,Vmin 分别代表粒子速度边界。

速度更新：由于粒子受个体惯性运动方向，个体最优运动方向和群体最优运动方向影响，速度更新由三者进行矢量叠加。

V^{i+1}=wV^{i}+C_1(P^{i}\_{best}-X^{i})+C_2(G^{i}\_{best}-X^{i})

其中，w,C_1,C_2 分别代表权重系数和学习因子（实际意义对下一次运动方向的权重）。

位置更新：粒子第i+1 次位置受第i 次位置和第i+1 次速度影响。

X^{i+1}=X^{i}+V^{i+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

优化问题:

min f(x_1,x_2)=x_1^2+x_2^2

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);