一、簡介

1 粒子群演算法的概念
粒子群優化演算法(PSO：Particle swarm optimization) 是一種進化計算技術（evolutionary computation）。源於對鳥群捕食的行為研究。粒子群優化演算法的基本思想：是通過群體中個體之間的協作和資訊共享來尋找最優解．
PSO的優勢：在於簡單容易實現並且沒有許多引數的調節。目前已被廣泛應用於函式優化、神經網路訓練、模糊系統控制以及其他遺傳演算法的應用領域。

2 粒子群演算法分析
2.1基本思想
粒子群演算法通過設計一種無質量的粒子來模擬鳥群中的鳥，粒子僅具有兩個屬性：速度和位置，速度代表移動的快慢，位置代表移動的方向。每個粒子在搜尋空間中單獨的搜尋最優解，並將其記為當前個體極值，並將個體極值與整個粒子群裡的其他粒子共享，找到最優的那個個體極值作為整個粒子群的當前全域性最優解，粒子群中的所有粒子根據自己找到的當前個體極值和整個粒子群共享的當前全域性最優解來調整自己的速度和位置。下面的動圖很形象地展示了PSO演算法的過程：

2 更新規則
PSO初始化為一群隨機粒子(隨機解)。然後通過迭代找到最優解。在每一次的迭代中，粒子通過跟蹤兩個“極值”(pbest，gbest)來更新自己。在找到這兩個最優值後，粒子通過下面的公式來更新自己的速度和位置。

公式(1)的第一部分稱為【記憶項】，表示上次速度大小和方向的影響；公式(1)的第二部分稱為【自身認知項】，是從當前點指向粒子自身最好點的一個向量，表示粒子的動作來源於自己經驗的部分；公式(1)的第三部分稱為【群體認知項】，是一個從當前點指向種群最好點的向量，反映了粒子間的協同合作和知識共享。粒子就是通過自己的經驗和同伴中最好的經驗來決定下一步的運動。以上面兩個公式為基礎，形成了PSO的標準形式。

公式(2)和 公式(3)被視為標準PSO演算法。
3 PSO演算法的流程和虛擬碼

二、原始碼

tic % 計時器
%% 清空環境，準備資料
close all
clear
clc
format compact
% 載入測試資料wine,其中包含的資料類別數為3；wine:178*13的矩陣,wine_labes:178*1的列向量
load wine
% 選定訓練集和測試集
% 將第一類的1-30,第二類的60-95,第三類的131-153做為訓練集
train_wine = [wine(1:30,:);wine(60:95,:);wine(131:153,:)];
% 相應的訓練集的標籤也要分離出來
train_wine_labels = [wine_labels(1:30);wine_labels(60:95);wine_labels(131:153)];
% 將第一類的31-59,第二類的96-130,第三類的154-178做為測試集
test_wine = [wine(31:59,:);wine(96:130,:);wine(154:178,:)];
% 相應的測試集的標籤也要分離出來
test_wine_labels = [wine_labels(31:59);wine_labels(96:130);wine_labels(154:178)];
% 資料預處理
% 資料預處理,將訓練集和測試集歸一化到[0,1]區間
[mtrain,ntrain] = size(train_wine);
[mtest,ntest] = size(test_wine);
dataset = [train_wine;test_wine];

