粒子群优化（PSO）算法
一.算法思想
粒子群算法（ Particle Swarm Optimization, PSO）最早是由Eberhart和Kennedy于1995年提出，它的基本概念源于对鸟群觅食行为的研究。设想这样一个场景：一群鸟在随机搜寻食物，在这个区域里只有一块食物，所有的鸟都不知道食物在哪里，但是它们知道当前的位置离食物还有多远。最简单有效的策略？寻找鸟群中离食物最近的个体来进行搜素。PSO算法就从这种生物种群行为特性中得到启发并用于求解优化问题。用一种粒子来模拟上述的鸟类个体，每个粒子可视为N维搜索空间中的一个搜索个体，粒子的当前位置即为对应优化问题的一个候选解，粒子的飞行过程即为该个体的搜索过程．粒子的飞行速度可根据粒子历史最优位置和种群历史最优位置进行动态调整．粒子仅具有两个属性：速度和位置，速度代表移动的快慢，位置代表移动的方向。每个粒子单独搜寻的最优解叫做个体极值，粒子群中最优的个体极值作为当前全局最优解。不断迭代，更新速度和位置。最终得到满足终止条件的最优解。
二.算法介绍
①问题抽象
鸟被抽象为没有质量和体积的微粒(点)，并延伸到N维空间，粒子i在N维空间的位置表示为矢量Xi＝(x1，x2，…，xN)，飞行速度表示为矢量Vi＝(v1，v2，…，vN)。每个粒子都有一个由目标函数决定的适应值(fitness value)，并且知道自己到目前为止发现的最好位置(pbest)和现在的位置Xi。这个可以看作是粒子自己的飞行经验。除此之外，每个粒子还知道到目前为止整个群体中所有粒子发现的最好位置(gbest)(gbest是pbest中的最好值)，这个可以看作是粒子同伴的经验。粒子就是通过自己的经验和同伴中最好的经验来决定下一步的运动。
②更新规则
③算法流程
1、初始化
首先，我们设置最大迭代次数，目标函数的自变量个数，粒子的最大速度，位置信息为整个搜索空间，我们在速度区间和搜索空间上随机初始化速度和位置，设置粒子群规模为M，每个粒子随机初始化一个飞翔速度。
2、 个体极值与全局最优解
定义适应度函数，个体极值为每个粒子找到的最优解，从这些最优解找到一个全局值，叫做本次全局最优解。与历史全局最优比较，进行更新。
3、 更新速度和位置的公式
4、 终止条件
（1）达到设定迭代次数；（2）代数之间的差值满足最小界限
以上就是最基本的一个标准PSO算法流程。和其它群智能算法一样，PSO算法在优化过程中，种群的多样性和算法的收敛速度之间始终存在着矛盾．对标准PSO算法的改进，无论是参数的选取、小生境技术的采用或是其他技术与PSO的融合，其目的都是希望在加强算法局部搜索能力的同时，保持种群的多样性，防止算法在快速收敛的同时出现早熟收敛。
④参数分析
参数：群体规模N，惯性因子 ，学习因子c1和c2，最大速度Vmax，最大迭代次数Gk。
群体规模N：一般取20～40，对较难或特定类别的问题可以取到100～200。
最大速度Vmax：决定当前位置与最好位置之间的区域的分辨率(或精度)。如果太快，则粒子有可能越过极小点；如果太慢，则粒子不能在局部极小点之外进行足够的探索，会陷入到局部极值区域内。这种限制可以达到防止计算溢出、决定问题空间搜索的粒度的目的。
权重因子：包括惯性因子和学习因子c1和c2。使粒子保持着运动惯性，使其具有扩展搜索空间的趋势，有能力探索新的区域。c1和c2代表将每个粒子推向pbest和gbest位置的统计加速项的权值。较低的值允许粒子在被拉回之前可以在目标区域外徘徊，较高的值导致粒子突然地冲向或越过目标区域。
⑤参数设置：
参数w,c1,c2的选择分别关系粒子速度的3个部分：惯性部分、社会部分和自身部分在搜索中的作用。如何选择、优化和调整参数，使得算法既能避免早熟又能比较快的收敛，对工程实践有着重要意义。
1. 惯性权重w描述粒子上一代速度对当前代速度的影响。w值较大，全局寻优能力强，局部寻优能力弱；反之，则局部寻优能力强。当问题空间较大时，为了在搜索速度和搜索精度之间达到平衡，通常做法是使算法在前期有较高的全局搜索能力以得到合适的种子，而在后期有较高的局部搜索能力以提高收敛精度。所以w不宜为一个固定的常数。
wmax最大惯性权重，wmin最小惯性权重，run当前迭代次数，runmax为算法迭代总次数。较大的w有较好的全局收敛能力，较小的w则有较强的局部收敛能力。因此，随着迭代次数的增加，惯性权重w应不断减少，从而使得粒子群算法在初期具有较强的全局收敛能力，而晚期具有较强的局部收敛能力。
2. 学习因子c2=0称为自我认识型粒子群算法，即“只有自我，没有社会”，完全没有信息的社会共享，导致算法收敛速度缓慢；学习因子c1=0称为无私型粒子群算法，即“只有社会，没有自我”，会迅速丧失群体多样性，容易陷入局部最优解而无法跳出；c1,c2都不为0，称为完全型粒子群算法，完全型粒子群算法更容易保持收敛速度和搜索效果的均衡，是较好的选择。
3. 群体大小m是一个整数，m很小时陷入局部最优解的可能性很大；m很大时PSO的优化能力很好，但是当群体数目增长至一定水平时，再增长将不再有显著作用，而且数目越大计算量也越大。群体规模m 一般取20～40，对较难或特定类别的问题 可以取到100～200。
4. 粒子群的最大速度Vmax对维护算法的探索能力与开发能力的平衡很重要，Vmax较大时，探索能力强，但粒子容易飞过最优解；Vmax较小时，开发能力强，但是容易陷入局部最优解。Vmax一般设为每维变量变化范围的10%-20%
恰当的选取算法的参数值可以改善算法的性能。
三.参数组合设置对算法影响的研究
、单独调整粒子群算法的权重
1.线性递减权重
Shi.Y认为较大的权重惯性有利于全局搜索，较小的权重有利于局部搜索，他提出的线性递减权重刚好满足这样的需求。一般而言，，
k为当前迭代次数，T为最大迭代次数，此时c1=c2=1.494475
主函数代码：
%% 清空环境
clc
clear
%% 参数初始化
%粒子群算法中的三个参数
c1 = 1.49445;%加速因子
c2 = 1.49445;
ws=0.9;
we=0.2;
%w=0.8 %惯性权重
maxgen=1000; % 进化次s数
sizepop=200; %种群规模
Vmax=1; %限制速度围
Vmin=-1;
popmax=5; %变量取值范围
popmin=-5;
dim=10; %适应度函数维数
func=1; %选择待优化的函数，1为Rastrigin,2为Schaffer,3为Griewank
Drawfunc(func);%画出待优化的函数，只画出二维情况作为可视化输出
%% 产生初始粒子和速度
for i=1:sizepop
%随机产生一个种群
pop(i,:)=popmax*rands(1,dim); %初始种群
V(i,:)=Vmax*rands(1,dim); %初始化速度
%计算适应度
fitness(i)=fun(pop(i,:),func); %粒子的适应度
end
%% 个体极值和群体极值
[bestfitness bestindex]=min(fitness);
gbest=pop(bestindex,:); %全局最佳
pbest=pop; %个体最佳
fitnesspbest=fitness; %个体最佳适应度值
fitnessgbest=bestfitness; %全局最佳适应度值
%% 迭代寻优
for i=1:maxgen
fprintf('第%d代，',i);
fprintf('最优适应度%f\n',fitnessgbest);
w=ws-(ws-we)/(i*maxgen);
for j=1:sizepop
%速度更新
V(j,:) = w*V(j,:) + c1*rand*(pbest(j,:) - pop(j,:)) + c2*rand*(gbest - pop(j,:)); %根据个体最优pbest和群体最优gbest计算下一时刻速度
V(j,find(V(j,:)>Vmax))=Vmax; %限制速度不能太大
V(j,find(V(j,:)<Vmin))=Vmin;
%种群更新
pop(j,:)=pop(j,:)+0.5*V(j,:); %位置更新
pop(j,find(pop(j,:)>popmax))=popmax;%坐标不能超出范围
pop(j,find(pop(j,:)<popmin))=popmin;
if rand>0.98 %加入变异种子，用于跳出局部最优值
pop(j,:)=rands(1,dim);
end
%更新第j个粒子的适应度值
fitness(j)=fun(pop(j,:),func);
end
for j=1:sizepop
%个体最优更新
if fitness(j) < fitnesspbest(j)
pbest(j,:) = pop(j,:);
fitnesspbest(j) = fitness(j);
end
%群体最优更新
if fitness(j) < fitnessgbest
gbest = pop(j,:);
fitnessgbest = fitness(j);
end
end
yy(i)=fitnessgbest;
end
%% 结果分析
figure;
plot(yy)
title('最优个体适应度','fontsize',12);
xlabel('进化代数','fontsize',12);ylabel('适应度','fontsize',12);
结果：由于PSO算法运行速度较快，结果相对比较随机，所以运行10次求平均值来作为最终结果
最优适应度=（4.974796+3.979836+2.984878+2.984883+3.979837+1.989918+3.979838+1.989924+3.979940+6.964174）/10=3.7808024
结论分析： 惯性权重w描述粒子上一代速度对当前代速度的影响。w值较大，全局寻优能力强，局部寻优能力弱；反之，则局部寻优能力强。当问题空间较大时，为了在搜索速度和搜索精度之间达到平衡，通常做法是使算法在前期有较高的全局搜索能力以得到合适的种子，而在后期有较高的局部搜索能力以提高收敛精度。所以w不宜为一个固定的常数。wmax最大惯性权重，wmin最小惯性权重，run当前迭代次数，runmax为算法迭代总次数。较大的w有较好的全局收敛能力，较小的w则有较强的局部收敛能力。因此，随着迭代次数的增加，惯性权重w应不断减少，从而使得粒子群算法在初期具有较强的全局收敛能力，而晚期具有较强的局部收敛能力。
②单独调整异步学习因子
此时w=0.8
代码：
%% 清空环境
clc
clear
%% 参数初始化
%粒子群算法中的三个参数
c1 = 1.49445;%加速因子
c2 = 1.49445;
%ws=0.9;
%we=0.2;
w=0.8 %惯性权重
cstart1=2.5;
cend1=0.5;
cstart2=0.5;
cend2=2.5;
maxgen=1000; % 进化次s数
sizepop=200; %种群规模
Vmax=1; %限制速度围
Vmin=-1;
popmax=5; %变量取值范围
popmin=-5;
dim=10; %适应度函数维数
func=1; %选择待优化的函数，1为Rastrigin,2为Schaffer,3为Griewank
Drawfunc(func);%画出待优化的函数，只画出二维情况作为可视化输出
%% 产生初始粒子和速度
for i=1:sizepop
%随机产生一个种群
pop(i,:)=popmax*rands(1,dim); %初始种群
V(i,:)=Vmax*rands(1,dim); %初始化速度
%计算适应度
fitness(i)=fun(pop(i,:),func); %粒子的适应度
end
%% 个体极值和群体极值
[bestfitness bestindex]=min(fitness);
gbest=pop(bestindex,:); %全局最佳
pbest=pop; %个体最佳
fitnesspbest=fitness; %个体最佳适应度值
fitnessgbest=bestfitness; %全局最佳适应度值
%% 迭代寻优
for i=1:maxgen
fprintf('第%d代，',i);
fprintf('最优适应度%f\n',fitnessgbest);
c1=cstart1+(cend1-cstart1)*(i/maxgen);
c2=cstart2+(cend2-cstart2)*(i/maxgen);
%w=ws-(ws-we)/(i*maxgen);
for j=1:sizepop
%速度更新
V(j,:) = w*V(j,:) + c1*rand*(pbest(j,:) - pop(j,:)) + c2*rand*(gbest - pop(j,:)); %根据个体最优pbest和群体最优gbest计算下一时刻速度
V(j,find(V(j,:)>Vmax))=Vmax; %限制速度不能太大
V(j,find(V(j,:)<Vmin))=Vmin;
%种群更新
pop(j,:)=pop(j,:)+0.5*V(j,:); %位置更新
pop(j,find(pop(j,:)>popmax))=popmax;%坐标不能超出范围
pop(j,find(pop(j,:)<popmin))=popmin;
if rand>0.98 %加入变异种子，用于跳出局部最优值
pop(j,:)=rands(1,dim);
end
%更新第j个粒子的适应度值
fitness(j)=fun(pop(j,:),func);
end
for j=1:sizepop
%个体最优更新
if fitness(j) < fitnesspbest(j)
pbest(j,:) = pop(j,:);
fitnesspbest(j) = fitness(j);
end
%群体最优更新
if fitness(j) < fitnessgbest
gbest = pop(j,:);
fitnessgbest = fitness(j);
end
end
yy(i)=fitnessgbest;
end
%% 结果分析
figure;
plot(yy)
title('最优个体适应度','fontsize',12);
xlabel('进化代数','fontsize',12);ylabel('适应度','fontsize',12);
结果：最优适应度=（3.979836+0+3.979836+2.984877+3.979836+0+5.969754+5.969754+3.979836+0）/10=3.0843729
结论分析：c1为“认知”部分，反映了粒子对自身经历或经验的记忆或回忆，代表粒子有向自身历史最佳位置逼近的趋势；c2为“社会”部分，反映了粒子间协同合作与知识共享的群体历史经验，代表粒子有向群体或领域历史最佳位置逼近的趋势。因为，刚开始迭代，社会经验不足，主要依靠自身经验来学习，所以，c1是线性递减，c2线性递增，两者互相配合来共同调整群体最优适应度。
③单独调整同步学习因子
此时w=0.8
主函数代码：
%% 清空环境
clc
clear
%% 参数初始化
%粒子群算法中的三个参数
%c1 = 1.49445;%加速因子
%c2 = 1.49445;
%ws=0.9;
%we=0.2;
w=0.8 %惯性权重
%cstart1=2.5;
%cend1=0.5;
%cstart2=0.5;
%cend2=2.5;
cmax=2.1;
cmin=0.8;
maxgen=1000; % 进化次s数
sizepop=200; %种群规模
Vmax=1; %限制速度围
Vmin=-1;
popmax=5; %变量取值范围
popmin=-5;
dim=10; %适应度函数维数
func=1; %选择待优化的函数，1为Rastrigin,2为Schaffer,3为Griewank
Drawfunc(func);%画出待优化的函数，只画出二维情况作为可视化输出
%% 产生初始粒子和速度
for i=1:sizepop
%随机产生一个种群
pop(i,:)=popmax*rands(1,dim); %初始种群
V(i,:)=Vmax*rands(1,dim); %初始化速度
%计算适应度
fitness(i)=fun(pop(i,:),func); %粒子的适应度
end
%% 个体极值和群体极值
[bestfitness bestindex]=min(fitness);
gbest=pop(bestindex,:); %全局最佳
pbest=pop; %个体最佳
fitnesspbest=fitness; %个体最佳适应度值
fitnessgbest=bestfitness; %全局最佳适应度值
%% 迭代寻优
for i=1:maxgen
fprintf('第%d代，',i);
fprintf('最优适应度%f\n',fitnessgbest);
c1=cmax-(cmax-cmin)*(i/maxgen);
c2=c1;
%c1=cstart1+(cend1-cstart1)*(i/maxgen);
%c2=cstart2+(cend2-cstart2)*(i/maxgen);
%w=ws-(ws-we)/(i*maxgen);
for j=1:sizepop
%速度更新
V(j,:) = w*V(j,:) + c1*rand*(pbest(j,:) - pop(j,:)) + c2*rand*(gbest - pop(j,:)); %根据个体最优pbest和群体最优gbest计算下一时刻速度
V(j,find(V(j,:)>Vmax))=Vmax; %限制速度不能太大
V(j,find(V(j,:)<Vmin))=Vmin;
%种群更新
pop(j,:)=pop(j,:)+0.5*V(j,:); %位置更新
pop(j,find(pop(j,:)>popmax))=popmax;%坐标不能超出范围
pop(j,find(pop(j,:)<popmin))=popmin;
if rand>0.98 %加入变异种子，用于跳出局部最优值
pop(j,:)=rands(1,dim);
end
%更新第j个粒子的适应度值
fitness(j)=fun(pop(j,:),func);
end
for j=1:sizepop
%个体最优更新
if fitness(j) < fitnesspbest(j)
pbest(j,:) = pop(j,:);
fitnesspbest(j) = fitness(j);
end
%群体最优更新
if fitness(j) < fitnessgbest
gbest = pop(j,:);
fitnessgbest = fitness(j);
end
end
yy(i)=fitnessgbest;
end
%% 结果分析
figure;
plot(yy)
title('最优个体适应度','fontsize',12);
xlabel('进化代数','fontsize',12);ylabel('适应度','fontsize',12);
结果：最优适应度=（6.074091+0.995041+1.192545+2.276951+2.984877+1.784098+1.991193+1.989918+4.974795+5.969754）/10=3.0232903
结论分析：由上面的②解释可知 ，c1,c2都为学习因子，所以可以把它们置成相同的，代表群体自身学习和社会学习同步。
④线性调整权重和异步学习因子组合
代码：
%% 清空环境
clc
clear
%% 参数初始化
%粒子群算法中的三个参数
%c1 = 1.49445;%加速因子
%c2 = 1.49445;
ws=0.9;
we=0.2;
%w=0.8 %惯性权重
cstart1=2.5;
cend1=0.5;
cstart2=0.5;
cend2=2.5;
%cmax=2.1;
%cmin=0.8;
maxgen=1000; % 进化次s数
sizepop=200; %种群规模
Vmax=1; %限制速度围
Vmin=-1;
popmax=5; %变量取值范围
popmin=-5;
dim=10; %适应度函数维数
func=1; %选择待优化的函数，1为Rastrigin,2为Schaffer,3为Griewank
Drawfunc(func);%画出待优化的函数，只画出二维情况作为可视化输出
%% 产生初始粒子和速度
for i=1:sizepop
%随机产生一个种群
pop(i,:)=popmax*rands(1,dim); %初始种群
V(i,:)=Vmax*rands(1,dim); %初始化速度
%计算适应度
fitness(i)=fun(pop(i,:),func); %粒子的适应度
end
%% 个体极值和群体极值
[bestfitness bestindex]=min(fitness);
gbest=pop(bestindex,:); %全局最佳
pbest=pop; %个体最佳
fitnesspbest=fitness; %个体最佳适应度值
fitnessgbest=bestfitness; %全局最佳适应度值
%% 迭代寻优
for i=1:maxgen
fprintf('第%d代，',i);
fprintf('最优适应度%f\n',fitnessgbest);
%c1=cmax-(cmax-cmin)*(i/maxgen);
%c2=c1;
c1=cstart1+(cend1-cstart1)*(i/maxgen);
c2=cstart2+(cend2-cstart2)*(i/maxgen);
w=ws-(ws-we)/(i*maxgen);
for j=1:sizepop
%速度更新
V(j,:) = w*V(j,:) + c1*rand*(pbest(j,:) - pop(j,:)) + c2*rand*(gbest - pop(j,:)); %根据个体最优pbest和群体最优gbest计算下一时刻速度
V(j,find(V(j,:)>Vmax))=Vmax; %限制速度不能太大
V(j,find(V(j,:)<Vmin))=Vmin;
%种群更新
pop(j,:)=pop(j,:)+0.5*V(j,:); %位置更新
pop(j,find(pop(j,:)>popmax))=popmax;%坐标不能超出范围
pop(j,find(pop(j,:)<popmin))=popmin;
if rand>0.98 %加入变异种子，用于跳出局部最优值
pop(j,:)=rands(1,dim);
end
%更新第j个粒子的适应度值
fitness(j)=fun(pop(j,:),func);
end
for j=1:sizepop
%个体最优更新
if fitness(j) < fitnesspbest(j)
pbest(j,:) = pop(j,:);
fitnesspbest(j) = fitness(j);
end
%群体最优更新
if fitness(j) < fitnessgbest
gbest = pop(j,:);
fitnessgbest = fitness(j);
end
end
yy(i)=fitnessgbest;
end
%% 结果分析
figure;
plot(yy)
title('最优个体适应度','fontsize',12);
xlabel('进化代数','fontsize',12);ylabel('适应度','fontsize',12);
结果：最优适应度=（2.005228+2.984878+4.974796+0.994959+1.989950+4.974796+1.989918+0.000173+2.984883+4.974798）/10=2.7874379
结论分析：对比以上三组实验，可以发现这组实验结果的最优适应度是最小的，最靠近0的，说明这组实验结果较以上三组更准确。表明了，线性递减权重和异步调整学习因子两者共同变化，维持粒子对局部和全部寻找最优适应度的搜索平衡的效果最好。也就是说权重和学习因子都不能在迭代过程一成不变，而是随着迭代次数逐渐增加，往着有利于群体搜索群体最优适应的方向改变。
四.实验总结
通过此次关于粒子群算法的实验我发现了这个算法的一些特点：
（1）它是一类不确定算法。不确定性体现了自然界生物的生物机制，并且在求解某些特定问题方面优于确定性算法。
（2）是一类概率型的全局优化算法。非确定算法的优点在于算法能有更多机会求解全局最优解。
（3）不依赖于优化问题本身的严格数学性质。
（4）是一种基于多个智能体的仿生优化算法。粒子群算法中的各个智能体之间通过相互协作来更好的适应环境，表现出与环境交互的能力.
（5）具有自组织和进化性以及记忆功能，所有粒子都保存优解的相关知识。
（6）都具有稳健性。稳健性是指在不同条件和环境下算法的实用性和有效性，但是现在粒子群算法的数学理论基础还不够牢固，算法的收敛性还需要讨论。
当然粒子群算法最大的优点就在于简单容易实现并且没有许多参数的调节。目前已被广泛应用于函数优化、神经网络训练、模糊系统控制以及其他遗传算法的应用领域。