现代算法(二) PSO 算法

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

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

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

 1 function res=PSO(N,M,w,c1,c2,Vmax)
 2 %w = 0.9; c1 = 2; c2 = 2; Vmax = 2;
 3 %P(:,1:2) = 200*rand(M,2)-100;
 4 P(:,1:2) = 3*rand(M,2);
 5 P(:,3:4) = 2*Vmax*rand(M,2)-Vmax;
 6 [g,pbest]= pb(P);
 7 j = 0;
 8 res = [];
 9 tic 
10 while j < N     
11     for i=1:M         % exchange message between each other         
12         if(fitness(P(i,1),P(i,2)) < fitness(pbest(i,2),pbest(i,3)))             
13             pbest(i,2:3) = P(i,1:2);             
14             pbest(i,1) = fitness(P(i,1),P(i,2));         
15         end
16         g = pbest(i,:);           % update g by Piself         
17         g = betterNei(pbest,i,g); % update g vs neighbour
18  
19         % update position and speed of Pi         % update Px         
20         P(i,3:4) = w*P(i,3:4)+c1*rand*(pbest(i,1:2)-P(i,1:2)) ...
21             + c2*rand*(g(2:3)-P(i,1:2));         
22         P(i,3:4) = P2Vmax(P(i,3:4),Vmax);         % update Pv         
23         P(i,1:2) = P(i,1:2) + P(i,3:4);     
24     end
25     j = j + 1; 
26     res(j) = g(1);
27 end
28 time = toc;
29 result(g,time);
30 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

1 function P = P2Vmax(P,Vmax)
2 % > Vmax 
3 P(P > Vmax) = Vmax; 
4 % < -Vmax 
5 P(P < -Vmax) = -Vmax;
6 end

1 function [g,pbes]=pb(P)
2 pbes(:,1) = fitness(P(:,1),P(:,2));
3 pbes(:,2:3) = P(:,1:2);
4 [~,temp1] = min(pbes);%找出列向量最小的位置
5 temp1 = min(temp1);%最小中的最小
6 g(1) = pbes(temp1);
7 g(2:3) = P(temp1,1:2);
8 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