PSO学习笔记

1.仿生思想

粒子群算法模拟鸟集群飞行觅食的行为,通过鸟之间的集体协作使群体达到目的。相对于其他演化方法,它具有原理简单,易于编程实现等特点。PSO数学描述如下:
初始化一群 随机粒子(随机解)。然后通过迭代找到最优解。在每一次的迭代中,粒子通过跟踪两个“极值”(Pbest,gbest)来更新自己。
在找到这两个最优值后,粒子通过下面的公式来更新自己的速度和位置。
\[\begin{array}{rl}{V_{t, j}(t+1)=w} & {V_{t, j}(t)+c_{1} \times r \text { and }() \times\left(p b e s t_{i, j}-x_{t, j}(t)\right)} \\ {} & {+c_{2} \times r a n d() \times\left(g b e s t_{i . j}-x_{i . j}(t)\right)}\end{array}\] （1）
\[ x_{i, j}(t+1)=x_{i, j}(t)+V_{i, j}(t+1), j=1,2, \ldots, d \]（2）
\[ i=1,2, \ldots, M . M \]
M是该群体中粒子的总数,\[V_{i}\]是粒子的速度;pbest和gbest分别定义为粒子最佳位置和种群最佳位置。rand是介于 (0 , 1 )之间的随机数;\[x_{i}\]是粒子的当前位置。cl和c2是学习因子,通常取c1=c2=1.4962,w为惯性权因子。

基本粒子群算法代码

%% Problem Definition

CostFunction=@(z) Fun(z);        % Cost Function
%Ufun=@(o) Ufun
nVar=4;            % Number of Decision Variables

VarSize=[1 nVar];   % Size of Decision Variables Matrix

VarMin=0;         % Lower Bound of Variables
VarMax=10;         % Upper Bound of Variables


%% PSO Parameters

MaxIt=500;      % Maximum Number of Iterations

nPop=30;        % Population Size (Swarm Size)

% PSO Parameters
w=1;            % Inertia Weight
wdamp=0.99;     % Inertia Weight Damping Ratio
c1=1.5;         % Personal Learning Coefficient
c2=2.0;         % Global Learning Coefficient

% If you would like to use Constriction Coefficients for PSO,
% uncomment the following block and comment the above set of parameters.

% % Constriction Coefficients
% phi1=2.05;
% phi2=2.05;
% phi=phi1+phi2;
% chi=2/(phi-2+sqrt(phi^2-4*phi));
% w=chi;          % Inertia Weight
% wdamp=1;        % Inertia Weight Damping Ratio
% c1=chi*phi1;    % Personal Learning Coefficient
% c2=chi*phi2;    % Global Learning Coefficient

% Velocity Limits
VelMax=0.1*(VarMax-VarMin);
VelMin=-VelMax;

%% Initialization

empty_particle.Position=[];
empty_particle.Cost=[];
empty_particle.Velocity=[];
empty_particle.Best.Position=[];
empty_particle.Best.Cost=[];

particle=repmat(empty_particle,nPop,1);

GlobalBest.Cost=inf;

for i=1:nPop
    
    % Initialize Position
    particle(i).Position=unifrnd(VarMin,VarMax,VarSize);
    
    % Initialize Velocity
    particle(i).Velocity=zeros(VarSize);
    
    % Evaluation
    particle(i).Cost=CostFunction(particle(i).Position);
    
    % Update Personal Best
    particle(i).Best.Position=particle(i).Position;
    particle(i).Best.Cost=particle(i).Cost;
    
    % Update Global Best
    if particle(i).Best.Cost<GlobalBest.Cost
        
        GlobalBest=particle(i).Best;
        
    end
    
end

BestCost=zeros(MaxIt,1);

%% PSO Main Loop

