粒子群算法基本原理

浅说一下思路：

一群鸟找吃的，假设有几只鸟同时都找到了好吃的，就判断谁找到的更好，相当于就是当前的全局最优，全体的鸟都考虑想这边靠拢。但是其他几只找到的鸟，也会把自己找到的好吃的作为考虑，选择最近的吃的，这个属于是局部最优。总之，就是一群鸟互相合作，使得利益最大化，找到全局最优

基本原理

先列公式

假设在一个D维的目标搜索空间中，有N个粒子组成一个群落，其中第i个例子表示为一个D维的向量：

$$X_i = (x_{i1} , x_{i2}, \cdots , x_{iD}), i = 1, 2, \cdots ,N$$

第i个粒子的“飞行”速度也是一个D维第向量，记为：

$$V_i = (v_{i1} , v_{i2} , \cdots ,v_{iD}) , i = 1, 2, \cdots , N$$

在第 t 代的第 i 个粒子向第 t+1 代进化时，根据如下式子更新：

$$v_{ij}(t+1) = wv_{ij}(t) + c_1r_1(t)[p_{ij}(t)-x_{ij}(t)] + c_2r_2(t)[p_{gj}(t)-x_{ij}(t)]$$ $$x_{ij}(t+1) = x_{ij}(t) + v_{ij}(t+1)$$

公式解释

第一个公式其实就是和坐标一样。假设是2维平面空间，$x_i = (1,2)$就是在平面直角坐标系的$x=1, y = 2$的点上嘛。多加几维也就是表示各个维度的点的位置了。所以就是N个粒子组成的群体，每个粒子都有自己的一个坐标位置。

第二个公式也就是利用向量表示速度。速度有大小和方向，这个其实也就是物理知识。同样的放到二维来理解，之后再提升维度。

第三个公式是迭代下一代的方向。公式主要是三部分：$wv_{ij}(t)$ 、 $c_1r_1(t)[p_{ij}(t)-x_{ij}(t)]$ 、 $c_2r_2(t)[p_{gj}(t)-x_{ij}(t)]$

• $wv_{ij}(t)$ ：表示自身的一个惯性，也就是保持原有运动的趋势。 $v_{ij}(t)$ 就是上一代第 $i$ 个粒子的方向向量的第 $j$ 列。其实就是上一代的方向。$w$是惯性权重。最终理解就是，每个粒子对于保持原有运动方向的执着程度。其中 $w$ 是全局参数。
• $c_1r_1(t)[p_{ij}(t)-x_{ij}(t)]$ ：表示自己发现的食物对自己的吸引，也就是有向局部最优解靠拢的趋势。 $c_i$ 是全局参数，是自我学习因子，或者理解为自我发现最优解的权重，权重越大表示我越相信这个粒子可以独立找到最优解。 $r_i(t)$ 是一个随机数，用来模拟食物的好坏。 $p_{ij}(t)$ 表示局部最优所在的坐标位置。 $x_{ij}(t)$ 就是自己当前的位置。 $[p_{ij}(t)-x_{ij}(t)]$ 是局部最优和自己当前位置之差，也就决定了靠近局部最优的方向和速度。
• $c_2r_2(t)[p_{gj}(t)-x_{ij}(t)]$ ：这个和上一个式子很接近。 $c_2$ 也是全局参数，表示群体学习因子，权重越大表示我越相信这个群体当前找到的全局最优就是最终真正的最优解。 $p_{gj}(t)$ 就是当前这一代群体找到的全局最优解的坐标。

三个部分讲解完毕。合起来就是考虑三个向量矢量相加为最终方向即可。这里基本也就可以理解各个权重的作用了。主要就是 $w$ 、 $c_1$ 、 $c_2$ 三个权重。

第四个公式就是位置迭代的公式了。就是当前位置和下一代移动的方向矢量相加，结果就是该粒子下一代的位置。

算法流程

全局参数

全局参数的选取取决于优化函数。

• N：群体粒子个数

• D：可行解的维数

• $w$：惯性权重

