现代算法(二) PSO 算法

　　概念不复述了,网上都有,这里讨论标准PSO

　　w = 0.9; c1 = 2; c2 = 2; Vmax = 3;

　　种群大小 M = 20;迭代次数 N = 1000

function res=PSO(N,M,w,c1,c2,Vmax)
%w = 0.9; c1 = 2; c2 = 2; Vmax = 2;
%P(:,1:2) = 200*rand(M,2)-100;
P(:,1:2) = 3*rand(M,2);
P(:,3:4) = 2*Vmax*rand(M,2)-Vmax;
[g,pbest]= pb(P);
j = 0;
res = [];
tic 
while j < N     
    for i=1:M         % exchange message between each other         
        if(fitness(P(i,1),P(i,2)) < fitness(pbest(i,2),pbest(i,3)))             
            pbest(i,2:3) = P(i,1:2);             
            pbest(i,1) = fitness(P(i,1),P(i,2));         
        end
        g = pbest(i,:);           % update g by Piself         
        g = betterNei(pbest,i,g); % update g vs neighbour
 
        % update position and speed of Pi         % update Px         
        P(i,3:4) = w*P(i,3:4)+c1*rand*(pbest(i,1:2)-P(i,1:2)) ...
            + c2*rand*(g(2:3)-P(i,1:2));         
        P(i,3:4) = P2Vmax(P(i,3:4),Vmax);         % update Pv         
        P(i,1:2) = P(i,1:2) + P(i,3:4);     
    end
    j = j + 1; 
    res(j) = g(1);
end
time = toc;
result(g,time);
end

 　　Fitness

function f=fitness(x1,x2)
f = 2*x1.^2-4*x1.*x2+4*x2.^2-6*x1-3*x2....
    +1000*(max(x1+x2-3,0))+1000*(max(4*x1+x2-9,0));%1000 可以换成任意非常大数

 　　P2Max

function P = P2Vmax(P,Vmax)
% > Vmax 
P(P > Vmax) = Vmax; 
% < -Vmax 
P(P < -Vmax) = -Vmax;
end

　　pb

function [g,pbes]=pb(P)
pbes(:,1) = fitness(P(:,1),P(:,2));
pbes(:,2:3) = P(:,1:2);
[~,temp1] = min(pbes);%找出列向量最小的位置
temp1 = min(temp1);%最小中的最小
g(1) = pbes(temp1);
g(2:3) = P(temp1,1:2);
end

　　result

function result(g,time) 
fprintf('      x1: %10.12f
',g(2)); 
fprintf('      x2: %10.12f
',g(3)); 
fprintf('f(x1,x2): %10.12f
',g(1)); 
fprintf('    time: %10.12f s
',time); 
end

　　Better_Nei

% 在邻居中寻找更好的点代替
function g = betterNei(pbes,i,g)
M = size(pbes,1);
pbe(1:i-1,:) = pbes(1:i-1,2:3);
% i的附近
pbe(i:M-1,:) = pbes(i+1:M,2:3);
temp0 = fitness(pbe(:,1),pbe(:,2));
[~,temp] = min(temp0);
temp = min(temp);
g(1) = temp0(temp);
g(2:3) = pbe(temp,:);
end

　　可视化

res1 = PSO(1000,20,0.9,2,2,3);
res2 = PSO(1000,20,0.9,2,2,2);
res3 = PSO(1000,20,0.9,2,2,4);
plot(1:1000,res1,1:1000,res2,1:1000,res3)
axis([1,1000,-12,0])
legend('Vmax=3','Vmax=2','Vmax=4','Location','NorthEast')

结果:

完.hahaha