注：本人参考http://www.cnblogs.com/tiandsp/category/348031.html来实现的
算法步骤：
1.首先确定粒子个数与迭代次数。
2.对每个粒子随机初始化位置与速度。
3.采用如下公式更新每个粒子的位置与速度。
Px=Px+Pv*t; %位置更新公式
Pv=Pv+(c1*rand*(Gx-Px))+(c2*rand*(PBx-Px)); %速度更新公式
这里c1和c2是加速因子，和梯度下降算法那里的加速因子我感觉很类似。
Gx是粒子群中最佳粒子的位置，PBx为当前粒子最佳位置。
4.每次迭代，首先检查新粒子适应度是否高于原最优适应度，如果高于则对自己的位置和适应度进行更新。然后再判断此粒子适应度是否高于全局最优粒子，如果高于则更新全局最优粒子适应度和位置。
clc;
clear all;
close all;
%%
[x,y] = meshgrid(-100:100,-100:100);
sigma = 50;
img = (1/(2*pi*sigma^2))*exp(-(x.^2+y.^2)/(2*sigma^2));
mesh(img)
hold on
%% 初始化粒子群，定义结构体
%结构体中八个元素，分别是粒子坐标，粒子速度，粒子适应度，
% 粒子最佳适应度，粒子最佳坐标
n = 10; %初始例子的个数
pos = struct([]);
for i = 1:n
pos(i).x = -100 + 200*rand;
pos(i).y = -100 + 200*rand;
pos(i).vx = -1 + 2*rand;
pos(i).vy = -1 + 2*rand;
pos(i).fit = 0;
pos(i).bestfit = 0;
pos(i).bestx = pos(i).x;
pos(i).besty = pos(i).y;
end
pos_best = pos(1);
k = 100;
for i = 1:k
plot3(pos_best.x + 100,pos_best.y + 100,pos_best.fit,'k.','markersize',40);
for j = 1:n
[pos(j),pos_best] = update_pos(pos(j),pos_best);
end
drawnow
end
function [pos,pos_best] = update_pos(pos,pos_best)
pos.x = pos.x + pos.vx;
pos.y = pos.y + pos.vy;
pos.fit = compute_fit(pos);
%Pv=Pv+(c1*rand*(pos_best(群体).x-pos.x))+(c2*rand*(pos.bestx-pos.x))
%这里c1,c2为加速因子
%Gx为具有最佳适应度粒子的位置
%PBx为当前粒子的最佳位置
c1 = 1;
c2 = 1;
pos.vx = pos.vx + (c1*rand*(pos_best.x -pos.x)) + (c2*rand*(pos.bestx-pos.x));
pos.vy = pos.vy + (c1*rand*(pos_best.y -pos.y)) + (c2*rand*(pos.besty-pos.y));
if pos.fit > pos.bestfit
pos.bestfit = pos.fit;
pos.bestx = pos.x;
pos.besty = pos.y;
if pos_best.fit < pos.bestfit
pos_best.fit = pos.bestfit;
pos_best.x = pos.bestx ;
pos_best.y = pos.besty ;
end
end
end
function pos_fit = compute_fit(pos)
x = pos.x;
y = pos.y;
sigma = 50;
if x <-100 || x > 100 || y <-100 || y >100
pos_fit = 0;
else
pos_fit = (1/(2*pi*sigma^2))*exp(-(x.^2+y.^2)/(2*sigma^2));
end