粒子群算法(8)---混合粒子群算法的实现

混合粒子群算法将全局粒子群算法与局部粒子群算法结合,其速度更新采用公式:


其中G(k+1)是全局版本的速度更新公式,而L(k+1)是局部版本的速度更新公式,混合粒子群算法采用H(k+1)的公式。
位置更新公式:

                                         

因为是局部版本与全局版本相结合,所以,粒子群的初始化函数应该与局部版本的相同,这里就不列出了,参看粒子群算法(7)中的LocalInitSwarm函数。
关键还是混合粒子群算法的单步更新函数,函数名为HybridStepPso
代码如下:

function [ParSwarm,OptSwarm]=HybridStepPso(ParSwarm,OptSwarm,AdaptFunc,ParticleScope,MaxW,MinW,LoopCount,CurCount)
% 功能描述:混合粒子群算法。将全局版本与局部版本相混合。基本的粒子群算法的单步更新位置,速度的算法
%
%[ParSwarm,OptSwarm]=LocalStepPsoByCircle(ParSwarm,OptSwarm,AdaptFunc,ParticleScope,MaxW,MinW,LoopCount,CurCount)
%
%算法思想:全局版本的速度更新公式:vq(k+1)=w*v(k)+c1*r1*(pg-w)+c2*r2*(pq-w)
%pg为个体历史最优,pq为全局最优
%局部版本的速度更新公式:vl(k+1)=w*v(k)+c1*r1*(pg-w)+c2*r2*(pl-w) pl为邻域最优
%现在速度更新公式vh=n*vq+(1-n)*vl;n属于0到1的一个数
% 输入参数:ParSwarm:粒子群矩阵,包含粒子的位置,速度与当前的目标函数值
%输入参数:OptSwarm:包含粒子群个体最优解与全局最优解的矩阵
%输入参数:ParticleScope:一个粒子在运算中各维的范围;
% 输入参数:AdaptFunc:适应度函数
%输入参数:LoopCount:迭代的总次数
%输入参数:CurCount:当前迭代的次数
%
% 返回值:含意同输入的同名参数
%
%用法:[ParSwarm,OptSwarm]=LocalStepPsoByCircle(ParSwarm,OptSwarm,AdaptFunc,ParticleScope,MaxW,MinW,LoopCount,CurCount)
%
% 异常:首先保证该文件在Matlab的搜索路径中,然后查看相关的提示信息。
%
%编制人:XXX
%编制时间:2010.5.6
%参考文献:XXX
%参考文献:XXX
%
%修改记录
%----------------------------------------------------------------
%2010.5.6
%修改人:XXX
% 添加2*unifrnd(0,1).*SubTract1(row,:)中的unifrnd(0,1)随机数,使性能大为提高
%修改混合因子,从小到大,开始从CurCount/LoopCount开始
%修改C1=2.05,C2=2.05
%修改速度的范围时区间每维范围的0.5(一半)
%参照基于MATLAB的粒子群优化算法程序设计
%
% 总体评价:使用这个版本的调节系数,效果比较好
%

%容错控制
if nargin~=8
    error('输入的参数个数错误。')
end
if nargout~=2
    error('输出的个数太少,不能保证循环迭代。')
end

%开始单步更新的操作

%*********************************************
%***** 更改下面的代码,可以更改惯性因子的变化*****
%---------------------------------------------------------------------
% 线形递减策略
w=MaxW-CurCount*((MaxW-MinW)/LoopCount);
%---------------------------------------------------------------------
%w 固定不变策略
%w=0.7;
%---------------------------------------------------------------------
% 参考文献:陈贵敏,贾建援,韩琪,粒子群优化算法的惯性权值递减策略研究,西安交通大学学报,2006,1
%w 非线形递减,以凹函数递减
%w=(MaxW-MinW)*(CurCount/LoopCount)^2+(MinW-MaxW)*(2*CurCount/LoopCount)+MaxW;
%---------------------------------------------------------------------
%w 非线形递减,以凹函数递减
%w=MinW*(MaxW/MinW)^(1/(1+10*CurCount/LoopCount));
%*****更改上面的代码,可以更改惯性因子的变化*****
%*********************************************
%更改下面代码可以改变混合因子的取值
%-----------------------------------------------
Hybrid=CurCount/LoopCount;
%-----------------------------------------------