• $c_1$：自我学习因子
• $c_2$：群体学习因子
• limit：位置的限制，别跑出地图了
• vlimit：速度的限制，别飞得太猛，直接闪现到另一个点，错过最优解

流程

• Step1：初始化种群的位置 $X$

  x = limit(1) + (limit(2) -  limit(1)) .* rand(N, d);

• Step2：计算个体适应度

  f = @(x) 函数表达式
  f(x)

• Step3：更新粒子速度和位置

  要记得检查边界。

  v = w*v + c1*rand*(xm-x) + c2*rand*(repmat(ym, N, 1)-x);  % 更新速度
  % 边界速度处理
  v(v > vlimit(2)) = vlimit(2);                             % 大于最大速度时取最大速度
  v(v < vlimit(1)) = vlimit(1);                             % 小于最大速度时取最小速度
  x = x + v;                                                % 位置更新
  % 边界位置处理  
  x(x > limit(2)) = limit(2);
  x(x < limit(1)) = limit(1);

• Step4：计算新位置的适应度，新位置适应度更高才移动向该点移动，否则不移动。（属于是固守当前的最优解，不贪便宜，比较稳）

  for i = 1 : N
      if fx(i)  < fxm(i)
          fxm(i) = fx(i);     % 更新个体历史最佳适应度  
          xm(i,:) = x(i,:);   % 更新个体历史最佳位置(取值)  
      end
  end

• Step5：继续执行Step2 - Step4直到满足终止条件，也就是满足迭代次数。或者有 k 代都没让适应度变化。

算法评价

• 优点

  原理简单，实现容易，参数较少，需要条件少

• 缺点

  很容易过早收敛到局部最优，而无法达到全局最优

问题案例

一维案例

这里的可行解就是x的取值。绘制的图像是二维图像。

求解$f = x \sin x \cos(2x) - 2x\sin(3x) + 3x*\sin(4x)$的最大值和最小值

最大值

%% 初始化种群  
clc, clear, close all;
f= @(x) x .* sin(x) .* cos(2 * x) - 2 * x .* sin(3 * x) +3 * x .* sin(4 * x); % 求这个函数的最大值

N = 20;                         % 初始种群个数  
d = 1;                          % 可行解维数  
G = 100;                        % 迭代次数       
limit = [0, 50];                % 设置位置参数限制  
vlimit = [-10, 10];             % 设置速度限制  
w = 0.8;                        % 惯性权重  
c1 = 0.5;                       % 自我学习因子  
c2 = 0.5;                       % 群体学习因子 

figure(1); ezplot(f, [0, 0.01, limit(2)]);   % 曲线

x = limit(1) + (limit(2)-limit(1)) .* rand(N, d); % 初始种群的位置  

v = rand(N, d);                  % 初始种群的速度
xm = x;                          % 每个个体的历史最佳位置
ym = zeros(1, d);                % 种群的历史最佳位置
fxm = ones(N, 1);                % 每个个体的历史最佳适应度
fym = 0;                         % 种群历史最佳适应度
hold on;
plot(xm, f(xm), 'ro'); title('初始状态图');
figure(2)
record = zeros(G, 1);            % 记录器
for t = 1 : G
    fx = f(x);                   % 个体当前适应度
    for i = 1 : N
        if fx(i) > fxm(i)
            fxm(i) = fx(i);      % 更新个体历史最佳适应度
            xm(i,:) = x(i,:);    % 更新个体历史最佳位置(取值)
        end
    end
    if  min(fxm) > fym 
        [fym, nmin] = max(fxm);  % 更新群体历史最佳适应度
        ym = xm(nmin, :);        % 更新群体历史最佳位置
    end
    v = w*v + c1*rand*(xm-x) + c2*rand*(repmat(ym, N, 1)-x);% 速度更新
    % 边界速度处理
    v(v > vlimit(2)) = vlimit(2);
    v(v < vlimit(1)) = vlimit(1);
    x = x + v; % 位置更新
    % 边界位置处理
    x(x > limit(2)) = limit(2);
    x(x < limit(1)) = limit(1);
    
    record(t) = fym; % 最大值记录
    x0 = 0 : 0.01 : limit(2);
    subplot(1,2,1); plot(x0, f(x0), 'b-', x, f(x), 'ro'); title('状态位置变化')
    subplot(1,2,2); plot(record(1:t)); title('最优适应度进化过程')
    pause(0.01) 
