basesteppso.m

粒子群算法的详细实现以及在大肠杆菌代谢过程中参数优化的应用。

function [ParSwarm,OptSwarm]=BaseStepPso(ParSwarm,OptSwarm,AdaptFunc,ParticleScope,MaxW,MinW,LoopCount,CurCount)
%function description: global version: basic update of one step in pso,
%position, velocity
%
%[ParSwarm,OptSwarm]=BaseStepPso(ParSwarm,OptSwarm,AdaptFunc,ParticleScope,MaxW,MinW,LoopCount,CurCount)
%
%inputparameter：ParSwarm:particle matrix，including position,velocity and
%current aim function
%inputparameter：OptSwarm：matrix which includes particle p-best and g-best
%inputparameter：ParticleScope:the scope of the dimension；
%inputparameter：AdaptFunc：fitness function
%inputparameter：LoopCount：iteration time
%inputparameter：CurCount：current iteration time
%
%returen value: same with the one_name inputparameter
%
%usage：[ParSwarm,OptSwarm]=BaseStepPso(ParSwarm,OptSwarm,AdaptFunc,Particle
%Scope,MaxW,MinW,LoopCount,CurCount)
%----------------------------------------------------------------
% add 2*unifrnd(0,1).*SubTract1(row,:) unifrnd(0,1) random number，largely
% improve the performance
%

%fault control
if nargin~=8
    error('fault number of input parameter')
end
if nargout~=2
    error('number of input is too small, can not do iteration')
end

%begin one_time iteratiom

%*********************************************
%*****change below, can do inertia gene change*****
%---------------------------------------------------------------------
%linear descending strategy
%w=MaxW-CurCount*((MaxW-MinW)/LoopCount);
%---------------------------------------------------------------------
%w fixed strategy
w=0.68;
%---------------------------------------------------------------------
%w nonlinear descending strategy
%w=(MaxW-MinW)*(CurCount/LoopCount)^2+(MinW-MaxW)*(2*CurCount/LoopCount)+MaxW;
%---------------------------------------------------------------------
%w nonlinear descending strategy
%w=MinW*(MaxW/MinW)^(1/(1+10*CurCount/LoopCount));
%*****change above, change of inertia gene*****
%*********************************************

%get size of particles and information of dimension for singal particle
[ParRow,ParCol]=size(ParSwarm);
%get the dimension
ParCol=(ParCol-1)/2;
SubTract1=OptSwarm(1:ParRow,:)-ParSwarm(:,1:ParCol);

%*********************************************
%*****change below, change c1,c2*****
c1=2;
c2=2;
%---------------------------------------------------------------------
%con=1;
%c1=4-exp(-con*abs(mean(ParSwarm(:,2*ParCol+1))-AdaptFunc(OptSwarm(ParRow+1,:))));
%c2=4-c1;
%----------------------------------------------------------------------
%*****change above, change c1,c2*****
%*********************************************
for row=1:ParRow
   SubTract2=OptSwarm(ParRow+1,:)-ParSwarm(row,1:ParCol);
   TempV=w.*ParSwarm(row,ParCol+1:2*ParCol)+2*unifrnd(0,1).*SubTract1(row,:)+2*unifrnd(0,1).*SubTract2;
   %velocity restrict code
   for h=1:ParCol
       if TempV(:,h)>ParticleScope(h,2)
           TempV(:,h)=ParticleScope(h,2);
       end
       if TempV(:,h)<-ParticleScope(h,2)
           TempV(:,h)=-ParticleScope(h,2)+1e-10;
           %plus 1e-10 prevent to be divided by 0;
       end
   end  
   
   %update velocity
   ParSwarm(row,ParCol+1:2*ParCol)=TempV;
   
   %*********************************************
   %*****change below, change restrict gene*****
   %---------------------------------------------------------------------
   %a=1;
   %---------------------------------------------------------------------
   a=0.729;
   %*****change above, change restrict gene*****
   %*********************************************
   
   %restrict position scope
   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

   %update position
   ParSwarm(row,1:ParCol)=TempPos;
   
   %calculate new fitness value
   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 loop end

%look for min solution in the matrix position, do g-best change 
%[maxValue,row]=max(ParSwarm(:,2*ParCol+1));
[minValue,row]=min(ParSwarm(:,2*ParCol+1));
%if AdaptFunc(ParSwarm(row,1:ParCol))>AdaptFunc(OptSwarm(ParRow+1,:))
if AdaptFunc(ParSwarm(row,1:ParCol))<AdaptFunc(OptSwarm(ParRow+1,:))
    OptSwarm(ParRow+1,:)=ParSwarm(row,1:ParCol);
end