计算智能----粒子群优化算法

    大自然是我们的老师，生物的进化过程、群体智能活动为我们设计一个又一个优化算法提供了灵感的源泉。粒子群优化算法（PSO）就是仿生算法的一个著名代表。它是一种群体智能的随机搜索算法。

            粒子群算法的两个重要公式分别是速度更新公式和位置更新公式。每个粒子在进化的过程中需要维护两个向量，一个是速度向量V=[v1,v2,....vn]，一个是位置向量X=[x1,x2,.....xn]，其中n代表问题的维数。粒子的速度是矢量，决定了粒子运动的方向和速率，而粒子的位置体现了所代表的解在空间中的位置，是评估该解质量的基础。算法还要求每个粒子维护自身的一个历史最优位置向量，用pbest表示，群体还要维护一个全局的最优位置，记为gbest，gbest从pbest中选出，是pbest中最优的那个。这个全局最优向量和pbest引导粒子向最优或者局部最优收敛。粒子群优化算法和遗传算法相比，没有了选择算子、交配算子和变异算子，仅仅通过速度更新公式和位置更新公式不断的进化得到最优解，算法相对简单，运行效率高。

             简要步骤如下：

             （1）初始化所有的个体，初始化速度和位置，并将个体的pbest设置为当前位置，而群体中最优的个体当作当前的gbest。

             （2）在每一代进化中，计算各个粒子的适应度函数，适应度一般与目标函数相关。

             （3）比较粒子的当前与历史适应度，如果当前更好则更新pbest，如果历史更好则不更新pbest。

             （4）同样的，在pbest中选出一个待比较的gbest与历史gbest相比较，若当前gbest更好则更新，反之则保持。

             （5）对每个粒子按照如下公式更新其速度和位置向量

               V(t+1) = W * V(t) + c1*rand*(pbest - Xi) + c2*rand*(gbest - Xi);

               X(t+1) = X(t) + V(t);

           W为惯性系数，一般初始化为0.9，随着迭代次数增加逐渐线性降低至0.4；rand为0~1之间的随机数；c1，c2为加速系数也成为学习因子，一般取c1=c2=2；同时在实际问题中还要对速度向量V和位置向量X进行限制，限制在可行域内。

              伪代码如下：

procedure PSO

 for each particle i

      initialize V and X for each particle i

      evaluate practice i and set pbest = X

end

gbest = pbest(min(f(pbest)))

while not stop

    for i = 1 to N

        update the V and X

        evaluate

        if fit(Xi) < fit(pbest)

            pbest = Xi

        if fit(pbest) < fit(gbest)

            gbest = pbest

end

end

 下面给出一个例子，

min (x1-1)^2 + (x2+4)^2

S.T： x1^2+x2^2 < 4

% draw picture
close all;
clear
clc
tic;
theta = 0 : pi/200 : 2*pi;
a = 2*cos(theta);
b = 2*sin(theta);
plot(a,b,'k');
grid on;
hold on;
plot(1,-4,'r*')
% initial index
w = 0.9;
c1 = 2;
c2 = 2;
popsize = 5;
% max step
max_step = 20;
%initial population
x1=rand(popsize,1)*4-2;
temp=sqrt(4-x1.^2);
x2 = 2.*temp.*rand(popsize,1) - temp;
plot(x1,x2,'o');
%initial pbest & gbest & velocity
pbest1 = rand(popsize,1)*4-2;
pbest2 = 2 .* sqrt(4-x1.^2) .* rand(popsize,1) - sqrt(4-x1.^2);
v1 = rand(popsize,1);
v2 = rand(popsize,1);
gbest1 = pbest1;
gbest2 = pbest2;

% count variable
count = 0;
y = (x1-1).*(x1-1)+(x2+4).*(x2+4);
f = (pbest1 - 1).^2 + (pbest2 + 4).^2;
pro_f = y;
while count < max_step
    count=count+1;
    w = max(w - (count/max_step)*w,0.4);
    % update x1 & x2 & velocity 
    v1 = w .* v1 + c1.*rand(popsize,1).*(pbest1 - x1) + c2 .* rand(popsize,1).*(gbest1 - x1);
    v2 = w .* v2 + c1.*rand(popsize,1).*(pbest2 - x2) + c2 .* rand(popsize,1).*(gbest2 - x2);
    x1 = max(-2,min(x1 + v1 , 2));
    bound = sqrt(4-x1.^2);
    x2 = max(-bound , min(x2 + v2 , bound));
    % computing cost function
    f = (x1 - 1).^2 + (x2 + 4).^2;
    % find pbest1 & pbest2
    for i = 1:1:popsize
        if(f(i) < pro_f)
            pbest1(i) = x1(i);
            pbest2(i) = x2(i); 
        end
    end
    pro_f = f;
    % update gbest 
   num = find(f==min(f));
   gbest1 = pbest1(num(1));
   gbest2 = pbest2(num(1));
   repmat(gbest1,popsize,1);
   repmat(gbest2,popsize,1);
   plot(pbest1,pbest2,'+','markeredgecolor','m');
   plot(gbest1,gbest2,'s','markeredgecolor','b');
   hold on;
   pause(0.1);
   
   % judge stop condition
%    if abs(pro_f - f) < 0.001
%        break;
%    else
%        pro_f = f;
%    end
end
X = gbest1;
Y = gbest2;
plot(X,Y,'*','markeredgecolor','r');
[X,Y]
min(f)
toc;
disp(['运行时间: ',num2str(toc)]);