for it=1:MaxIt
    
    for i=1:nPop
        
        % Update Velocity
        particle(i).Velocity = w*particle(i).Velocity ...
            +c1*rand(VarSize).*(particle(i).Best.Position-particle(i).Position) ...
            +c2*rand(VarSize).*(GlobalBest.Position-particle(i).Position);
        
        % Apply Velocity Limits
        particle(i).Velocity = max(particle(i).Velocity,VelMin);
        particle(i).Velocity = min(particle(i).Velocity,VelMax);
        
        % Update Position
        particle(i).Position = particle(i).Position + particle(i).Velocity;
        
        % Velocity Mirror Effect
        IsOutside=(particle(i).Position<VarMin | particle(i).Position>VarMax);
        particle(i).Velocity(IsOutside)=-particle(i).Velocity(IsOutside);
        
        % Apply Position Limits
        particle(i).Position = max(particle(i).Position,VarMin);
        particle(i).Position = min(particle(i).Position,VarMax);
        
        % Evaluation
        particle(i).Cost = CostFunction(particle(i).Position);
        
        % Update Personal Best
        if particle(i).Cost<particle(i).Best.Cost
            
            particle(i).Best.Position=particle(i).Position;
            particle(i).Best.Cost=particle(i).Cost;
            
            % Update Global Best
            if particle(i).Best.Cost<GlobalBest.Cost
                
                GlobalBest=particle(i).Best;
                
            end
            
        end
        
    end
    
    BestCost(it)=GlobalBest.Cost;
    
    disp(['Iteration ' num2str(it) ': Best Cost = ' num2str(BestCost(it))]);
    
    w=w*wdamp;

end

BestSol = GlobalBest;
cg_curve = BestCost;
%% Results

%Draw objective space
subplot(1,2,2);
semilogy((1:50:MaxIt),cg_curve(1:50:MaxIt),'k-d','markersize',10,'Color','g','LineWidth',1.5)
title('Convergence curve','FontSize',16)
xlabel('Iteration','FontSize',16);
ylabel('Best score obtained so far','FontSize',16);
set(gca,'FontSize',16);

hold on
axis tight
grid off
box on
legend('PSO')
set(legend,'FontSize',12);

下面这个代码是基准测试函数，用来检验PSO的收敛性能的

function z=Fun(u)

% F1 Sphere [-100,100]
%z=sum(u.^2);

%F2 Schwel [-10,10]
 %z=sum(abs(u))+prod(abs(u));

%F3 [-100,100]
%dim=size(u,2);
%z=0;
%for i=1:dim
  %z=z+sum(u(1:i)-1)^2;
 % end


% F4 [-100,100]
%z=max(abs(u));

%F5 Rosenbrock  [-30,30]
%dim = 30;
%z=sum(100*(u(2:dim)-(u(1:dim-1).^2)).^2+(u(1:dim-1)-1).^2);


%F6 [-100,100]
%z=sum(abs((u+.5)).^2);

%F7 [-1.28,1.28]
%dim=size(u,2);
%z=sum([1:dim].*(u.^4))+rand;

%F8 [-500,500]
%z=sum(-u.*sin(sqrt(abs(u))));

%F9 Rastrign [-5.12,5.12]
%dim=size(u,2);
%z=sum(u.^2-10*cos(2*pi.*u))+10*dim;

%F10 Ackley [-32,32]
%dim=size(u,2);
%z=-20*exp(-.2*sqrt(sum(u.^2)/dim))-exp(sum(cos(2*pi.*u))/dim)+20+exp(1);

%F11 Griewank [-600,600]
%dim=size(u,2);
%z=sum(u.^2)/4000-prod(cos(u./sqrt([1:dim])))+1;

%F12  %[-50,50]
%dim=size(u,2);
%z=(pi/dim)*(10*((sin(pi*(1+(u(1)+1)/4)))^2)+sum((((u(1:dim-1)+1)./4).^2).*...
%(1+10.*((sin(pi.*(1+(u(2:dim)+1)./4)))).^2))+((u(dim)+1)/4)^2)+sum(Ufun(u,10,100,4));
%function o=Ufun(u,a,k,m)
%o=k.*((u-a).^m).*(u>a)+k.*((-u-a).^m).*(u<(-a));


%F13  %[-50,50]
%dim=size(u,2);
%z=.1*((sin(3*pi*u(1)))^2+sum((u(1:dim-1)-1).^2.*(1+(sin(3.*pi.*u(2:dim))).^2))+...
%((u(dim)-1)^2)*(1+(sin(2*pi*u(dim)))^2))+sum(Ufun(u,5,100,4));
%function o=Ufun(u,a,k,m)
%o=k.*((u-a).^m).*(u>a)+k.*((-u-a).^m).*(u<(-a));

%F14  [-65.536,-65.536]
%aS=[-32 -16 0 16 32 -32 -16 0 16 32 -32 -16 0 16 32 -32 -16 0 16 32 -32 -16 0 16 32;,...
   % -32 -32 -32 -32 -32 -16 -16 -16 -16 -16 0 0 0 0 0 16 16 16 16 16 32 32 32 32 32];
