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