% 得到粒子群群体大小以及一个粒子维数的信息
[ParRow,ParCol]=size(ParSwarm);

%得到粒子的维数
ParCol=(ParCol-1)/2;
GlobleSubTract1=OptSwarm(1:ParRow,:)-ParSwarm(:,1:ParCol);%粒子自身历史最优解位置减去粒子当前位置
LocalSubTract1=OptSwarm(1:ParRow,:)-ParSwarm(:,1:ParCol);%粒子自身历史最优解位置减去粒子当前位置
LocalSubTract2=OptSwarm(ParRow+1:end-1,:)-ParSwarm(:,1:ParCol);%粒子邻域最优解位置减去粒子当前位置
%*********************************************
%***** 更改下面的代码,可以更改c1,c2的变化*****
c1=2.05;
c2=2.05;
%---------------------------------------------------------------------
%con=1;
%c1=4-exp(-con*abs(mean(ParSwarm(:,2*ParCol+1))-AdaptFunc(OptSwarm(ParRow+1,:))));
%c2=4-c1;
%----------------------------------------------------------------------
%***** 更改上面的代码,可以更改c1,c2的变化*****
%*********************************************
for row=1:ParRow 
       GlobleSubTract2=OptSwarm(ParRow*2+1,:)-ParSwarm(row,1:ParCol);%全局最优的位置减去每个粒子当前的位置   
       LocalTempV=w.*ParSwarm(row,ParCol+1:2*ParCol)+c1*unifrnd(0,1).*LocalSubTract1(row,:)+c2*unifrnd(0,1).*LocalSubTract2(row,:);
       GlobleTempV=w.*ParSwarm(row,ParCol+1:2*ParCol)+c1*unifrnd(0,1).*GlobleSubTract1(row,:)+c2*unifrnd(0,1).*GlobleSubTract2;
       TempV=Hybrid.*GlobleTempV+(1-Hybrid).*LocalTempV;
   %限制速度的代码
   for h=1:ParCol
       if TempV(:,h)>ParticleScope(h,2)/2.0
           TempV(:,h)=ParticleScope(h,2)/2.0;
       end
       if TempV(:,h)<-ParticleScope(h,2)/2.0
           TempV(:,h)=(-ParticleScope(h,2)+1e-10)/2.0;
           %加1e-10防止适应度函数被零除
       end
   end  
   
   % 更新速度
   ParSwarm(row,ParCol+1:2*ParCol)=TempV;
   
   %*********************************************
   %***** 更改下面的代码,可以更改约束因子的变化*****
   %---------------------------------------------------------------------
   %a=1;
   %---------------------------------------------------------------------
   a=0.729;
   %***** 更改上面的代码,可以更改约束因子的变化*****
   %*********************************************
   
   % 限制位置的范围
   TempPos=ParSwarm(row,1:ParCol)+a*TempV;
   for h=1:ParCol
       if TempPos(:,h)>ParticleScope(h,2)
           TempPos(:,h)=ParticleScope(h,2);
       end
       if TempPos(:,h)<=ParticleScope(h,1)
           TempPos(:,h)=ParticleScope(h,1)+1e-10;           
       end
   end

   %更新位置 
   ParSwarm(row,1:ParCol)=TempPos;
   
   % 计算每个粒子的新的适应度值
   ParSwarm(row,2*ParCol+1)=AdaptFunc(ParSwarm(row,1:ParCol));
   if ParSwarm(row,2*ParCol+1)>AdaptFunc(OptSwarm(row,1:ParCol))
       OptSwarm(row,1:ParCol)=ParSwarm(row,1:ParCol);
   end
