PSO算法原理

粒子群优化( Particle Swarm Optimization，PSO) 算法是 Kennedy 和 Eberhart〔1〕在 1995 年首次提出一种模拟鸟群和鱼群等动物寻找食物的社会行为的群智能优化算法。
在粒子群算法中，每个粒子拥有两个特征：速度和位置，如下：
在这里插入图片描述

进化操作

粒子的位置变化公式：
在这里插入图片描述
其中：
c1 和 c2是学习因子，本文定义c1=1.6，c2=1.8

第一个公式是速度更新的公式：
该公式的第一项是惯性部分、第二项是认知部分、第三项是社会部分。
认知部分即：向着自身在历史中所处的最好位置的方向移动。
社会部分即：向着群体中最好的位置的方向移动。

为了改善PSO的性能，Shi和Eberhart引入了惯性权重omiga。速度更新公式即为：
在这里插入图片描述
引入omiga的目的是为了较好地平衡PSO算法的全局搜索能力和局部搜索能力，这样的PSO算法也被称为标准PSO算法。
本文定义omiga的取值范围为
①：权重向量线性递减，Shi和Eberhart将omiga的取值范围定义为[0.4,0.9]

②：也有些研究将权重向量线性递增，公式如下，本文验证部分测试集的时候效果不如权重向量线性递减策略。

在这里插入图片描述

算法流程图

在这里插入图片描述

matlab代码实现

main函数部分

clear all;
close all;
clc;
N=300;
Dim=30;
ub=100;
lb=-100;
T=500;
% 学习因子
c1=1.6;
c2=1.8;
% 速度系数
vel=0.01;
[P_gbest,G_best]=PSO(Dim,ub,lb,vel,N,T,c1,c2);
figure,
plot(G_best,'Color','red');
xlim([1,500]);

适应度函数部分

function my_fitness=my_function(X)
%F1测试问题[-100,100] min=0
% o=sum(X.^2);
% my_fitness=o;

%F2测试问题[-10,10] min=0;
% o=sum(abs(X))+prod(abs(X));
% my_fitness=o;


%F3测试问题——结果极差 min=0
% dim=size(X,2);
% o=0;
% for i=1:dim
%     o=o+sum(X(1:i))^2;
% end
% my_fitness=o;

%F4测试问题——结果极差[-100,100] min=0
% o=max(abs(X));
% my_fitness=o;

%F5测试问题[-30,30] min=0
% dim=size(X,2);
% my_fitness=sum(100*(X(2:dim)-(X(1:dim-1).^2)).^2+(X(1:dim-1)-1).^2);

%F6测试问题[-100,100] min=0
o=sum(abs((X+.5)).^2);
my_fitness=o;
end

PSO算法主体

function [P_gbest,G_best]=PSO(Dim,ub,lb,vel,N,T,c1,c2)
%% 初始化种群
Pop=zeros(N,Dim);
Pop_vel=zeros(N,Dim);
for i=1:N
    for j=1:Dim
        Pop(i,j)=lb+rand().*(ub-lb);%初始化个体
        Pop_vel(i,j)=vel*(lb+rand().*(ub-lb));%初始化个体速度
    end
end
% 计算个体的适应度值
for i=1:N
    fitness(i)=my_function(Pop(i,:));      
end
[fit_gbest,fit_gbest_index]=min(fitness);
P_gbest=Pop(fit_gbest_index,:);%全局中最好的个体
% 初始时 每个个体当前位置即历史中最好的位置
P_pbest=Pop;
%% 进化部分
for t=1:T
    %omiga=2-(2-0.4)*(t/T);%定义范围为[0.4,2],随着迭代次数线性减少
    omiga=(0.9-0.4)*(1-(t/T))+0.4;%惯性权重线性递减策略
    %omiga=(0.4-0.9)*(1-(t/T))+0.9;%惯性权重递增策略——效果不佳
    % 遍历所有个体，更新所有个体的位置
    for i=1:N
        for j=1:Dim
            r1=rand();
            r2=rand();
            Pop_vel(i,j)=omiga*Pop_vel(i,j)+c1*r1*(P_pbest(i,j)-Pop(i,j))+c2*r2*(P_gbest(j)-Pop(i,j));
            Pop(i,j)=Pop(i,j)+Pop_vel(i,j);
        end
    end
    % 遍历个体计算适应度值
    for i=1:N
        Flag4Up=Pop(i,:)>ub;
        Flag4Lp=Pop(i,:)<lb;
        Pop(i,:)=(Pop(i,:).*(~(Flag4Up+Flag4Lp)))+Flag4Up.*ub+Flag4Lp.*lb;%修正粒子群的位置

        FlagUp=Pop_vel(i,:)>vel*ub;
        FlagLp=Pop_vel(i,:)<vel*lb;
        Pop_vel(i,:)=(Pop_vel(i,:).*(~(FlagUp+FlagLp)))+FlagUp.*vel.*ub+FlagLp.*lb.*vel;%修正粒子群的速度

        y=my_function(Pop(i,:));
        if y<fitness(i)
            P_pbest(i,:)=Pop(i,:);%更新在历史中每个粒子最好的位置
            fitness(i)=y;
        end
    end
    % 更新最好个体的位置
    [fit_best,fit_best_index]=min(fitness);
    G_best(t)=fit_best;
    P_gbest=Pop(fit_best_index,:);
end
end