粒子群优化算法（PSO）
粒子群优化算法是一种进化计算技术，由Eberhart博士和Kennedy博士发明，源于对鸟群捕食的行为研究。PSO算法同遗传算法类似，是一种基于迭代的优化工具。系统初始化为一组随机解，通过迭代搜寻最优值。但是并没有遗传算法用的交叉以及变异，而是粒子在解空间追随最优的粒子进行搜索。
在PSO算法中，每个优化问题的解都是搜索空间中的一只鸟，被抽象为没有质量和体积的微粒，并将其延伸到N维空间。粒子 在
维空间里的位置表示为一个矢量，每个粒子的飞行速度也表示为一个矢量。
所有的粒子都有一个由被优化的函数决定的【适应值】（fittness），每个粒子还有一个速度，决定它们飞翔的力向和距离。粒子们知道【自己到目前为止发现的最好位置】（pbest）和现在的位置，这个可以看做是粒子自己的飞行经验。除此之外，每个粒子还知道【到目前为止整个群体中所有粒子发现的最好位置】（gbest，是pbest中的最好值），这个可以看做是粒子同伴的经验。粒子就是通过自己的经验和同伴中最好的经验来决定下一步的行动。
PSO算法首先初始化一群随机粒子（随机解），然后粒子们就追随当前的最优粒子在解空间中搜索，即通过迭代找到最优解。假设 维搜索空间中的第
个粒子的位置和速度分别为
和
，在每一次迭代中，粒子通过跟踪两个最优解来更新自己，第一个就是粒子本身所找到的最优解，即个体极值pbest，
；另一个是整个种群目前找到的最优解，即全局最优解gbest，
。在找到这两个最优值时，粒子根据如下公式来更新自己的速度和新的位置。
其中， 为惯性权重，
和
为正的学习因子，
和
为0到1之间均匀分布的随机数。
粒子群算法的性能很大程度上取决于算法的控制参数，例如粒子数、最大速度、学习因子、惯性权重等，各个参数的选取规则如下：
粒子数：粒子数的多少根据问题的复杂程度自行决定。对于一般的优化问题，取20-40个粒子就完全可以得到很好的结果；对于比较简单的问题，10个粒子已经足够可以取得好的结果；对于比较复杂的问题或者特定类别的问题，粒子数可以取到100以上。
粒子的维度：这是由优化问题决定，就是问题解的维度。
粒子的范围：由优化问题决定，每一维可设定不同的范围。
最大速度Vmax：决定粒子在一个循环中最大的移动距离，通常设定为粒子的范围宽度。
学习因子：学习因子使粒子具有自我总结和向群体中优秀个体学习的能力，从而向群体内或邻域内最优点靠近，通常取
和
为2，但也有其他的取值，一般
惯性权重：决定了对粒子当前速度继承的多少，合适的选择可以使粒子具有均衡的探索能力和开发能力，惯性权重的取法一般有常数法、线性递减法、自适应法等。
基本粒子群算法
原理
基本粒子群算法采用常数学习因子 和
及常惯性权重
，粒子根据如下的公式来更新自己的速度和新的位置。
算法步骤
基本粒子群算法的基本步骤如下：
【1】随机初始化种群中各微粒的位置和速度；
【2】评价每个微粒的适应度，将当前各微粒的位置和适应值存储在各微粒的pbest中，将所有pbest中适应值最优个体的位置和适应值存储于gbest中；
【3】用下式更新粒子的速度和位移：
【4】对每个微粒，将其适应值与其经历过的最好位置作比较，如果较好，则将其作为当前的最好位置；
【5】比较当前所有pbest和gbest的值，更新gbest；
【6】若满足停止条件（通常为预设的运算精度或迭代次数），搜索停止，输出结果，否则返回【3】继续搜索。
Matlab代码与试算
采用基本粒子群算法求取Sphere Model函数 的最小值。
Particle_Swarm_Optimization_Method.m
function [x_optimization,f_optimization]=Particle_Swarm_Optimization_Method(fitness,N,c1,c2,w,M,D)
format long;
% f：待优化的目标函数
% N：粒子数目
% c1：学习因子1
% c2：学习因子2
% w：惯性权重
% M：最大迭代次数
% D：自变量的个数
% x_optimization：目标函数取最小值时的自变量值
% f_optimization：目标函数的最小值
%------初始化种群的个体------------
for i = 1:N
for j = 1:D % D维搜索空间
x(i,j) = randn; %随机初始化位置 randn：正态分布的随机数
v(i,j) = randn; %随机初始化速度
end
end
%------先计算各个粒子的适应度，并初始化Pi和Pg----------------------
for i = 1:N
p(i) = fitness(x(i,:)); % 个体极值
pi(i,:) = x(i,:);
end
Pg = x(N,:); %Pg为全局最优
for i = 1:(N-1)
if fitness(x(i,:)) < fitness(Pg) % 将适应值最优的个体存储
Pg = x(i,:);
end
end
%------进入主要循环，按照公式依次迭代------------
for t = 1:M
for i = 1:N
v(i,:) = w*v(i,:) + c1*rand*(pi(i,:) - x(i,:)) + c2*rand*(Pg - x(i,:));
x(i,:) = x(i,:) + v(i,:);
if fitness(x(i,:)) < p(i) % 更新个体极值
p(i) = fitness(x(i,:));
pi(i,:) = x(i,:);
end
if p(i) < fitness(Pg) % 更新全局极值
Pg = pi(i,:);
end
end
Pbest(t)=fitness(Pg);
end
x_optimization = Pg';
f_optimization = fitness(Pg);
解：本例中函数最小点为 ，最小值为0，现在用PSO算法求最小值，首先看看不同迭代步数对结果的影响，粒子群规模为40，学习因子都为2，惯性权重取0.5，迭代步数分别取1000，5000，10000。
首先建立目标函数文件fitness.m。
fitness.m
function F = fitness(x)
F = 0;
for i = 1:30
F = F + x(i)^2;
end
test.m
[x_optimization,f_optimization] = Particle_Swarm_Optimization_Method(@fitness,40,2,2,0.5,1000,30); x_optimization = double(x_optimization); f_optimization = double(f_optimization); x_optimization f_optimization [x_optimization,f_optimization] = Particle_Swarm_Optimization_Method(@fitness,40,2,2,0.5,5000,30); x_optimization = double(x_optimization); f_optimization = double(f_optimization); x_optimization f_optimization [x_optimization,f_optimization] = Particle_Swarm_Optimization_Method(@fitness,40,2,2,0.5,10000,30); x_optimization = double(x_optimization); f_optimization = double(f_optimization); x_optimization f_optimization
命令行窗口
x_optimization = -0.203643806779367 0.004002658504391 -0.074941760623031 0.202828222081087 -0.000524266138323 -0.103282496122516 0.069940940915824 0.224712848311454 0.009116556529499 0.041747672665697 -0.115744530014599 0.008324066893113 -0.022343703252904 0.006402636234152 -0.231688976225837 0.105157069032689 0.019036253581153 -0.003145132913665 -0.035976401745438 -0.045138518377071 -0.077820436024370 -0.025756407633560 -0.159193629330254 0.278237342057476 -0.054883001208617 -0.231413816601646 0.166511326667179 -0.066557614071586 -0.030319312815608 -0.046405622613089 f_optimization = 0.439842896379014 x_optimization = 0.033911123724507 -0.008147918069209 0.017400189665077 -0.061229253488577 -0.123615521693440 -0.000936577773137 -0.096360645602730 -0.026363975956553 -0.106348809431333 0.006567647391543 -0.072915168393561 0.070642516969111 -0.031982133303148 0.052805456865892 -0.018763102557558 -0.049897180300623 0.051058695822665 0.099637809979345 0.038350127166670 -0.014412752178912 -0.084432915352984 -0.088508395148359 0.032420678214677 0.060300023200814 -0.003245103240019 0.051213192850383 -0.151210567142661 0.003972696340481 -0.046250488846935 0.159703206183376 f_optimization = 0.145864445774178 x_optimization = 0.004240164300054 -0.023527149367452 0.008715853267022 0.010757849454739 -0.129276052404614 0.024563617151546 0.012667589873705 -0.002475120619650 0.031828379323105 -0.019966653146839 -0.048641610586738 0.020712666845974 0.018611086412345 0.030191385786884 -0.004277974342346 0.012197935289109 -0.046489632116562 -0.018770500813800 -0.011646839137421 0.020356913351388 0.059918199922185 -0.080422636675985 -0.005354456376674 0.005028386980940 -0.040197783963151 -0.072035425652768 -0.010230701402834 0.017110069408012 0.012346401009822 -0.006123407845671 f_optimization = 0.044445115342608
在其他参数不变的情况下，一般迭代步数越大，求得的解的精度越高，但这并不是绝对的，因为PSO算法本质上也是一种随机算法，即使用同样的参数，每一次求解也可能得出不同的结果，同时如果对于多峰函数，PSO算法还有可能陷入局部最优点。
下面看粒子群规模对结果的影响，学习因子都为2，惯性权重为0.5，迭代步数都为10000，粒子群规模分别取50，60和80。
test.m
[x_optimization,f_optimization] = Particle_Swarm_Optimization_Method(@fitness,50,2,2,0.5,10000,30); x_optimization = double(x_optimization); f_optimization = double(f_optimization); x_optimization f_optimization [x_optimization,f_optimization] = Particle_Swarm_Optimization_Method(@fitness,60,2,2,0.5,10000,30); x_optimization = double(x_optimization); f_optimization = double(f_optimization); x_optimization f_optimization [x_optimization,f_optimization] = Particle_Swarm_Optimization_Method(@fitness,80,2,2,0.5,10000,30); x_optimization = double(x_optimization); f_optimization = double(f_optimization); x_optimization f_optimization
命令行窗口
x_optimization = -0.073611380522001 -0.021908541053983 -0.010810473749682 -0.042714310599363 -0.070682636079655 -0.017109091404758 0.049270845876185 0.031245077524827 0.083780249899494 -0.026065691289565 -0.023970672784154 -0.021291206278738 -0.091336531009053 0.000349365739912 0.016359094667683 -0.044224071171242 -0.070403279713934 -0.032001755957378 -0.024121811241049 -0.075618445900566 -0.061313274310489 -0.075184985641426 0.001548313670636 -0.010444997045968 0.012462497120872 0.002872591839675 0.082681063240601 -0.009230455495692 0.070214795579745 -0.070799535364695 f_optimization = 0.074656914017208 x_optimization = -0.094410028494356 -0.004589819770237 -0.007366364476676 0.076992652418444 -0.101145516026344 0.036864206183417 0.012220405060786 -0.039118880084193 -0.114145667260480 -0.124536379419096 0.026158906710145 0.078751071623844 -0.015334122007902 0.086005738304570 -0.017118781458174 0.019194143770836 0.017973083667948 0.042497351297172 -0.009329139701721 0.091028771109835 0.021201738445990 -0.061943723354998 -0.070485354240977 0.114457580587634 -0.052689011966822 0.083869616525665 -0.054902075056660 -0.138650587818845 -0.097705454809928 -0.251093283455119 f_optimization = 0.209403863870569 x_optimization = -0.003220432246500 -0.000076816207993 -0.005314899939919 0.014522092407788 -0.014696692331831 -0.007989286300044 0.001861681233914 -0.001638663900564 -0.005774603946167 -0.000464771967507 -0.006137818370138 -0.001398927980209 -0.004465238679096 0.005913867711140 -0.006720148483339 -0.000172177952597 -0.013951478158435 -0.003534508147095 -0.014222172034220 -0.004949109936742 -0.002065563091123 -0.005694020228356 0.016250856162974 -0.000572994144839 -0.004261331573720 -0.002280055059197 -0.015097904173360 0.001463793107481 0.002448043817841 -0.000206975356649 f_optimization = 0.001703273236163
粒子群规模不是越大越好，关键是各个参数之间的搭配，才能求得比较好的结果。
2.权重改进的粒子群算法
在微粒群算法的可调整参数中，惯性权重 是最重要的参数，较大的
有利于提高算法的全局搜索能力，而较小的
会增强算法的局部搜索能力，根据不同的权重变化公式，可得到不同的PSO算法，常见的有线性递减权重法、自适应权重法、随机权重法。
线性递减权重法
原理
由于较大的惯性因子有利于跳出局部极小点，便于全局搜索，而较小的惯性因子则有利于对当前的搜索区域进行精确局部搜索，以利于算法收敛，因此针对PSO算法容易早熟以及算法后期易在全局最优解附近产生振荡现象，可以采用线性变化的权重，让惯性权重从最大值 线性减小到最小值
，
随算法迭代次数的变化公式为：
其中， 、
分别表示
的最大值和最小值，
表示当前迭代步数，
表示最大迭代步数，通常取
。
算法步骤
线性递减粒子群算法的基本步骤如下：
【1】随机初始化种群中各微粒的位置和速度；
【2】评价每个微粒的适应度，将当前各微粒的位置和适应值存储在各微粒的pbest中，将所有pbest中适应值最优个体的位置和适应值存储于gbest中；
【3】用下式更新粒子的速度和位移：
【4】更新权重 ；
【5】对每个微粒，将其适应值与其经历过的最好位置作比较，如果较好，则将其作为当前的最好位置；比较当前所有pbest和gbest的值，更新gbest；
【6】若满足停止条件（通常为预设的运算精度或迭代次数），搜索停止，输出结果，否则返回【3】继续搜索。
Matlab代码与试算
用线性递减权重的粒子群算法求函数 的最小值。其中粒子数取40，学习因子都取2，最大权重取0.9，最小权重取0.4，迭代步数取10000。
解：此函数的最小点为 ，最小值为
。
建立目标函数文件fitness.m
fitness.m
function F = fitness(x) F = 100*(x(1)^2 - x(2))^2 + (1 - x(1))^2;
test.m
[x_optimization,f_optimization] = LDW_PSO(@fitness,40,2,2,0.9,0.4,10000,2); x_optimization f_optimization
LDW_PSO.m
function [x_optimization,f_optimization] = LDW_PSO(fitness,N,c1,c2,wmax,wmin,M,D) % f：待优化的目标函数 % N：粒子数目 % c1：学习因子1 % c2：学习因子2 % wmax：最大惯性权重 % wmin：最小惯性权重 % M：最大迭代次数 % D：自变量的个数 % x_optimization：目标函数取最小值时的自变量值 % f_optimization：目标函数的最小值 format long; %------初始化种群的个体------------ for i = 1:N for j = 1:D % D维搜索空间 x(i,j) = randn; %随机初始化位置 randn：正态分布的随机数 v(i,j) = randn; %随机初始化速度 end end %------先计算各个粒子的适应度，并初始化Pi和Pg---------------------- for i = 1:N p(i) = fitness(x(i,:)); % 个体极值 pi(i,:) = x(i,:); end Pg = x(N,:); %Pg为全局最优 for i = 1:(N-1) if fitness(x(i,:)) < fitness(Pg) % 将适应值最优的个体存储 Pg = x(i,:); end end %------进入主要循环，按照公式依次迭代------------ for t = 1:M for i = 1:N w = wmax - (t-1)*(wmax-wmin)/(M-1); % 权重线性递减 v(i,:) = w*(v(i,:) + c1*rand*(pi(i,:) - x(i,:)) + c2*rand*(Pg - x(i,:))); x(i,:) = x(i,:) + v(i,:); if fitness(x(i,:)) < p(i) % 更新个体极值 p(i) = fitness(x(i,:)); pi(i,:) = x(i,:); end if p(i) < fitness(Pg) % 更新全局极值 Pg = pi(i,:); end end end x_optimization = Pg'; f_optimization = fitness(Pg);
命令行窗口
x_optimization = 1 1 f_optimization = 0
对于本例题中的函数而言，用线性递减权重的粒子群算法求得了非常精确的最优点。但是在实际问题中，对于不同问题，其每次迭代所需的比例关系并不相同，所以 的线性递减只对某些问题很有效。
此外，如果在进化初期搜索不到最优点，随着 的逐渐减小，算法局部收敛能力加强，容易陷入局部最优；如果在进化初期探测到次好点，这时
的相对取小就可使算法很快搜索到最优点，而
的线性递减降低了算法的收敛速度。
自适应权重法
原理
为了平衡PSO算法的全局搜索能力和局部改良能力，还可采用非线性的动态惯性权重系数公式，其表达式如下：
其中 分别表示
的最大值和最小值，
表示粒子当前的目标函数值，
和
分别表示当前所有微粒的平均目标值和最小目标值。在上式中，惯性权重随着微粒的目标函数值而自动改变，因此称为自适应权重。
当各微粒的目标值趋于一致或者趋于局部最优时，将使惯性权重增加，而各微粒的目标值比较分散时，将使惯性权重减小，同时对于目标函数值优于平均目标值的微粒，其对应的惯性权重因子较小，从而保护了该微粒，反之对于目标函数值差于平均目标值的微粒，其对应的惯性权重因子较大，使得该微粒向较好的搜索区域靠拢。
算法步骤
自适应权重粒子群算法的基本步骤如下：
【1】随机初始化种群中各微粒的位置和速度；
【2】评价每个微粒的适应度，将当前各微粒的位置和适应值存储在各微粒的pbest中，将所有pbest中适应值最优个体的位置和适应值存储于gbest中；
【3】用下式更新粒子的速度和位移：
【4】更新权重
【5】对每个微粒，将其适应值与其经历过的最好位置作比较，如果较好，则将其作为当前的最好位置；比较当前所有pbest和gbest的值，更新gbest；
【6】若满足停止条件（通常为预设的运算精度或迭代次数），搜索停止，输出结果，否则返回【3】继续搜索。
Matlab代码与试算
用自适应权重的粒子群算法求函数 的最小值。取粒子数为40，学习因子都取2，最大惯性权重取0.9，最小惯性权重取0.6，迭代步数取10000。
解：建立fitness.m文件
fitness.m
function F = fitness(x) F = 0.5 + (sin(x(1)^2 + x(2)^2)^2 - 0.5)/(1.0 + 0.001*(x(1)^2 + x(2)^2))^2;
test.m
[x_optimization,f_optimization] = AW_PSO(@fitness,40,2,2,0.9,0.6,10000,2); x_optimization f_optimization
AW_PSO.m
function [x_optimization,f_optimization] = AW_PSO(fitness,N,c1,c2,wmax,wmin,M,D) % f：待优化的目标函数 % N：粒子数目 % c1：学习因子1 % c2：学习因子2 % wmax：最大惯性权重 % wmin：最小惯性权重 % M：最大迭代次数 % D：自变量的个数 % x_optimization：目标函数取最小值时的自变量值 % f_optimization：目标函数的最小值 format long; %------初始化种群的个体------------ for i = 1:N for j = 1:D % D维搜索空间 x(i,j) = randn; %随机初始化位置 randn：正态分布的随机数 v(i,j) = randn; %随机初始化速度 end end %------先计算各个粒子的适应度，并初始化Pi和Pg---------------------- for i = 1:N p(i) = fitness(x(i,:)); % 个体极值 pi(i,:) = x(i,:); end Pg = x(N,:); %Pg为全局最优 for i = 1:(N-1) if fitness(x(i,:)) < fitness(Pg) % 将适应值最优的个体存储 Pg = x(i,:); end end %------进入主要循环，按照公式依次迭代------------ for t = 1:M for j = 1:N f_optimization(j) = fitness(x(j,:)); end f_avg = sum(f_optimization)/N; % 适应度平均值 f_min = min(f_optimization); % 适应度最小值 for i = 1:N if f_optimization(i) <= f_avg w = wmin - (f_optimization(i) - f_min)*(wmax - wmin)/(f_avg - f_min); else w = wmax; end v(i,:) = w*(v(i,:) + c1*rand*(pi(i,:) - x(i,:)) + c2*rand*(Pg - x(i,:))); x(i,:) = x(i,:) + v(i,:); if fitness(x(i,:)) < p(i) % 更新个体极值 p(i) = fitness(x(i,:)); pi(i,:) = x(i,:); end if p(i) < fitness(Pg) % 更新全局极值 Pg = pi(i,:); end end end x_optimization = Pg'; f_optimization = fitness(Pg);
命令行窗口
x_optimization = 1.0e-06 * 0.172784775909731 -0.008767853535804 f_optimization = 0
本例题中的函数的最小点为 ，最小值为0，自适应权重的粒子群算法求得的结果的精度还是可以的。
随机权重法
原理
将PSO算法中设定 为服从某种随机分布的随机数，这样一定程度上可从两方面来克服
的线性递减所带来的不足。
首先，如果在进化初期接近最好点，随机 可能产生相对小的
值，加快算法的收敛速度，另外，如果在算法初期找不到最好点，
的线性递减，使得算法最终收敛不到此最好点，而
的随机生成可以克服这种局限。
计算公式如下：
其中 表示标准正态分布的随机数，
表示0到1之间的随机数。
算法步骤
随机权重粒子群算法的基本步骤如下：
【1】随机初始化种群中各微粒的位置和速度；
【2】评价每个微粒的适应度，将当前各微粒的位置和适应值存储在各微粒的pbest中，将所有pbest中适应值最优个体的位置和适应值存储于gbest中；
【3】用下式更新粒子的速度和位移：
【4】更新权重
【5】对每个微粒，将其适应值与其经历过的最好位置作比较，如果较好，则将其作为当前的最好位置；比较当前所有pbest和gbest的值，更新gbest；
【6】若满足停止条件（通常为预设的运算精度或迭代次数），搜索停止，输出结果，否则返回【3】继续搜索。
Matlab代码与试算
求下面函数的最小值
取粒子数为40，学习因子都取2，随机权重平均值的最大值取0.8，随机权重平均值的最小值取0.5，随机权重平均值的方差取0.2，迭代步数取10000。
解：首先建立目标函数文件fitness.m文件
fitness.m
function F = fitness(x) F = 4*x(1)^2 - 2.1*x(1)^4 + x(1)^6/3 + x(1)*x(2) - 4*x(2)^2 + 4*x(2)^4;
test.m
[x_optimization,f_optimization] = RW_PSO(@fitness,40,2,2,0.8,0.5,0.2,10000,2); x_optimization f_optimization
RW_PSO.m
function [x_optimization,f_optimization] = RW_PSO(fitness,N,c1,c2,mu_max,mu_min,sigma,M,D) % f：待优化的目标函数 % N：粒子数目 % c1：学习因子1 % c2：学习因子2 % mu_max：随机权重平均值最大值 % mu_min：随机权重平均值最小值 % sigma：随机权重的方差 % M：最大迭代次数 % D：自变量的个数 % x_optimization：目标函数取最小值时的自变量值 % f_optimization：目标函数的最小值 format long; %------初始化种群的个体------------ for i = 1:N for j = 1:D % D维搜索空间 x(i,j) = randn; %随机初始化位置 randn：正态分布的随机数 v(i,j) = randn; %随机初始化速度 end end %------先计算各个粒子的适应度，并初始化Pi和Pg---------------------- for i = 1:N p(i) = fitness(x(i,:)); % 个体极值 pi(i,:) = x(i,:); end Pg = x(N,:); %Pg为全局最优 for i = 1:(N-1) if fitness(x(i,:)) < fitness(Pg) % 将适应值最优的个体存储 Pg = x(i,:); end end %------进入主要循环，按照公式依次迭代------------ for t = 1:M for i = 1:N mu = mu_min + (mu_max - mu_min)*rand(); % 随机权重的平均值 w = mu + sigma*randn(); % 随机权重 v(i,:) = w*(v(i,:) + c1*rand*(pi(i,:) - x(i,:)) + c2*rand*(Pg - x(i,:))); x(i,:) = x(i,:) + v(i,:); if fitness(x(i,:)) < p(i) % 更新个体极值 p(i) = fitness(x(i,:)); pi(i,:) = x(i,:); end if p(i) < fitness(Pg) % 更新全局极值 Pg = pi(i,:); end end end x_optimization = Pg'; f_optimization = fitness(Pg);
命令行窗口
x_optimization = -0.089842011145412 0.712656402204421 f_optimization = -1.031628453489878
本题中的函数的理论最小点有两个，分别为 和
，最小值为
，随机权重的粒子群算法求得的结果精度还可以。