end
%for循环结束
%确定邻域的范围
linyurange=fix(ParRow/2);
%确定当前迭代的邻域范围
jiange=ceil(LoopCount/linyurange);
linyu=ceil(CurCount/jiange);
    for row=1:ParRow
        if row-linyu>0&&row+linyu<=ParRow
            tempM =[ParSwarm(row-linyu:row-1,:);ParSwarm(row+1:row+linyu,:)];
            
            [maxValue,linyurow]=max(tempM(:,2*ParCol+1));
            if maxValue>AdaptFunc(OptSwarm(ParRow+row,:))
                OptSwarm(ParRow+row,:)=tempM(linyurow,1:ParCol);
            end
        else
            if row-linyu<=0
                %该行上面的部分突出了边界,下面绝对不会突破边界
                if row==1
                    tempM=[ParSwarm(ParRow+row-linyu:end,:);ParSwarm(row+1:row+linyu,:)];
                    [maxValue,linyurow]=max(tempM(:,2*ParCol+1));
                     if maxValue>AdaptFunc(OptSwarm(ParRow+row,:))
                        OptSwarm(ParRow+row,:)=tempM(linyurow,1:ParCol);
                     end
                else
                    tempM=[ParSwarm(1:row-1,:);ParSwarm(ParRow+row-linyu:end,:);ParSwarm(row+1:row+linyu,:)];
                    [maxValue,linyurow]=max(tempM(:,2*ParCol+1));
                     if maxValue>AdaptFunc(OptSwarm(ParRow+row,:))
                        OptSwarm(ParRow+row,:)=tempM(linyurow,1:ParCol);
                     end
                end
            else
                %该行下面的部分突出了边界,上面绝对不会突破边界
                if row==ParRow
                    tempM=[ParSwarm(ParRow-linyu:row-1,:);ParSwarm(1:linyu,:)];
                    [maxValue,linyurow]=max(tempM(:,2*ParCol+1));
                     if maxValue>AdaptFunc(OptSwarm(ParRow+row,:))
                        OptSwarm(ParRow+row,:)=tempM(linyurow,1:ParCol);
                     end
                else
                    tempM=[ParSwarm(row-linyu:row-1,:);ParSwarm(row+1:end,:);ParSwarm(1:linyu-(ParRow-row),:)];  
                    [maxValue,linyurow]=max(tempM(:,2*ParCol+1));
                     if maxValue>AdaptFunc(OptSwarm(ParRow+row,:))
                        OptSwarm(ParRow+row,:)=tempM(linyurow,1:ParCol);
                     end
                end
            end
        end
    end%for
    %寻找适应度函数值最大的解在矩阵中的位置(行数),进行全局最优的改变 
[maxValue,row]=max(ParSwarm(:,2*ParCol+1));
if AdaptFunc(ParSwarm(row,1:ParCol))>AdaptFunc(OptSwarm(ParRow*2+1,:))
    OptSwarm(ParRow*2+1,:)=ParSwarm(row,1:ParCol);
end
  
注意代码的91行到96行,这几行就是混合粒子群速度更新公式,其他部分基本与前面的实现一样。
最后还是一个把这两个函数组装在一起的函数,同样采用LocalPsoProcessByCircle函数,详细见粒子群算法(7)的内容,最后还是给出一个应用实例。
Scope=[-10 10;-10 10;-10 10;-10 10;-10 10;-10 10;-10 10;-10 10;-10 10;-10 10];//粒子每维的限制范围
qun=20;//粒子群的规模
lizi=10;//每个粒子的维数
[Result,OnLine,OffLine,MinMaxMeanAdapt,BestofStep]=LocalPsoProcessByCircle(qun,lizi,Scope,@localinitswarm,@Hybridsteppso,@Rastrigin,0,0,1000,0);

注意:在这个LocalPsoProcessByCircle函数中,使用HybridStepPso作为单步更新的函数,其余基本与局部粒子群算法相同。
经过本人的实际测试,运行条件相同,最好的是局部版本的PSO,混合的PSO并不像有些文献上说的那么好,也许是我实现的不对,如果有那个大侠实现的效果更好,可以给我联系,我们可以共享代码。
同时也希望那些砖家、叫兽们共享你们的效果非常好的代码。

本人已经实现了一个PSO的工具箱,不过效果不好,本人水平低劣,又需要的可以联系我。


不知道优快云能不能做链接下载,如果可以,请告诉我,我做个链接,大家可以随便下载,共同交流。

附注:本文为转载文章

出处:http://blog.youkuaiyun.com/niuyongjie/article/details/5602104

作者:niuyongjie

评论 6
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值