CLPSO源码逐行注解-CSDN博客

本文链接：https://blog.csdn.net/qq_45731845/article/details/136258012

文章目录

前言
CLPSO源码注解版
总结

前言

	最近在学习CLPSO的时候，发现作者的源码写的可读性非常差，网上也没有相关的注解，所以在这里我花了点时间写了相关的注解在里面，发出来供大家学习，我也可能会有理解不到位的地方，可以讨论交流一下

CLPSO源码注解版

function [gbest,gbestval,fitcount]= CLPSO_new_func(fhd,Max_Gen,Max_FES,Particle_Number,Dimension,VRmin,VRmax,varargin)
%[gbest,gbestval,fitcount]= CLPSO_new_func('f8',3500,200000,30,30,-5.12,5.12)
rand('state',sum(100*clock));
%最大迭代次数
me=Max_Gen;
%种群数
ps=Particle_Number;
%D 是维数的意思
D=Dimension;
%加速度常数
cc=[1 1];   %acceleration constants
%首先在0-1上生成一个间隔为（ps-1)的变量
%再把这个变量映射到0-5上
t=0:1/(ps-1):1;t=5.*t;
%论文中提到的交叉概率 --这样就和论文中的公式对上了
Pc=0.0+(0.5-0.0).*(exp(t)-exp(t(1)))./(exp(t(ps))-exp(t(1)));
% Pc=0.5.*ones(1,ps);
%确定局部搜索中选择的维度数量
%m中的每个元素对应一个粒子表明该粒子在局部搜索中被选择的维度数量
m=0.*ones(ps,1);
%惯性权重
iwt=0.9-(1:me)*(0.7/me);
% iwt=0.729-(1:me)*(0.0/me);
cc=[1.49445 1.49445];
%假如你只给了一个数 就给他扩展一下 方便后续操作
if length(VRmin)==1
    VRmin=repmat(VRmin,1,D);
    VRmax=repmat(VRmax,1,D);
end
%最大位置范围
mv=0.2*(VRmax-VRmin);
%然后再扩展到n个粒子的
VRmin=repmat(VRmin,ps,1);
VRmax=repmat(VRmax,ps,1);
%这个才是后续的限制吗？刚开始很大 然后再把速度限制小了在位置0.2倍以内
Vmin=repmat(-mv,ps,1);
Vmax=-Vmin;
%初始化位置
pos=VRmin+(VRmax-VRmin).*rand(ps,D);

%e(i,1)保存适应度的值 feval---传值估计
for i=1:ps;
    e(i,1)=feval(fhd,pos(i,:),varargin{:});
end
%fitcount成了粒子数量了
fitcount=ps;
%初始化速度
vel=Vmin+2.*Vmax.*rand(ps,D);%initialize the velocity of the particles
%pbset：当前个人最优位置
pbest=pos;
%当前个人最优适应度
pbestval=e; %initialize the pbest and the pbest's fitness value
%gbestval:全局最佳适应度
[gbestval,gbestid]=min(pbestval);
%全局最佳个人
gbest=pbest(gbestid,:);%initialize the gbest and the gbest's fitness value
%将维度进行广播方便后续计算
gbestrep=repmat(gbest,ps,1);
%刷新时间
stay_num=zeros(ps,1);

%这里ai干嘛不太清楚
ai=zeros(ps,D);
%生成了一个矩阵是40行5列 每一列是1-40
%这个初始化有意思哎
%它的里面可以打印出来其实就是比如第一个粒子 第一行就是 1 1 1 1 ，第二个就是
%2 2 2 2..就是对应论文里面说的当PC和生成随机数相比更小的时候，要保留原来的最佳
f_pbest=1:ps;f_pbest=repmat(f_pbest',1,D);
for k=1:ps
    ar=randperm(D);
    %ai 矩阵的作用是指导每个粒子在搜索过程中选择局部搜索和全局搜索的维度
    %在循环中，对于每个粒子 k，通过随机置换 D 维的顺序（ar=randperm(D)），
    %然后根据 m(k) 的值，将前 m(k) 个维度标记为1，表示采用全局搜索策略，其余维度标记为0，表示采用局部搜索策略
    %在代码中，m 矩阵的初始值是全0，因此所有粒子都会选择局部搜索策略
    ai(k,ar(1:m(k)))=1;      %%%大家仔细看这里，看似作者希望加上全局最优项，实际上m并未更新，全局最优也一直没有使用
    fi1=ceil(ps*rand(1,D));
    fi2=ceil(ps*rand(1,D));
    fi=(pbestval(fi1)<pbestval(fi2))'.*fi1+(pbestval(fi1)>=pbestval(fi2))'.*fi2;
    %量用于控制局部最优解和全局最优解之间的权衡
    %这个rand(1,D)是生成（0-1)之间的数一共有D个
    bi=ceil(rand(1,D)-1+Pc(k));
    %如果bi == 0 则随机保证至少一个维度保证有被选中的
    %所以bi肯定在0-1之间
    if bi==zeros(1,D),rc=randperm(D);bi(rc(1))=1;end
    %更新每个粒子的最优解。 bi等于0 说明相信原来初始化的
    %bi等于1 说明要的是目前生成的
    f_pbest(k,:)=bi.*fi+(1-bi).*f_pbest(k,:);