[dataset_scale,ps] = mapminmax(dataset',0,1);
dataset_scale = dataset_scale';

train_wine = dataset_scale(1:mtrain,:);
test_wine = dataset_scale( (mtrain+1):(mtrain+mtest),: );
%% PSO引數尋優
% 粒子位置和速度更新引數
c1=1.5;
c2=1.5; 
maxgen=50; % 迭代次數
sizepop=20; % 種群規模（粒子數）
% 目標函式資訊
dim=2; % 待優化引數個數
xub=[100,100]; % 引數取值上界
xlb=[0.01,0.01]; % 引數取值下界
vub=[5,5]; % 速度上界
vlb=[-5,-5]; % 速度下界
%% 初始化
% 粒子位置初始化
xRange=repmat((xub-xlb),[sizepop,1]);
xLower=repmat(xlb,[sizepop,1]);
pop=rand(sizepop,dim).*xRange+xLower;
% 粒子速度初始化
vRange=repmat((vub-vlb),[sizepop,1]);
vLower=repmat(vlb,[sizepop,1]);
V=rand(sizepop,dim).*vRange+vLower;
% 計算初始化的目標函式值（適應度函式值）
fitness=ones(sizepop,1);
for k=1:sizepop
    fitness(k)=objfun(pop(k,:),train_wine_labels,train_wine,test_wine_labels,test_wine);
end
%% 尋找初始極值
[bestfitness,bestindex]=min(fitness);
zbest=pop(bestindex,:); % 全域性最優解的位置
gbest=pop; % 個體位置
fitnessgbest=fitness; % 個體適應度值
fitnesszbest=bestfitness; % 全域性最優適應度值
%% 迭代尋優
yy=ones(sizepop,1); % 用於儲存每次迭代的最優目標函式值
for k=1:maxgen
    % 粒子位置和速度更新
    for m=1:sizepop
        % 粒子速度更新
        sol=V(m,:)+c1*rand*(gbest(m,:)-pop(m,:))+c2*rand*(zbest-pop(m,:));
        % 確保粒子速度取值範圍不越界
        ind=find(sol<vlb);
        sol(ind)=vlb(ind);
        ind=find(sol>vub);
        sol(ind)=vub(ind);
        V(m,:)=sol;
        % 粒子位置更新
        sol=pop(m,:)+0.5*V(m,:);
        % 確保粒子位置取值範圍不越界
        ind=find(sol<xlb);
        sol(ind)=xlb(ind);
        ind=find(sol>xub);
        sol(ind)=xub(ind);
        pop(m,:)=sol;
        % 更新粒子適應度值
        fitness(m)=objfun(pop(m,:),train_wine_labels,train_wine,test_wine_labels,test_wine);
    end
    
    % 個體極值及位置和群體極值及位置更新
    for m=1:sizepop
        % 個體極值及其位置更新
        if fitness(m)<fitnessgbest(m)
            gbest(m,:)=pop(m,:);
            fitnessgbest(m)=fitness(m);
        end
        % 群體極值及其位置更新
        if fitness(m)<fitnesszbest
            zbest=pop(m,:);
            fitnesszbest=fitness(m);
        end
    end
    
    % 記錄每代最優值
    yy(k)=fitnesszbest;
end
%% 列印引數選擇結果
bestc=zbest(1);
bestg=zbest(2);
disp('列印選擇結果');
str=sprintf('Best c = %g，Best g = %g',bestc,bestg);
disp(str)
%% 利用最佳的引數進行SVM網路訓練
cmd_gwosvm = ['-c ',num2str(bestc),' -g ',num2str(bestg)];
model_gwosvm = svmtrain(train_wine_labels,train_wine,cmd_gwosvm);
%% SVM網路預測
[predict_label,accuracy] = svmpredict(test_wine_labels,test_wine,model_gwosvm);
% 列印測試集分類準確率
total = length(test_wine_labels);
right = sum(predict_label == test_wine_labels);
disp('列印測試集分類準確率');
 
%% 結果分析
% 測試集的實際分類和預測分類圖
figure;
hold on;
plot(test_wine_labels,'o');
plot(predict_label,'r*');
xlabel('測試集樣本','FontSize',12);
ylabel('類別標籤','FontSize',12);
legend('實際測試集分類','預測測試集分類');
title('測試集的實際分類和預測分類圖','FontSize',12);
grid on
snapnow
%% 顯示程式執行時間
toc
% 粒子速度初始化
vRange=repmat((vub-vlb),[sizepop,1]);
vLower=repmat(vlb,[sizepop,1]);
V=rand(sizepop,dim).*vRange+vLower;
% 計算初始化的目標函式值（適應度函式值）
fitness=ones(sizepop,1);
for k=1:sizepop
    fitness(k)=fun(pop(k,:));
end
%% 尋找初始極值
[bestfitness,bestindex]=min(fitness);
zbest=pop(bestindex,:); % 全域性最優解的位置
gbest=pop; % 個體位置
fitnessgbest=fitness; % 個體適應度值
fitnesszbest=bestfitness; % 全域性最優適應度值
%% 迭代尋優
yy=ones(sizepop,1); % 用於儲存每次迭代的最優目標函式值
for k=1:maxgen
    % 粒子位置和速度更新
    for m=1:sizepop
        % 粒子速度更新
        sol=V(m,:)+c1*rand*(gbest(m,:)-pop(m,:))+c2*rand*(zbest-pop(m,:));
        % 確保粒子速度取值範圍不越界
        ind=find(sol<vlb);
        sol(ind)=vlb(ind);
        ind=find(sol>vub);
        sol(ind)=vub(ind);
        V(m,:)=sol;
        % 粒子位置更新
        sol=pop(m,:)+0.5*V(m,:);
        % 確保粒子位置取值範圍不越界
        ind=find(sol<xlb);
        sol(ind)=xlb(ind);
        ind=find(sol>xub);
        sol(ind)=xub(ind);
        pop(m,:)=sol;
        % 更新粒子適應度值
        fitness(m)=fun(pop(m,:));
    end
    
    % 個體極值及位置和群體極值及位置更新
    for m=1:sizepop
        % 個體極值及其位置更新
        if fitness(m)<fitnessgbest(m)
            gbest(m,:)=pop(m,:);
            fitnessgbest(m)=fitness(m);
        end
        % 群體極值及其位置更新
        if fitness(m)<fitnesszbest
            zbest=pop(m,:);
            fitnesszbest=fitness(m);
        end
    end
    
    % 記錄每代最優值
    yy(k)=fitnesszbest;
end
 

三、執行結果

四、備註

版本：2014a