粒子群算法(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的工具箱，不过效果不好，本人水平低劣，又需要的可以联系我。



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

附注：本文为转载文章

出处：http://blog.csdn.net/niuyongjie/article/details/5602104

作者：niuyongjie