PSO学习笔记
1.仿生思想
粒子群算法模拟鸟集群飞行觅食的行为,通过鸟之间的集体协作使群体达到目的。相对于其他演化方法,它具有原理简单,易于编程实现等特点。PSO数学描述如下:
初始化一群 随机粒子(随机解)。然后通过迭代找到最优解。在每一次的迭代中,粒子通过跟踪两个“极值”(Pbest,gbest)来更新自己。
在找到这两个最优值后,粒子通过下面的公式来更新自己的速度和位置。
[egin{array}{rl}{V_{t, j}(t+1)=w} & {V_{t, j}(t)+c_{1} imes r ext { and }() imesleft(p b e s t_{i, j}-x_{t, j}(t)
ight)} \ {} & {+c_{2} imes r a n d() imesleft(g b e s t_{i . j}-x_{i . j}(t)
ight)}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为惯性权因子。
#基本粒子群算法代码
```{r}
%% 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的收敛性能的
```{r}
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，需要用测试函数就取消注释好了]