end

x0 = 0 : 0.01 : limit(2);
figure(4); plot(x0, f(x0), 'b-', x, f(x), 'ro'); title('最终状态位置')
disp(['最优值：',num2str(fym)]);
disp(['变量取值：',num2str(ym)]);

最小值

%% 初始化种群  
clc, clear, close all;
f= @(x) x .* sin(x) .* cos(2 * x) - 2 * x .* sin(3 * x) +3 * x .* sin(4 * x); % 求这个函数的最小值

N = 20;                         % 初始种群个数
d = 1;                          % 可行解维数
G = 100;                        % 迭代次数
limit = [0, 50];                % 设置位置参数限制
vlimit = [-10, 10];             % 设置速度限制
w = 0.8;                        % 惯性权重
c1 = 0.5;                       % 自我学习因子
c2 = 0.5;                       % 群体学习因

figure(1);ezplot(f,[0,0.01,limit(2)]);   % 曲线

x = limit(1) + (limit(2) -  limit(1)) .* rand(N, d); % 初始种群的位置

v = rand(N, d);                  % 初始种群的速度
xm = x;                          % 每个个体的历史最佳位置
ym = zeros(1, d);                % 种群的历史最佳位置
fxm = ones(N, 1)*inf;            % 每个个体的历史最佳适应度
fym = inf;                       % 种群历史最佳适应度
hold on;
plot(xm, f(xm), 'ro');title('初始状态图');
figure(2)
record = zeros(G, 1);          % 记录器
for t = 1 : G
    fx = f(x) ;                % 个体当前适应度
    for i = 1 : N
        if fx(i) < fxm(i)
            fxm(i) = fx(i);     % 更新个体历史最佳适应度
            xm(i,:) = x(i,:);   % 更新个体历史最佳位置(取值)
        end
    end
    if  min(fxm) < fym
        [fym, nmin] = min(fxm);   % 更新群体历史最佳适应度
        ym = xm(nmin, :);      % 更新群体历史最佳位置
    end
    v = w*v + c1*rand*(xm-x) + c2*rand*(repmat(ym, N, 1)-x);% 速度更新
    % 边界速度处理
    v(v > vlimit(2)) = vlimit(2);
    v(v < vlimit(1)) = vlimit(1);
    x = x + v; % 位置更新
    % 边界位置处理
    x(x > limit(2)) = limit(2);  
    x(x < limit(1)) = limit(1);  
    
    record(t) = fym; % 最大值记录
    x0 = 0 : 0.01 : limit(2);
    subplot(1,2,1) plot(x0, f(x0), 'b-', x, f(x), 'ro'); title('状态位置变化')
    subplot(1,2,2); plot(record(1:t)); title('最优适应度进化过程')
    pause(0.01)
end

x0 = 0 : 0.01 : limit(2);
figure(4); plot(x0, f(x0), 'b-', x, f(x), 'ro'); title('最终状态位置')
disp(['最优值：',num2str(fym)]);
disp(['变量取值：',num2str(ym)]);

可以发现其他地方都没有改，主要就是改变了两个符号，这个在实际问题中可以讨论。就是下面代码都大于符号换成小于符号

for i = 1 : N
    if fx(i) < fxm(i)
        fxm(i) = fx(i);     % 更新个体历史最佳适应度
        xm(i,:) = x(i,:);   % 更新个体历史最佳位置(取值)
    end
end
if  min(fxm) < fym
    [fym, nmin] = min(fxm);   % 更新群体历史最佳适应度
    ym = xm(nmin, :);      % 更新群体历史最佳位置
end

