matlab练习程序（粒子群优化PSO）

算法没有和图像处理直接相关，不过对于图像分类中的模式识别相关算法，也许会用到这个优化算法。

算法步骤：

1.首先确定粒子个数与迭代次数。

2.对每个粒子随机初始化位置与速度。

3.采用如下公式更新每个粒子的位置与速度。

Px=Px+Pv*t;　%位置更新公式　

Pv=Pv+(c1*rand*(Gx-Px))+(c2*rand*(PBx-Px)); %速度更新公式

这里c1和c2是加速因子，和梯度下降算法那里的加速因子我感觉很类似。

Gx是粒子群中最佳粒子的位置，PBx为当前粒子最佳位置。

4.每次迭代，首先检查新粒子适应度是否高于原最优适应度，如果高于则对自己的位置和适应度进行更新。然后再判断此粒子适应度是否高于全局最优粒子，如果高于则更新全局最优粒子适应度和位置。

因为自己不是主要研究这方面算法的，所以还有一些疑问（自问自答？）。

1.算法需要目标函数，如果没有目标函数怎么办。也许就不用这个算法了，或者其他什么算法先求出了目标函数了。

2.既然给了目标函数，那么直接遍历所有值再max()应该就能求得最佳位置。而PSO算法是不是只是为了减少运算量，比如我这里200*200的矩阵，本来需要计算40000次函数，而PSO只计算了100次函数就得到近似最优解了。

难怪叫优化算法，反正我暂时只能这样理解了，其他细节代码注释的很清楚了。

下图展示了一个PSO的运行结果，目标函数是高斯函数，绿点代表最佳粒子的位置：

matlab代码如下：

main.m

clear all;close all;clc;

[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;   %粒子群粒子个数

%初始化粒子群，定义结构体
%结构体中八个元素，分别是粒子坐标，粒子速度，粒子适应度，粒子最佳适应度，粒子最佳坐标
par=struct([]);　　　　　　　　　　
for i=1:n
    par(i).x=-100+200*rand();   %[-100 100]对x位置随机初始化
    par(i).y=-100+200*rand();   %[-100 100]对y位置随机初始化
    par(i).vx=-1+2*rand();      %[-1 1]对vx速度随机初始化
    par(i).vy=-1+2*rand();      %[-1 1]对vy速度随机初始化
    par(i).fit=0;               %粒子适应度为0初始化
    par(i).bestfit=0;           %粒子最佳适应度为0初始化
    par(i).bestx=par(i).x;      %粒子x最佳位置初始化
    par(i).besty=par(i).y;      %粒子y最佳位置初始化
end
par_best=par(1);    %初始化粒子群中最佳粒子

for k=1:10    
    plot3(par_best.x+100,par_best.y+100,par_best.fit,'g*'); %画出最佳粒子的位置，+100为相对偏移
    for p=1:n
        [par(p) par_best]=update_par(par(p),par_best);  %更新每个粒子信息         
    end  
end

update_par.m

function [par par_best]=update_par(par,par_best)
    
    %Px=Px+Pv*t,这里t=1,Px为当前粒子的位置，Pv为当前粒子的速度
    par.x=par.x+par.vx;   
    par.y=par.x+par.vy;   
    
    par.fit=compute_fit(par);    %计算当前粒子适应度
    
    %Pv=Pv+(c1*rand*(Gx-Px))+(c2*rand*(PBx-Px))
    %这里c1,c2为加速因子
    %Gx为具有最佳适应度粒子的位置
    %PBx为当前粒子的最佳位置
    c1=1;
    c2=1;
    par.vx=par.vx+c1*rand()*(par_best.x-par.x)+c2*rand()*(par.bestx-par.x);   
    par.vy=par.vy+c1*rand()*(par_best.y-par.y)+c2*rand()*(par.besty-par.y);
 
    if par.fit>par.bestfit      %如果当前粒子适应度要好于当前粒子最佳适应度
        par.bestfit=par.fit;    %则更新当前粒子最佳适应度
        par.bestx=par.x;        %更新当前粒子最佳位置
        par.besty=par.y;
        if par.bestfit>par_best.fit     %如果当前粒子最佳适应度好于最佳粒子适应度
            par_best.fit=par.bestfit;   %则更新最佳粒子适应度
            par_best.x=par.x;           %更新最佳粒子位置
            par_best.y=par.y;
        end
    end

end

compute_fit.m

function re=compute_fit(par)
    x=par.x;
    y=par.y;
    sigma=50;
    if x<-100 || x>100 || y<-100 || y>100
        re=0;        %超出范围适应度为0
    else            %否则适应度按目标函数求解
        re= (1/(2*pi*sigma^2))*exp(-(x.^2+y.^2)/(2*sigma^2)); 
    end
end