1 简介
1.1 关于速度和位置
粒子群算法通过设计一种无质量的粒子来模拟鸟群中的鸟，粒子仅具有两个属性：速度和位置，速度代表移动的快慢，位置代表移动的方向。
鸟被抽象为没有质量和体积的微粒(点)，并延伸到N维空间，粒子i在N维空间的位置表示为矢量Xi＝(x1，x2，…，xN)，飞行速度表示为矢量Vi＝(v1，v2，…，vN)。每个粒子都有一个由目标函数决定的适应值(fitness value)，并且知道自己到目前为止发现的最好位置(pbest)和现在的位置Xi。这个可以看作是粒子自己的飞行经验。除此之外，每个粒子还知道到目前为止整个群体中所有粒子发现的最好位置(gbest)(gbest是pbest中的最好值)，这个可以看作是粒子同伴的经验。粒子就是通过自己的经验和同伴中最好的经验来决定下一步的运动。
2.2 速度和位置的更新
PSO初始化为一群随机粒子(随机解)。然后通过迭代找到最优解。在每一次的迭代中，粒子通过跟踪两个“极值”(pbest，gbest)来更新自己。在找到这两个最优值后，粒子通过下面的公式来更新自己的速度和位置。
对于公式(1)：
公式(1)的第①部分称为【记忆项】，表示上次速度大小和方向的影响；
公式(1)的第②部分称为【自身认知项】，是从当前点指向粒子自身最好点的一个矢量，表示粒子的动作来源于自己经验的部分；
公式(1)的第③部分称为【群体认知项】，是一个从当前点指向种群最好点的矢量，反映了粒子间的协同合作和知识共享。粒子就是通过自己的经验和同伴中最好的经验来决定下一步的运动。
以上面两个公式为基础，再来看一个公式：
公式(2)和 公式(3)被视为标准PSO算法。
1.3 标准PSO算法的流程
1）初始化一群微粒(群体规模为N)，包括随机位置和速度；
2）评价每个微粒的适应度；
3）对每个微粒，将其适应值与其经过的最好位置pbest作比较，如果较好，则将其作为当前的最好位置pbest；
4）对每个微粒，将其适应值与其经过的最好位置gbest作比较，如果较好，则将其作为当前的最好位置gbest；
5）根据公式(2)、(3)调整微粒速度和位置；
6）未达到结束条件则转第2）步。
迭代终止条件根据具体问题一般选为最大迭代次数Gk或(和)微粒群迄今为止搜索到的最优位置满足预定最小适应阈值。
PSO流程图解
2 部分代码
%-------------------------------------------------------------------------%
% P A R T I C L E S W A R M O P T I M I Z A T I O N %
%-------------------------------------------------------------------------%
% Function "PSO" runs a Particle Optimization Algorithm (PSO) for optimi- %
% zing a given 2D function. %
% %
% Input parameters: %
% * param: Struct with pop_s, w, C1, C2, type, max_gen, xy_lim, %
% n and err parameters. %
% * h_cont: Handle of the contour graph. %
% * h_fit: Handle of the fitness variation graph. %
% %
% Output variables: %
% * pos: Optimized position. %
% * f: Optimized fitness. %
% * g: Last computed generation. %
%-------------------------------------------------------------------------%
% - Author: V韈tor Mart韓ez Cagigal %
% - Date: 05/12/2015 %
% - Version: 1.0 %
% - E-mail: vicmarcag (dot) gmail (dot) com %
%-------------------------------------------------------------------------%
function [pos, f, g] = pso(fun, param, h_cont, h_fit)
% Storing parameters
pop_s = param.pop_s; % Population size
w = param.w; % Inertia constant (
C1 = param.C1; % Individual confidence factor
C2 = param.C2; % Swarm confidence factor
type = param.type; % Minimization or maximization?
max_gen = param.max_gen; % Maximum generation limit
xy_lim = param.xy_lim; % Recommended X and Y limits
n = param.n; % Recommended number of lines to plot
err = param.err; % Precision (error tolerance)
% Parameter correction: percentage of the region of interest limits
w = abs((w*xy_lim(2))/100);
C1 = abs((C1*xy_lim(2))/100);
C2 = abs((C2*xy_lim(2))/100);
% Initialization: S contains the whole swarm
S = NaN(pop_s, 4, 2);
% Position initialization in S(:,1,[x y])
S(:,1,1) = xy_lim(1) + (xy_lim(2)-xy_lim(1)).*rand(pop_s,1);
S(:,1,2) = xy_lim(1) + (xy_lim(2)-xy_lim(1)).*rand(pop_s,1);
% Velocity initialization in S(:,2,[vx,vy])
S(:,2,1:2) = 0.10.*S(:,1,1:2);
% Best particle position so far init in S(:,3,[best_x,best_y])
S(:,3,1:2) = S(:,1,1:2);
% Best particle fitness value so far init in S(:,4,best_f)
S(:,4,1) = fun(S(:,1,1), S(:,1,2));
% Best global value and position ever found initialization
if(type == 1), [gBest(3),index] = min(S(:,4,1)); end
if(type == 2), [gBest(3),index] = max(S(:,4,1)); end
gBest(1:2) = S(index,3,1:2);
% Plotting the swarm
% Particles position
hold(h_cont,'on');
h_par = plot(h_cont, S(:,1,1), S(:,1,2),'.k');
axis(h_cont, [xy_lim(1) xy_lim(2) xy_lim(1) xy_lim(2)]);
hold(h_cont,'on');
% Best global and generation fitness init
fGlo = []; fGen = []; h_bes = [];
% Generations' iteration loop
tol = err * 100;
gen = 0;
while ((gen <= max_gen) && (tol > err))
% New generation
gen = gen + 1;
% Updating the particles velocity
S(:,2,1:2) = w.*squeeze(S(:,2,1:2)) + ...
(C1.*repmat(rand(pop_s,1),1,2).*(squeeze(S(:,3,1:2)-S(:,1,1:2)))) + ...
(C2.*repmat(rand(pop_s,1),1,2).*(repmat(gBest(1:2),pop_s,1)-squeeze(S(:,1,1:2))));
% Updating the particles position
S(:,1,1:2) = S(:,1,1:2) + S(:,2,1:2);
% Evaluating the fitness of each particle
f = fun(squeeze(S(:,1,1)), squeeze(S(:,1,2)));
if(type == 1) % Minimization
S(f<S(:,4,1),3,1:2) = S(f<S(:,4,1),1,1:2); % Updating p_best position
S(f<S(:,4,1),4,1) = f(f<S(:,4,1)); % Updating p_best fitness
[genBest,index] = min(S(:,4,1)); % Best generation global fitness
if(genBest < gBest(3))
tol = gBest(3) - genBest; % Tolerance
gBest(3) = genBest; % Updating g_best fitness
gBest(1:2) = S(index,3,1:2); % Updating g_best position
end
end
if(type == 2) % Maximization
S(f>S(:,4,1),3,1:2) = S(f>S(:,4,1),1,1:2); % Updating p_best position
S(f>S(:,4,1),4,1) = f(f>S(:,4,1)); % Updating p_best fitness
[genBest,index] = max(S(:,4,1)); % Best generation global fitness
if(genBest > gBest(3))
tol = genBest - gBest(3); % Tolerance
prevBest = gBest(3); % Previous g_best for tolerance
gBest(3) = genBest; % Updating g_best fitness
gBest(1:2) = S(index,3,1:2); % Updating g_best position
end
end
% Plotting the swarm
% Particles position
delete(h_par); delete(h_bes)
h_par = plot(h_cont, S(:,1,1), S(:,1,2),'.k'); hold(h_cont,'on');
h_bes = plot(h_cont, gBest(1), gBest(2),'+r'); hold(h_cont,'on');
axis(h_cont, [xy_lim(1) xy_lim(2) xy_lim(1) xy_lim(2)]);
% Best global and generation fitness
if(type == 1), fGlo(length(fGlo)+1) = 1/gBest(3); else fGlo(length(fGlo)+1) = gBest(3); end
if(type == 1), fGen(length(fGen)+1) = 1/min(f); else fGen(length(fGen)+1) = min(f); end
plot(h_fit, 1:gen, fGlo, 'm', 'LineWidth', 2); hold(h_fit,'on');
plot(h_fit, 1:gen, fGen, 'k'); hold(h_fit,'on');
xlim(h_fit, [1 max_gen]); xlabel(h_fit, 'Generations','FontSize',9);
title(h_fit,['Generation : ' num2str(gen) ' - Fitness: ' num2str(gBest(3))]);
legend(h_fit,'f_{global}','f_{gen}','Location','SouthEast');
if(type == 1), ylabel(h_fit, '1/f','FontSize',9); else ylabel(h_fit, 'f','FontSize',9); end
set(gca,'FontSize',10);
drawnow;
end
hold(h_fit,'off'); hold(h_cont,'off');
% Returning the sub-optimal minima/maxima
pos = gBest(1:2);
f = gBest(3);
g = gen;
end
3 仿真结果
4 参考文献
[1]李浩. 基于粒子群优化算法的车间调度系统的研究与设计[D]. 宁夏大学, 2018.
博主简介：擅长智能优化算法、神经网络预测、信号处理、元胞自动机、图像处理、路径规划、无人机等多种领域的Matlab仿真，相关matlab代码问题可私信交流。
部分理论引用网络文献，若有侵权联系博主删除。
