混合粒子群算法的实现

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

 

其中G（k+1）是全局版本的速度更新公式，而L（k+1）是局部版本的速度更新公式，混合粒子群算法采用H（k+1）的公式。

位置更新公式

 

因为是局部版本与全局版本相结合，所以，粒子群的初始化函数应该与局部版本的相同，这里就不列出了，参看粒子群算法(7）中的LocalInitSwarm函数。

关键还是混合粒子群算法的单步更新函数，函数名为HybridStepPso

代码如下：

view plaincopy to clipboardprint?
01.function [ParSwarm,OptSwarm]=HybridStepPso(ParSwarm,OptSwarm,AdaptFunc,ParticleScope,MaxW,MinW,LoopCount,CurCount)  
02.% 功能描述：混合粒子群算法。将全局版本与局部版本相混合。基本的粒子群算法的单步更新位置,速度的算法  
03.%  
04.%[ParSwarm,OptSwarm]=LocalStepPsoByCircle(ParSwarm,OptSwarm,AdaptFunc,ParticleScope,MaxW,MinW,LoopCount,CurCount)  
05.%  
06.%算法思想：全局版本的速度更新公式：vq(k+1)=w*v(k)+c1*r1*(pg-w)+c2*r2*(pq-w)  
07.%pg为个体历史最优，pq为全局最优  
08.%局部版本的速度更新公式：vl(k+1)=w*v(k)+c1*r1*(pg-w)+c2*r2*(pl-w) pl为邻域最优  
09.%现在速度更新公式vh=n*vq+(1-n)*vl;n属于0到1的一个数  
10.% 输入参数：ParSwarm:粒子群矩阵，包含粒子的位置，速度与当前的目标函数值  
11.%输入参数：OptSwarm：包含粒子群个体最优解与全局最优解的矩阵  
12.%输入参数：ParticleScope:一个粒子在运算中各维的范围；  
13.% 输入参数：AdaptFunc：适应度函数  
14.%输入参数：LoopCount：迭代的总次数  
15.%输入参数：CurCount：当前迭代的次数  
16.%  
17.% 返回值：含意同输入的同名参数  
18.%  
19.%用法：[ParSwarm,OptSwarm]=LocalStepPsoByCircle(ParSwarm,OptSwarm,AdaptFunc,ParticleScope,MaxW,MinW,LoopCount,CurCount)  
20.%  
21.% 异常：首先保证该文件在Matlab的搜索路径中，然后查看相关的提示信息。  
22.%  
23.%编制人：XXX  
24.%编制时间：2010.5.6  
25.%参考文献：XXX  
26.%参考文献：XXX  
27.%  
28.%修改记录  
29.%----------------------------------------------------------------  
30.%2010.5.6  
31.%修改人：XXX  
32.% 添加2*unifrnd(0,1).*SubTract1(row,:)中的unifrnd(0,1)随机数，使性能大为提高  
33.%修改混合因子，从小到大，开始从CurCount/LoopCount开始  
34.%修改C1=2.05,C2=2.05  
35.%修改速度的范围时区间每维范围的0.5(一半)  
36.%参照基于MATLAB的粒子群优化算法程序设计  
37.%  
38.% 总体评价：使用这个版本的调节系数，效果比较好  
39.%  
40. 
41.%容错控制  
42.if nargin~=8  
43.    error('输入的参数个数错误。')  
44.end  
45.if nargout~=2  
46.    error('输出的个数太少，不能保证循环迭代。')  
47.end  
48. 
49.%开始单步更新的操作  
50. 
51.%*********************************************  
52.%***** 更改下面的代码，可以更改惯性因子的变化*****  
53.%---------------------------------------------------------------------  
54.% 线形递减策略  
55.w=MaxW-CurCount*((MaxW-MinW)/LoopCount);  
56.%---------------------------------------------------------------------  
57.%w 固定不变策略  
58.%w=0.7;  
59.%---------------------------------------------------------------------  
60.% 参考文献：陈贵敏，贾建援，韩琪，粒子群优化算法的惯性权值递减策略研究，西安交通大学学报，2006，1  
61.%w 非线形递减，以凹函数递减  
62.%w=(MaxW-MinW)*(CurCount/LoopCount)^2+(MinW-MaxW)*(2*CurCount/LoopCount)+MaxW;  
63.%---------------------------------------------------------------------  
64.%w 非线形递减，以凹函数递减  
65.%w=MinW*(MaxW/MinW)^(1/(1+10*CurCount/LoopCount));  
66.%*****更改上面的代码，可以更改惯性因子的变化*****  
67.%*********************************************  
68.%更改下面代码可以改变混合因子的取值  
69.%-----------------------------------------------  
70.Hybrid=CurCount/LoopCount;  
71.%-----------------------------------------------  
72. 
73.% 得到粒子群群体大小以及一个粒子维数的信息  
74.[ParRow,ParCol]=size(ParSwarm);  
75. 
76.%得到粒子的维数  
77.ParCol=(ParCol-1)/2;  
78.GlobleSubTract1=OptSwarm(1:ParRow,:)-ParSwarm(:,1:ParCol);%粒子自身历史最优解位置减去粒子当前位置  
79.LocalSubTract1=OptSwarm(1:ParRow,:)-ParSwarm(:,1:ParCol);%粒子自身历史最优解位置减去粒子当前位置  
80.LocalSubTract2=OptSwarm(ParRow+1:end-1,:)-ParSwarm(:,1:ParCol);%粒子邻域最优解位置减去粒子当前位置  
81.%*********************************************  
82.%***** 更改下面的代码，可以更改c1,c2的变化*****  
83.c1=2.05;  
84.c2=2.05;  
85.%---------------------------------------------------------------------  
86.%con=1;  
87.%c1=4-exp(-con*abs(mean(ParSwarm(:,2*ParCol+1))-AdaptFunc(OptSwarm(ParRow+1,:))));  
88.%c2=4-c1;  
89.%----------------------------------------------------------------------  
90.%***** 更改上面的代码，可以更改c1,c2的变化*****  
91.%*********************************************  
92.for row=1:ParRow   
93.       GlobleSubTract2=OptSwarm(ParRow*2+1,:)-ParSwarm(row,1:ParCol);%全局最优的位置减去每个粒子当前的位置     
94.       LocalTempV=w.*ParSwarm(row,ParCol+1:2*ParCol)+c1*unifrnd(0,1).*LocalSubTract1(row,:)+c2*unifrnd(0,1).*LocalSubTract2(row,:);  
95.       GlobleTempV=w.*ParSwarm(row,ParCol+1:2*ParCol)+c1*unifrnd(0,1).*GlobleSubTract1(row,:)+c2*unifrnd(0,1).*GlobleSubTract2;  
96.       TempV=Hybrid.*GlobleTempV+(1-Hybrid).*LocalTempV;  
97.   %限制速度的代码  
98.   for h=1:ParCol  
99.       if TempV(:,h)>ParticleScope(h,2)/2.0  
100.           TempV(:,h)=ParticleScope(h,2)/2.0;  
101.       end  
102.       if TempV(:,h)<-ParticleScope(h,2)/2.0  
103.           TempV(:,h)=(-ParticleScope(h,2)+1e-10)/2.0;  
104.           %加1e-10防止适应度函数被零除  
105.       end  
106.   end    
107.     
108.   % 更新速度  
109.   ParSwarm(row,ParCol+1:2*ParCol)=TempV;  
110.     
111.   %*********************************************  
112.   %***** 更改下面的代码，可以更改约束因子的变化*****  
113.   %---------------------------------------------------------------------  
114.   %a=1;  
115.   %---------------------------------------------------------------------  
116.   a=0.729;  
117.   %***** 更改上面的代码，可以更改约束因子的变化*****  
118.   %*********************************************  
119.     
120.   % 限制位置的范围  
121.   TempPos=ParSwarm(row,1:ParCol)+a*TempV;  
122.   for h=1:ParCol  
123.       if TempPos(:,h)>ParticleScope(h,2)  
124.           TempPos(:,h)=ParticleScope(h,2);  
125.       end  
126.       if TempPos(:,h)<=ParticleScope(h,1)  
127.           TempPos(:,h)=ParticleScope(h,1)+1e-10;             
128.       end  
129.   end  
130. 
131.   %更新位置   
132.   ParSwarm(row,1:ParCol)=TempPos;  
133.     
134.   % 计算每个粒子的新的适应度值  
135.   ParSwarm(row,2*ParCol+1)=AdaptFunc(ParSwarm(row,1:ParCol));  
136.   if ParSwarm(row,2*ParCol+1)>AdaptFunc(OptSwarm(row,1:ParCol))  
137.       OptSwarm(row,1:ParCol)=ParSwarm(row,1:ParCol);  
138.   end  
139.end  
140.%for循环结束  
141.%确定邻域的范围  
142.linyurange=fix(ParRow/2);  
143.%确定当前迭代的邻域范围  
144.jiange=ceil(LoopCount/linyurange);  
145.linyu=ceil(CurCount/jiange);  
146.    for row=1:ParRow  
147.        if row-linyu>0&&row+linyu<=ParRow  
148.            tempM =[ParSwarm(row-linyu:row-1,:);ParSwarm(row+1:row+linyu,:)];  
149.              
150.            [maxValue,linyurow]=max(tempM(:,2*ParCol+1));  
151.            if maxValue>AdaptFunc(OptSwarm(ParRow+row,:))  
152.                OptSwarm(ParRow+row,:)=tempM(linyurow,1:ParCol);  
153.            end  
154.        else 
155.            if row-linyu<=0  
156.                %该行上面的部分突出了边界，下面绝对不会突破边界  
157.                if row==1  
158.                    tempM=[ParSwarm(ParRow+row-linyu:end,:);ParSwarm(row+1:row+linyu,:)];  
159.                    [maxValue,linyurow]=max(tempM(:,2*ParCol+1));  
160.                     if maxValue>AdaptFunc(OptSwarm(ParRow+row,:))  
161.                        OptSwarm(ParRow+row,:)=tempM(linyurow,1:ParCol);  
162.                     end  
163.                else 
164.                    tempM=[ParSwarm(1:row-1,:);ParSwarm(ParRow+row-linyu:end,:);ParSwarm(row+1:row+linyu,:)];  
165.                    [maxValue,linyurow]=max(tempM(:,2*ParCol+1));  
166.                     if maxValue>AdaptFunc(OptSwarm(ParRow+row,:))  
167.                        OptSwarm(ParRow+row,:)=tempM(linyurow,1:ParCol);  
168.                     end  
169.                end  
170.            else 
171.                %该行下面的部分突出了边界，上面绝对不会突破边界  
172.                if row==ParRow  
173.                    tempM=[ParSwarm(ParRow-linyu:row-1,:);ParSwarm(1:linyu,:)];  
174.                    [maxValue,linyurow]=max(tempM(:,2*ParCol+1));  
175.                     if maxValue>AdaptFunc(OptSwarm(ParRow+row,:))  
176.                        OptSwarm(ParRow+row,:)=tempM(linyurow,1:ParCol);  
177.                     end  
178.                else 
179.                    tempM=[ParSwarm(row-linyu:row-1,:);ParSwarm(row+1:end,:);ParSwarm(1:linyu-(ParRow-row),:)];    
180.                    [maxValue,linyurow]=max(tempM(:,2*ParCol+1));  
181.                     if maxValue>AdaptFunc(OptSwarm(ParRow+row,:))  
182.                        OptSwarm(ParRow+row,:)=tempM(linyurow,1:ParCol);  
183.                     end  
184.                end  
185.            end  
186.        end  
187.    end%for 
188.    %寻找适应度函数值最大的解在矩阵中的位置(行数)，进行全局最优的改变   
189.[maxValue,row]=max(ParSwarm(:,2*ParCol+1));  
190.if AdaptFunc(ParSwarm(row,1:ParCol))>AdaptFunc(OptSwarm(ParRow*2+1,:))  
191.    OptSwarm(ParRow*2+1,:)=ParSwarm(row,1:ParCol);  
192.end  
193.   
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）的内容，最后还是给出一个应用实例。

view plaincopy to clipboardprint?
01.Scope=[-10 10;-10 10;-10 10;-10 10;-10 10;-10 10;-10 10;-10 10;-10 10;-10 10];//粒子每维的限制范围  
02.qun=20;//粒子群的规模  
03.lizi=10;//每个粒子的维数  
04.[Result,OnLine,OffLine,MinMaxMeanAdapt,BestofStep]=LocalPsoProcessByCircle(qun,lizi,Scope,@localinitswarm,@Hybridsteppso,@Rastrigin,0,0,1000,0); 
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能不能做链接下载，如果可以，请告诉我，我做个链接，大家可以随便下载，共同交流

本文来自CSDN博客，转载请标明出处：http://blog.csdn.net/niuyongjie/archive/2010/05/17/5602104.aspx