%for j=1:25
  % bS(j)=sum((u'-aS(:,j)).^6);
%end
%z=(1/500+sum(1./([1:25]+bS))).^(-1);

% F15   %[-5,5]
%aK=[.1957 .1947 .1735 .16 .0844 .0627 .0456 .0342 .0323 .0235 .0246];
%bK=[.25 .5 1 2 4 6 8 10 12 14 16];bK=1./bK;
%z=sum((aK-((u(1).*(bK.^2+u(2).*bK))./(bK.^2+u(3).*bK+u(4)))).^2);
%end

%F16    %[-5,5]
  %z=4*(u(1)^2)-2.1*(u(1)^4)+(u(1)^6)/3+u(1)*u(2)-4*(u(2)^2)+4*(u(2)^4);

  % F17   %[-5,5]
%z=(u(2)-(u(1)^2)*5.1/(4*(pi^2))+5/pi*u(1)-6)^2+10*(1-1/(8*pi))*cos(u(1))+10;



% F18    %[-5,5]
%z=(1+(u(1)+u(2)+1)^2*(19-14*u(1)+3*(u(1)^2)-14*u(2)+6*u(1)*u(2)+3*u(2)^2))*...
 % (30+(2*u(1)-3*u(2))^2*(18-32*u(1)+12*(u(1)^2)+48*u(2)-36*u(1)*u(2)+27*(u(2)^2)));


% F19     %[1,3]
%aH=[3 10 30;.1 10 35;3 10 30;.1 10 35];cH=[1 1.2 3 3.2];
%pH=[.3689 .117 .2673;.4699 .4387 .747;.1091 .8732 .5547;.03815 .5743 .8828];
%z=0;
%for i=1:4
    %z=z-cH(i)*exp(-(sum(aH(i,:).*((u-pH(i,:)).^2))));
%end


% F20       %[0,1]
%aH=[10 3 17 3.5 1.7 8;.05 10 17 .1 8 14;3 3.5 1.7 10 17 8;17 8 .05 10 .1 14];
%cH=[1 1.2 3 3.2];
%pH=[.1312 .1696 .5569 .0124 .8283 .5886;.2329 .4135 .8307 .3736 .1004 .9991;...
%.2348 .1415 .3522 .2883 .3047 .6650;.4047 .8828 .8732 .5743 .1091 .0381];
%z=0;
%for i=1:4
 % z=z-cH(i)*exp(-(sum(aH(i,:).*((u-pH(i,:)).^2))));
%end

% F21      %[0,10]
%aSH=[4 4 4 4;1 1 1 1;8 8 8 8;6 6 6 6;3 7 3 7;2 9 2 9;5 5 3 3;8 1 8 1;6 2 6 2;7 3.6 7 3.6];
%cSH=[.1 .2 .2 .4 .4 .6 .3 .7 .5 .5];
%z=0;
%for i=1:5
  %z=z-((u-aSH(i,:))*(u-aSH(i,:))'+cSH(i))^(-1);
%end


%F22      %[0,10]
aSH=[4 4 4 4;1 1 1 1;8 8 8 8;6 6 6 6;3 7 3 7;2 9 2 9;5 5 3 3;8 1 8 1;6 2 6 2;7 3.6 7 3.6];
cSH=[.1 .2 .2 .4 .4 .6 .3 .7 .5 .5];
z=0;
for i=1:7
    z=z-((u-aSH(i,:))*(u-aSH(i,:))'+cSH(i))^(-1);
end


% F23      %[0,10]
%aSH=[4 4 4 4;1 1 1 1;8 8 8 8;6 6 6 6;3 7 3 7;2 9 2 9;5 5 3 3;8 1 8 1;6 2 6 2;7 3.6 7 3.6];
%cSH=[.1 .2 .2 .4 .4 .6 .3 .7 .5 .5];
%z=0;
%for i=1:10
    %z=z-((u-aSH(i,:))*(u-aSH(i,:))'+cSH(i))^(-1);
%end




%F����[-10��10]

%x=u(1,1);  
%y=u(1,2);  
%temp=x^2+y^2;  
%z=0.5+(sin(sqrt(temp))^2-0.5)/(1+0.001*temp)^2;
%for i=1:30
  % dim=size(u,2);
   %temp= i-1/dim-1;
   %  z= sum((10.^(6*temp)).* u.^2);      
    
%end   

notes:

%这个符号在matlab中表示注释的意思快捷键是ctrl+r，取消注释是ctrl+l，需要用测试函数就取消注释好了

原文地址：https://www.cnblogs.com/gaowenxingxing/p/11866636.html