一、简介
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算法的流程和伪代码
二、源代码
%% 清空环境变量
function chapter_PSO
close all;
clear;
clc;
format compact;
%% 数据提取
% 画出测试数据的分维可视化图
figure
subplot(3,5,1);
hold on
for run = 1:178
plot(run,wine_labels(run),'*');
end
xlabel('样本','FontSize',10);
ylabel('类别标签','FontSize',10);
title('class','FontSize',10);
for run = 2:14
subplot(3,5,run);
hold on;
str = ['attrib ',num2str(run-1)];
for i = 1:178
plot(i,wine(i,run-1),'*');
end
xlabel('样本','FontSize',10);
ylabel('属性值','FontSize',10);
title(str,'FontSize',10);
end
% 选定训练集和测试集
% 将第一类的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];
% mapminmax为MATLAB自带的归一化函数
[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),: );
%% 选择最佳的SVM参数c&g
[bestacc,bestc,bestg] = psoSVMcgForClass(train_wine_labels,train_wine);
% 打印选择结果
disp('打印选择结果');
str = sprintf( 'Best Cross Validation Accuracy = %g%% Best c = %g Best g = %g',bestacc,bestc,bestg);
disp(str);
%% 利用最佳的参数进行SVM网络训练
cmd = ['-c ',num2str(bestc),' -g ',num2str(bestg)];
model = svmtrain(train_wine_labels,train_wine,cmd);
%% SVM网络预测
[predict_label,accuracy] = svmpredict(test_wine_labels,test_wine,model);
% 打印测试集分类准确率
total = length(test_wine_labels);
right = sum(predict_label == test_wine_labels);
disp('打印测试集分类准确率');
str = sprintf( 'Accuracy = %g%% (%d/%d)',accuracy(1),right,total);
disp(str);
%% 结果分析
% 测试集的实际分类和预测分类图
% 通过图可以看出只有三个测试样本是被错分的
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;
%% 子函数 psoSVMcgForClass.m
function [bestCVaccuarcy,bestc,bestg,pso_option] = psoSVMcgForClass(train_label,train,pso_option)
% psoSVMcgForClass
% 参数初始化
if nargin == 2
pso_option = struct('c1',1.5,'c2',1.7,'maxgen',200,'sizepop',20, ...
'k',0.6,'wV',1,'wP',1,'v',5, ...
'popcmax',10^2,'popcmin',10^(-1),'popgmax',10^3,'popgmin',10^(-2));
end
% c1:初始为1.5,pso参数局部搜索能力
% c2:初始为1.7,pso参数全局搜索能力
% maxgen:初始为200,最大进化数量
% sizepop:初始为20,种群最大数量
% k:初始为0.6(k belongs to [0.1,1.0]),速率和x的关系(V = kX)
% wV:初始为1(wV best belongs to [0.8,1.2]),速率更新公式中速度前面的弹性系数
% wP:初始为1,种群更新公式中速度前面的弹性系数
% v:初始为3,SVM Cross Validation参数
% popcmax:初始为100,SVM 参数c的变化的最大值.
% popcmin:初始为0.1,SVM 参数c的变化的最小值.
% popgmax:初始为1000,SVM 参数g的变化的最大值.
% popgmin:初始为0.01,SVM 参数c的变化的最小值.
Vcmax = pso_option.k*pso_option.popcmax;
Vcmin = -Vcmax ;
Vgmax = pso_option.k*pso_option.popgmax;
Vgmin = -Vgmax ;
eps = 10^(-3);
% 产生初始粒子和速度
for i=1:pso_option.sizepop
% 随机产生种群和速度
pop(i,1) = (pso_option.popcmax-pso_option.popcmin)*rand+pso_option.popcmin;
pop(i,2) = (pso_option.popgmax-pso_option.popgmin)*rand+pso_option.popgmin;
V(i,1)=Vcmax*rands(1,1);
V(i,2)=Vgmax*rands(1,1);
% 计算初始适应度
cmd = ['-v ',num2str(pso_option.v),' -c ',num2str( pop(i,1) ),' -g ',num2str( pop(i,2) )];
fitness(i) = svmtrain(train_label, train, cmd);
fitness(i) = -fitness(i);
end
% 找极值和极值点
[global_fitness bestindex]=min(fitness); % 全局极值
local_fitness=fitness; % 个体极值初始化
global_x=pop(bestindex,:); % 全局极值点
local_x=pop; % 个体极值点初始化
% 每一代种群的平均适应度
avgfitness_gen = zeros(1,pso_option.maxgen);
% 迭代寻优
for i=1:pso_option.maxgen
for j=1:pso_option.sizepop
%速度更新
V(j,:) = pso_option.wV*V(j,:) + pso_option.c1*rand*(local_x(j,:) - pop(j,:)) + pso_option.c2*rand*(global_x - pop(j,:));
if V(j,1) > Vcmax
V(j,1) = Vcmax;
end
if V(j,1) < Vcmin
V(j,1) = Vcmin;
end
if V(j,2) > Vgmax
V(j,2) = Vgmax;
end
if V(j,2) < Vgmin
V(j,2) = Vgmin;
end
%种群更新
pop(j,:)=pop(j,:) + pso_option.wP*V(j,:);
if pop(j,1) > pso_option.popcmax
pop(j,1) = pso_option.popcmax;
end
if pop(j,1) < pso_option.popcmin
pop(j,1) = pso_option.popcmin;
end
if pop(j,2) > pso_option.popgmax
pop(j,2) = pso_option.popgmax;
end
if pop(j,2) < pso_option.popgmin
pop(j,2) = pso_option.popgmin;
end
% 自适应粒子变异
if rand>0.5
k=ceil(2*rand);
if k == 1
pop(j,k) = (20-1)*rand+1;
end
if k == 2
pop(j,k) = (pso_option.popgmax-pso_option.popgmin)*rand + pso_option.popgmin;
end
end
%适应度值
cmd = ['-v ',num2str(pso_option.v),' -c ',num2str( pop(j,1) ),' -g ',num2str( pop(j,2) )];
fitness(j) = svmtrain(train_label, train, cmd);
fitness(j) = -fitness(j);
cmd_temp = ['-c ',num2str( pop(j,1) ),' -g ',num2str( pop(j,2) )];
model = svmtrain(train_label, train, cmd_temp);
if fitness(j) >= -65
continue;
end
三、运行结果
四、备注
版本：2014a