粒子群算法matlab实现（新版）

粒子群算法matlab实现

往期回顾
粒子群算法的基本原理： （初探粒子群优化算法（PSO））
曾经编写的源码：（粒子群算法matlab…）
本期动机
曾经写的源码从功能上来说还是比较完整的。但现在看来，当时我在编码的技巧上还过于幼稚，不会使用张量化的操作或运算。因此本期分享一个更高效的粒子群源码，并逐行注释。

function [gbest,gbestval,record,fitcount]= PSO_func(fhd,Dimension,Particle_Number,Max_Gen,VRmin,VRmax)
%#### input ####%
% fhd:目标函数
% Dimension:决策向量的维度
% Particle_Number:种群数量
% Max_Gen:最大迭代次数
% VRmin:决策向量最小值（1*D）
% VRmax:决策向量最大值（1*D）
%#### output ####%
% gbest:最优位置
% gbestval:最优函数值
% record:每代的全局最优（1*Max_Gen）
% fitcount:最终的评估次数

    rand('state',sum(100*clock)); % 更新随机数流的初试状态，增加随机性。
    me=Max_Gen; % 最大迭代次数
    ps=Particle_Number; % 种群规模（粒子数量）
    D=Dimension; % 维度
    cc=[2 2];   % 加速度参数（就是公式里的C1和C2）
    iwt=0.9-(1:me).*(0.5./me); % 惯性参数
    % iwt=0.5.*ones(1,me);
    
    % 所有决策变量的取值范围都相同时，只需要输入那一个值，而不是所有值重复的一个向量
    if length(VRmin)==1 
        VRmin=repmat(VRmin,1,D);
        VRmax=repmat(VRmax,1,D);
    end
    
    %#### 初始化 ####%
    % 这部分没有注释的行都是为了后面迭代方便而构成的拼接向量，有多少个粒子就拼接多少个向量。
    mv=0.5*(VRmax-VRmin); % 决策变量取值范围的平均值
    VRmin=repmat(VRmin,ps,1);
    VRmax=repmat(VRmax,ps,1);
    Vmin=repmat(-mv,ps,1); % 速度的最小值
    Vmax=-Vmin; % 速度的最大值
    pos=VRmin+(VRmax-VRmin).*rand(ps,D); % 初始化种群位置
    e=feval(fhd,pos'); % 计算适应度
    fitcount=ps; % 统计初始化的评估次数
    vel=Vmin+2.*Vmax.*rand(ps,D); % 初始化粒子的速度
    pbest=pos; % 个体最优适位置
    pbestval=e; % 个体最优适应度
    [gbestval,gbestid]=min(pbestval); % 全局最优适应度和全局最优粒子索引
    gbest=pbest(gbestid,:);% 全局最优位置
    record=repmat(gbestval,1,me); % 记录每代的全局最优适应度
    gbestrep=repmat(gbest,ps,1); 
    
    %#### 迭代过程 ####%
    for i=2:me % 这里从2开始是因为初始化算第一次迭代

        aa=cc(1).*rand(ps,D).*(pbest-pos)+cc(2).*rand(ps,D).*(gbestrep-pos); % 加速度更新
        vel=iwt(i).*vel+aa; % 速度更新
        vel=(vel>Vmax).*Vmax+(vel<=Vmax).*vel; % 控制粒子速度不超过最大速度
        vel=(vel<Vmin).*Vmin+(vel>=Vmin).*vel; % 控制粒子速度不低于最小速度
        pos=pos+vel; % 位置更新
        pos=((pos>=VRmin)&(pos<=VRmax)).*pos...
            +(pos<VRmin).*(VRmin+0.25.*(VRmax-VRmin).*rand(ps,D))+(pos>VRmax).*(VRmax-0.25.*(VRmax-VRmin).*rand(ps,D)); % 控制粒子位置不脱离搜索空间
        e=feval(fhd,pos'); % 计算适应度
        fitcount=fitcount+ps; % 更新评估次数
        tmp=(pbestval<e); % 粒子新的适应度是否最优，最优为1，否则为0
        temp=repmat(tmp',1,D); 
        pbest=temp.*pbest+(1-temp).*pos; % 更新个体最优位置
        pbestval=tmp.*pbestval+(1-tmp).*e; % 更新个体最优适应度
        [gbestval,tmp]=min(pbestval); % 更新全局最优适应度和粒子索引
        gbest=pbest(tmp,:); % 更新全局最优位置
        gbestrep=repmat(gbest,ps,1);
        record(i)=(gbestval>=record(i))*record(i)+(gbestval<record(i))*gbestval; % 记录每代的全局最优适应度
    end
end

代码原型为CEC数值优化竞赛的示例，略加修改并注释。