end

stop_num=0;
i=1;


while i<=me&fitcount<=Max_FES
    i=i+1;
    for k=1:ps
        %如果某个粒子的停留次数超过了阈值
        if stay_num(k)>=5
            %     if round(i/10)==i/10%|stay_num(k)>=5
            stay_num(k)=0;
            %重置ai矩阵（局部搜索矩阵）
            ai(k,:)=zeros(1,D);
            %重置这个粒子的综合学习矩阵位置？
            f_pbest(k,:)=k.*ones(1,D);
            ar=randperm(D);
            %选择有多少个给他置1进行全局搜索---实际上论文里没用到
            ai(k,ar(1:m(k)))=1;
            fi1=ceil(ps*rand(1,D));
            fi2=ceil(ps*rand(1,D));
            %随便选俩粒子看谁牛逼
            fi=(pbestval(fi1)<pbestval(fi2))'.*fi1+(pbestval(fi1)>=pbestval(fi2))'.*fi2;
            %bi对应论文里的选择问题---大于了选择CL竞争得到的 小于了选自己本身维度
            bi=ceil(rand(1,D)-1+Pc(k));
            if bi==zeros(1,D),rc=randperm(D);bi(rc(1))=1;end
            %一部分相信被人 一部分相信其他（其他因为他初始化的时就是记录的本维度）
            f_pbest(k,:)=bi.*fi+(1-bi).*f_pbest(k,:);
        end
        %所以最后的f_best就是一堆下标的合集 比如第一行就是 1 1 5 8..
        %表明第一个粒子的第一维向第一个粒子学习 第三维向第5个粒子的第三维学习
        for dimcnt=1:D
            %pbest(f_pbest(k,dimcnt),dimcnt)
            %--这个是根据p_bset的记录的（综合学习得到的需要从哪些粒子学习的下标）去更新个人最佳位置pbset：历史个人最佳的位置
            %得到一个新位置 pbest_f() 记录为新的个人最佳位置
            pbest_f(k,dimcnt)=pbest(f_pbest(k,dimcnt),dimcnt);
        end
        %这里的话因为ai为0，所以cc(2).*ai(k,:).*rand(1,D).*(gbestrep(k,:)-pos(k,:))
        %是没用的 和论文就对上了 没有使用全局最优位置
        aa(k,:)=cc(1).*(1-ai(k,:)).*rand(1,D).*(pbest_f(k,:)-pos(k,:))+cc(2).*ai(k,:).*rand(1,D).*(gbestrep(k,:)-pos(k,:));%~~~~~~~~~~~~~~~~~~~~~~
        %论文中提到的公式
        vel(k,:)=iwt(i).*vel(k,:)+aa(k,:);
        %更新速度 在上下界内
        vel(k,:)=(vel(k,:)>mv).*mv+(vel(k,:)<=mv).*vel(k,:);
        vel(k,:)=(vel(k,:)<(-mv)).*(-mv)+(vel(k,:)>=(-mv)).*vel(k,:);
        pos(k,:)=pos(k,:)+vel(k,:);
        
        %这一步的意思是，我只给那些位置满足要求的人计算适应度
        if (sum(pos(k,:)>VRmax(k,:))+sum(pos(k,:)<VRmin(k,:)))==0;
            e(k,1)=feval(fhd,pos(k,:),varargin{:});
            fitcount=fitcount+1;
            %更新完适应度发现不对 变小了--就要记仇了
            tmp=(pbestval(k)<=e(k));
            if tmp==1
                stay_num(k)=stay_num(k)+1;
            end
            %temp其实就是由全是 1和全是0组成的
            temp=repmat(tmp,1,D);
            %这一步其实可以看做是一个if
            %temp都是1 那就说明这次更新不行 个人最佳是原来的
            %否则个人最佳就是最新的
            pbest(k,:)=temp.*pbest(k,:)+(1-temp).*pos(k,:);
            %更新适应度的值同理
            pbestval(k)=tmp.*pbestval(k)+(1-tmp).*e(k);%update the pbest
            %更新全局最优-----这个虽然这里没用上 后面找最佳还是有用的
            if pbestval(k)<gbestval
                gbest=pbest(k,:);
                gbestval=pbestval(k);
                gbestrep=repmat(gbest,ps,1);%update the gbest
            end
        end
        %假如这次更新没满足个人要求 那个人最佳位置和适应度就不会变化咯        
    end
    
    % if round(i/100)==i/100
    %     plot(pos(:,D-1),pos(:,D),'b*');hold on;
    %     for k=1:floor(D/2)
    %         plot(gbest(:,2*k-1),gbest(:,2*k),'r*');
    %     end
    %     hold off
    %     title(['PSO: ',num2str(i),' generations, Gbestval=',num2str(gbestval)]);
    %     axis([VRmin(1,D-1),VRmax(1,D-1),VRmin(1,D),VRmax(1,D)])
    %     drawnow
    % end
    
    if fitcount>=Max_FES
        break;
    end
    if (i==me)&(fitcount<Max_FES)
        i=i-1;
    end
end
gbestval