粒子群算法（PSO）介绍及matlab实现

        PSO模拟鸟群的捕食行为。设想这样一个场景：一群鸟在随机搜索食物。在这个区域里只有一块食物。所有的鸟都不知道食物在那里。但是他们知道当前的位置离食物还有多远。那么找到食物的最优策略是什么呢。最简单有效的就是搜寻目前离食物最近的鸟的周围区域。

       由于每个鸟都不知道食物在哪里，但是却知道距离多远，所以每次可以知道整个鸟群哪只鸟距离目标最近，当然最近的距离也知道，这个时候这只鸟可以对整个群体发“信号”，然后整个鸟群就向着这只鸟方向飞行，多次迭代之后，整个鸟群就已经接近目标了

       这里的目标比较泛华，可以理解为求一个方程组在一定范围内的极值等等

       基本思想：

2.根据适应度算法计算粒子适应度值

3.根据粒子计算的适应度值寻找粒子当前搜索到的最优位置与整个粒子群当前的最优位置

       下面就利用粒子群来对一个方程组在给定范围内求最小值，以加深对粒子群算法的理解，当然该算法用处很多，远远不止这一个。

        方程组：0.4 * (x - 5)^2 + 0.3 * (y - 6)^2 - 0.2 ;    x,y的范围：[-100,100]， 这个时候目标就是求最小值了

D=2;%维数 一般是自变量个数 此时速度与位置都用两个变量表示
w = 1;
NP=12;%种群规模 D的五倍到十倍
%学习因子c1 c2 范围在0-4
c1 = 1.5; 
c2 = 1.5;  
gen_max=1000;%最大进化代数
bounds_p=ones(D,2);%定义位置的取值边界
bounds_p(:,1)=-100;%取值下界
bounds_p(:,2)=100;% 取值上界
bounds_v=ones(D,2);%定义速度的取值边界
bounds_v(:,1)=-100;%取值下界
bounds_v(:,2)=100;% 取值上界
x=ones(NP,1)*bounds_p(:,1)'+(ones(NP,1)*(bounds_p(:,2)'-bounds_p(:,1)')) .*(rand(NP,D));%初始种群位置信息，使位置都位于最大值与最小值间，最后形成NP行2列的矩阵
v=ones(NP,1)*bounds_v(:,1)'+(ones(NP,1)*(bounds_v(:,2)'-bounds_v(:,1)')) .*(rand(NP,D));%初始种群速度信息，使位置都位于最大值与最小值间，最后形成NP行2列的矩阵
count=1;%进化代数
cost=zeros(NP,1);%存放个体的适应值，初始化为NP行一列的零矩阵
cost(1)=fitnessrevisemodel(x(1,:));%计算第一个个体适应值
fitness_zbest=cost(1);%存放全局最优值，假设第一个个体便是最优值
zbest=x(1,:);%存放全局对应最优值的解

for i2=2:NP %计算初始种群中的最优值和最优值对应的解
    cost(i2)=fitnessrevisemodel(x(i2,:));
    if(cost(i2)<= fitness_zbest)
        fitness_zbest=cost(i2);
        zbest=x(i2,:);
    end
end

gbest = x; %存放个体最优值，因为此时个体仅仅产生了一次，所以就取当前这一次为最优值
fitness_gbest = cost; %存放个体最优适应值，因为此时个体仅仅产生了一次，所以就取当前这一次为最优适应值

%粒子群算法  
while(count<gen_max)%循环终止条件是达到最大的进化代数
    for i2=1:NP
        %速度更新
        v(i2,:) = w * v(i2, :) + c1 * rand * (gbest(i2, :) - x(i2, :)) + c2 * rand * (zbest - x(i2, :)); %速度更新公式，考虑个体最优值与总体最优值
        v(i2, find(v(i2, :) > 100)) = 100; %如果第i2个个体的速度（用一行两列矩阵表示，这一行任一列大于100均要改变）大于最大速度100，即改为100
        v(i2, find(v(i2, :) < -100)) = -100; %如果第i2个个体的速度（用一行两列矩阵表示，这一行任一列小于-100均要改变）小于最小速度100，即改为-100
        
        %位置更新
        x(i2, :) = x(i2, :) + v(i2, :); %速度更新公式
        x(i2, find(x(i2, :) > 100)) = 100;
        x(i2, find(x(i2, :) < -100)) = -100;
        
        %自适应变异（避免算法陷入局部最优）
        if rand > 0.8
            k = ceil(2 * rand);
            x(i2, k) = rand;
        end
        
        %因为上述的x的第i2个个体更新，所以第i2个个体适应值也要更新
        cost(i2) = fitnessrevisemodel(x(i2, :));
        
        %如果个体适应值小于即优于原来该最优个体适应值
        if cost(i2) < fitness_gbest(i2)
            gbest(i2, :) = x(i2, :); %更新该个体最优值
            fitness_gbest(i2) = cost(i2); %更新个体最优适应值
        end
        
        %如果个体适应值比总体最优适应值小，更新
        if cost(i2) < fitness_zbest
            zbest = cost(i2, :); %存放新的全局最优值
            fitness_zbest = cost(i2); %存放新的全局最优适应值
        end
    end
    yy(count) = fitness_zbest; %记录每一次迭代后的新的全局最优适应值，理论上迭代结束后，最后一次的yy(count)即为所求值
    count = count + 1;
end
yy