注意：如果值不是理想值，即道理局部最优解，可以设置粒子个数更大，达到全局最优解。

二维案例

可行解是平面坐标$(x, y)$ 的取值。绘制的图像是三维图像，图像比较向鸡蛋托盘。

求解：$f = 20 + x^2 + y^2 - 10 \cos(2 \pi x) - 10\cos(2\pi y)， x, y \in [-5.12, 5.12]$的最小值

最小值

%% 初始化种群  
clc, clear, close all;
f = @(x,y) 20 +  x.^2 + y.^2 - 10*cos(2*pi.*x) - 10*cos(2*pi.*y);
 
x0 = -5.12 : 0.05 : 5.12; % x的边界
y0 = x0;                  % y的边界
[X,Y] = meshgrid(x0, y0); % 绘制边界图
Z = f(X,Y);
figure(1); mesh(X,Y,Z);
colormap(parula(5));

N = 50;                         % 初始种群个数  
d = 2;                          % 可行解维数  
G = 100;                        % 迭代次数
limit = [-5.12,5.12];           % 设置位置参数限制
vlimit = [-.5, .5];             % 设置速度限制
w = 0.8;                        % 惯性权重
c1 = 0.5;                       % 自我学习因子
c2 = 0.5;                       % 群体学习因子

x = limit(1) + (limit(2) -  limit(1)) .* rand(N, d);%初始种群的位置  

v = rand(N, d);                  % 初始种群的速度
xm = x;                          % 每个个体的历史最佳位置
ym = zeros(1, d);                % 种群的历史最佳位置
fxm = ones(N, 1)*inf;            % 每个个体的历史最佳适应度
fym = inf;                       % 种群历史最佳适应度
hold on 
scatter3(x(:,1), x(:,2), f(x(:,1),x(:,2)), 'r*');
figure(2)

% 群体更新  
record = zeros(G,1);            % 记录器
for t = 1 : G
     fx = f(x(:,1), x(:,2));    % 个体当前适应度
     for i = 1 : N
        if  fx(i) < fxm(i)
            fxm(i) = fx(i);     % 更新个体历史最佳适应度
            xm(i,:) = x(i,:);   % 更新个体历史最佳位置(取值)
        end
     end
    if min(fxm) < fym
        [fym, nmin] = min(fxm); % 更新群体历史最佳适应度
        ym = xm(nmin, :);       % 更新群体历史最佳位置
    end
    v = v * w + c1 * rand * (xm - x) + c2 * rand * (repmat(ym, N, 1) - x); % 速度更新
    % 边界速度处理
    v(v > vlimit(2)) = vlimit(2);
    v(v < vlimit(1)) = vlimit(1);
    x = x + v; % 位置更新
    % 边界位置处理
    x(x > limit(2)) = limit(2);
    x(x < limit(1)) = limit(1);
    record(t) = fym; % 最优值记录
    subplot(1, 2, 1);
    mesh(X, Y, Z);
    hold on;
    scatter3(x(:,1), x(:,2), f(x(:,1),x(:,2)), 'r*'); title(['状态位置变化', '-迭代次数：', num2str(t)])
    subplot(1,2,2); plot(record); title('最优适应度进化过程')
    pause(0.01)
end

figure(4);mesh(X,Y,Z); hold on
scatter3(x(:,1), x(:,2), f(x(:,1), x(:,2)), 'r*'); title('最终状态位置')
disp(['最优值：', num2str(fym)]);
disp(['变量取值：', num2str(ym)]);

想要改为最大值其实就是改是否移动部分的代码

for i = 1 : N
    if  fx(i) < fxm(i)
        fxm(i) = fx(i);     % 更新个体历史最佳适应度
        xm(i,:) = x(i,:);   % 更新个体历史最佳位置(取值)
    end
end
if min(fxm) < fym
    [fym, nmin] = min(fxm); % 更新群体历史最佳适应度
    ym = xm(nmin, :);       % 更新群体历史最佳位置
end

优化求解的精确度：粒子群算法优化方法