一、简介

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