粒子群算法笔记

介绍：

粒子群是智能算法的一种，也称为鸟群算法，简单来说就是鸟根据惯性、自身经验、社会经验三者调整自己的每一步路线，最后找到最优解的方法。说白了，粒子群算法的实质依然是蒙特卡洛，只不过稍微优化了一些，它的核心算法在下面:

$$\begin{array}{cc}\text{ 符号 }& \text{ 含义 }\\ n& \text{ 粒子个数 }\\ c_1& \text{ 粒子的个体学习因子, 也称为个体加速因子 }\\ c_2& \text{ 粒子的社会学习因子，也称为社会加速因子 }\\ w& \text{ 速度的惯性权重 }\\ v_i^d& \text{ 第d次迭代时，第i个粒子的速度 }\\ x_i^d& \text{ 第d次迭代时，第i个粒子所在的位置 }\\ f(x)& \text{ 在位置x时的适应度值(一般取目标函数值) }\\ pbes{t}_i^d& \text{ 到第d次迭代为止，第i个粒子经过的最好的位置 }\\ gbes{t}^d& \text{ 到第d次迭代为止，所有粒子经过的最好的位置 }\end{array}$$

这只鸟第 $d$ 步的速度 $=$ 上一步自身的速度惯性 $＋$ 自我认知部分 $+$ 社会认知部分（每—步运动的时间 $t$ —般取 1）

$${\mathbf{v}}_i^d=\mathbf{w}{\mathbf{v}}_i^{d-1}+{\mathbf{c}}_1{\mathbf{r}}_1\left(\mathbf{p}\mathbf{b}\mathbf{e}\mathbf{s}{\mathbf{t}}_i^d-{\mathbf{x}}_i^d\right)+{\mathbf{c}}_2{\mathbf{r}}_2\left(\mathbf{g}\mathbf{b}\mathbf{e}\mathbf{s}{\mathbf{t}}^d-{\mathbf{x}}_i^d\right)$$

这只鸟第 $d$ 步所在的位置 $=$ 第 $d-1$ 步所在的位置 $＋$ 第 $d-1$ 步的速度 $\ast$ 运动的时间（三个部分的和）

$${\mathbf{x}}_i^{\mathbf{d}+1}={\mathbf{x}}_i^{\mathbf{d}}+{\mathbf{v}}_{\mathbf{i}}^d$$

【流程图】

粒子群算法的作用：主要是求解函数的最值，但也可以用于求解方程组、多元函数的拟合、拟合微分方程。不过后面三个属于粒子群的进阶应用了，我打算放在文章最后

改进：

一、惯性权重

① 线性递减惯性权重

解析：惯性权重随着迭代次数减少，有利于先搜索全局再搜索局部最优解

$w_1^d=w_{\text{start}}-\left(w_{\text{start}}-w_{\text{end}}\right)\times \frac{d}{K}$

其中 $d$ 是当前迭代的次数，$K$ 是迭代总次数，$w_{\text{start}}$ 一般取 0.9，$w_{\text{end}}$一般取 0.4

② 非线性递减惯性权重

解析：原理与上一个相似，只不过惯性权重是非线性递减的

$w_2^d=w_{\text{start}}-\left(w_{\text{start}}-w_{\text{end}}\right)\times {\left(\frac{d}{K}\right)}^2$

$w_3^d=w_{\text{start}}-\left(w_{\text{start}}-w_{\text{end}}\right)\times \left[\frac{2d}{K}-{\left(\frac{d}{K}\right)}^2\right]$

【线性递减与非线性递减图示】

③ 自适应惯性权重（以求最小值为例）

解析：适应度越小，说明距离最优解越近，此时更需要局部搜索，即需要一个较小的惯性权重；反之，适应度越大，说明距离最优解越远，此时更需要全局搜索，即需要一个较大的惯性权重

$$w_i^d=\{\begin{array}{c}{w}_{\min }+\left(w_{\max }-w_{\min }\right)\frac{f\left(x_i^d\right)-f_{\min }^d}{f_{\text{average}}^d-f_{\min }^d},f\left(x_i^d\right)\le f_{\text{average}}^d\\ w_{\max },f\left(x_i^d\right)>f_{\text{average}}^d\end{array}$$

其中：

$w_{\text{start}}$ 一般也取 0.9，$w_{\text{end}}$一般也取 0.4

$f_{\text{average}}^d={\sum }_{i=1}f\left(x_i^d\right)/n$ ，即第 $d$ 次迭代时所有粒子的平均适应度

$f_{\min }^d=\min \left\{f\left(x_1^d\right),f\left(x_2^d\right),\cdots ,f\left(x_n^d\right)\right\}$，即第 $d$ 次迭代时所有粒子的最小适应度

④ 随机惯性权重

解析：可以避免在迭代前期局部搜索能力的不足，也可以避免在迭代后期全局搜索能力的不足

$\omega ={\mu }_{\min }+\left({\mu }_{\max }-{\mu }_{\min }\right)\times \mathrm{rand}()+\sigma \times \mathrm{randn}()$

其中：

${\mu }_{\min }$是随机惯性权重的最小值，${\mu }_{\max }$是随机惯性权重的最大值

$rand()$为 [0,1] 均匀分布随机数，$randn()$为正态分布的随机数

$\sigma$ (方差) 用来度量随机变量权重 $\omega$ 与其数学期望之间的偏离程度，使实验误差服从正态分布

二、学习因子

① 压缩因子

解析：可确保粒子群算法的收敛性，并可取消对速度的边界限制

$c1=c2=2.05,C=c1+c2=4.1,$ 收缩因子 $\mathbf{\Phi }=\frac{2}{|\left(2-\mathbf{C}-\sqrt{{\mathbf{C}}^2-4\mathbf{C}}\right)|}$

速度更新的公式改为：

${\mathbf{v}}_i^d=\mathbf{\Phi }\times \left[\mathbf{w}{\mathbf{v}}_i^{d-1}+{\mathbf{c}}_1{\mathbf{r}}_1\left(\mathbf{p}\mathbf{b}\mathbf{e}\mathbf{s}{\mathbf{t}}_i^d-{\mathbf{x}}_i^d\right)+{\mathbf{c}}_2{\mathbf{r}}_2\left(\mathbf{g}\mathbf{b}\mathbf{e}\mathbf{s}{\mathbf{t}}^d-{\mathbf{x}}_i^d\right)\right]$

② 非对称学习因子

解析：在算法搜索初期采用较大的 $c_1$ 值和较小的 $c_2$ 值，使粒子尽量发散到搜索空间,增加全局搜索的能力。随着迭代次数的增加，使 $c_1$ 线性递减，$c_2$ 线性递增，从而加强了粒子向全局最优点的收敛能力，也就是局部搜索的能力

$$\begin{array}{rl}{c}_1^{ini}=2.5,c_1^{fin}& =0.5,c_2^{ini}=1,c_2^{fin}=2.25\\ c_1^d=c_1^{ini}+\left(c_1^{fin}-c_1^{ini}\right)& \times \frac{d}{K};c_2^d=c_2^{ini}+\left(c_2^{fin}-c_2^{ini}\right)\times \frac{d}{K}\end{array}$$

实践：

我选择使用自适应惯性权重和收缩因子做为改进的粒子群算法，并使用 $MATLAB$ 自带粒子群算法进行验算

（PS：代码不作讲解，毕竟这不是教学，只是笔记罢了）

【代码部分】

$MATLAB$ 自带粒子群算法函数代码：

code.m

close all
clear;clc
narvs = 2; % 变量个数
x_lb = [-3 4.1]; % x的下界(长度等于变量的个数，每个变量对应一个下界约束)
x_ub = [12.1 5.8]; % x的上界

tic
options = optimoptions('particleswarm','FunctionTolerance',1e-12,'MaxStallIterations',50,'MaxIterations',20000,'HybridFcn',@fmincon,'PlotFcn','pswplotbestf');
title('Best Function Value: -38.8503','Interpreter','none');
[x,fval] = particleswarm(@Obj_fun,narvs,x_lb,x_ub,options);
fval = -fval
toc

% 最优值：38.7328

Obj_fun.m

function y = Obj_fun(x)  
    y = -(21.5 + x(1)*sin(4*pi*x(1)) + x(2)*sin(20*pi*x(2)));
end

【结果图】

自适应惯性权重和收缩因子优化粒子群算法函数：

realcode.m

close all
clear;clc
n = 50; % 粒子数量
narvs = 2; % 变量个数
c1 = 2;  % 每个粒子的个体学习因子，也称为个体加速常数
c2 = 2;  % 每个粒子的社会学习因子，也称为社会加速常数
C = c1 + c2; % 因子之和
w_min = 0.4;  % 最小惯性权重
w_max = 0.9;  % 最小惯性权重
fai = 2/(2-C-(C^2-4*C)^0.5); % 收缩因子
K = 500;  % 迭代的次数
vmax = 1.2; % 粒子的最大速度
x_lb = [-3 4.1]; % x的下界(长度等于变量的个数，每个变量对应一个下界约束)
x_ub = [12.1 5.8]; % x的上界

for i = 1: narvs
    x(:,i) = x_lb(i) + (x_ub(i)-x_lb(i))*rand(n,1);    % 随机初始化粒子所在的位置在定义域内
end
v = -vmax + 2*vmax .* rand(n,narvs);  % 随机初始化粒子的速度（这里我们设置为[-vmax,vmax]）

%% 计算适应度
fit = zeros(n,1);  % 初始化这n个粒子的适应度全为0
for i = 1:n  % 循环整个粒子群，计算每一个粒子的适应度
    fit(i) = Obj_fun_2(x(i,:));   % 调用Obj_fun_2函数来计算适应度
end 
pbest = x;   % 初始化这n个粒子迄今为止找到的最佳位置（是一个n*narvs的向量）
ind = find(fit == max(fit), 1);  % 找到适应度最大的那个粒子的下标
gbest = x(ind,:);  % 定义所有粒子迄今为止找到的最佳位置（是一个1*narvs的向量）

%% 函数图像
ezsurf(@(x,y)(21.5 + x.*sin(4.*pi.*x) + y.*sin(20.*pi.*y)),[-10 10 -10 10]); % 以x和y代替x1和x2
hold on;

%% 在图上标上这n个粒子的位置用于演示
h = scatter3(x(:,1),x(:,2),fit,'*r');  % scatter3是绘制三维散点图的函数（这里返回h是为了得到图形的句柄，未来我们对其位置进行更新）

%% 迭代K次来更新速度与位置
fitnessbest = ones(K,1);  % 初始化每次迭代得到的最佳的适应度
for d = 1:K  % 开始迭代，一共迭代K次
    for i = 1:n   % 依次更新第i个粒子的速度与位置 % 自适应惯性权重
        f_i = fit(i);  % 取出第i个粒子的适应度
        f_avg = sum(fit)/n;  % 计算此时适应度的平均值
        f_max = max(fit); % 计算此时适应度的最大值
        if f_i >= f_avg
             if f_avg ~= f_max  % 如果分母为0，我们就干脆令w=w_min
                w = w_min + (w_max - w_min)*(f_max - f_i)/(f_max - f_avg);
            else
                w = w_min;
             end
        else
            w = w_max;
        end
        v(i,:) = fai * (w*v(i,:) + c1*rand(1)*(pbest(i,:) - x(i,:)) + c2*rand(1)*(gbest - x(i,:)));  % 更新第i个粒子的速度
        x(i,:) = x(i,:) + v(i,:); % 更新第i个粒子的位置
        % 如果粒子的位置超出了定义域，就对其进行调整
        for j = 1: narvs
            if x(i,j) < x_lb(j)
                x(i,j) = x_lb(j);
            elseif x(i,j) > x_ub(j)
                x(i,j) = x_ub(j);
            end
        end
        fit(i) = Obj_fun_2(x(i,:));  % 重新计算第i个粒子的适应度
        if fit(i) > Obj_fun_2(pbest(i,:))   % 如果第i个粒子的适应度大于这个粒子迄今为止找到的最佳位置对应的适应度
           pbest(i,:) = x(i,:);   % 那就更新第i个粒子迄今为止找到的最佳位置
        end
        if  fit(i) > Obj_fun_2(gbest)  % 如果第i个粒子的适应度大于所有的粒子迄今为止找到的最佳位置对应的适应度
            gbest = pbest(i,:);   % 那就更新所有粒子迄今为止找到的最佳位置
        end
    end
    fitnessbest(d) = Obj_fun_2(gbest);  % 更新第d次迭代得到的最佳的适应度
    pause(0.1)  % 暂停0.1s
    h.XData = x(:,1);  % 更新散点图句柄的x轴的数据（此时粒子的位置在图上发生了变化）
    h.YData = x(:,2);   % 更新散点图句柄的y轴的数据（此时粒子的位置在图上发生了变化）
    h.ZData = fit;  % 更新散点图句柄的z轴的数据（此时粒子的位置在图上发生了变化）
end

figure(2) 
plot(fitnessbest)  % 绘制出每次迭代最佳适应度的变化图
xlabel('迭代次数');
disp('最佳的位置是：'); disp(gbest)
disp('此时最优值是：'); disp(Obj_fun_2(gbest))

% 最优值：38.6247

Obj_fun_2.m

function y = Obj_fun_2(x)
    y = 21.5 + x(1)*sin(4*pi*x(1)) + x(2)*sin(20*pi*x(2));
end

【动态图】

【结果图】

进阶：

上面讲了粒子群求解函数最值问题，现在针对求解方程组、多元函数的拟合、拟合微分方程进行实践分析

一、求解方程组

这里可以使用三个方法进行求解方程组，vpasolve 函数、fsolve 函数和粒子群，前面两个方法需要指定初始值和搜索范围，而粒子群只用指定搜索范围，相比起来会更加简单。注意，粒子群和两个函数使用的方程组函数不一样（具体看下面代码）

【代码部分】

main.m

%% vpasolve 函数

clear; clc
syms x1 x2 x3
eqn =  [abs(x1+x2)-abs(x3)  == 0, x1*x2*x3+18 == 0, x1^2*x2+3*x3 == 0]
% [x1, x2, x3] = vpasolve(eqn, [x1, x2, x3], [0 0 0])  % 这个初始值不行
[x1, x2, x3] = vpasolve(eqn, [x1, x2, x3], [1 1 1])  % 换一个初始值试试

% x = 6.80305735323491        -0.414133803887221          6.38892354934769
% SSE：0


%% fsolve 函数

x0 = [0,0,0];  % 初始值
x = fsolve(@my_fun,x0)
% 换一个初始值试试
x0 = [1,1,1];  % 初始值
format long g  % 显示更多的小数位数
x = fsolve(@my_fun,x0)

% x = 6.80305735323491        -0.414133803887221          6.38892354934769
% SSE：7.105427357601e-15


%% 粒子群算法(尝试多次，找到最优解)

clear; clc
narvs = 3;
% 使用粒子群算法，不需要指定初始值，只需要给定一个搜索的范围
lb = -10*ones(1,3);  ub = 10*ones(1,3);  
options = optimoptions('particleswarm','FunctionTolerance',1e-12,'MaxStallIterations',100,'MaxIterations',20000,'SwarmSize',100);
[x, fval] = particleswarm(@Obj_fun,narvs,lb,ub,options) 

% x = 1.42058091501338         -4.34008181244093          2.91950089742755
% SSE：2.22044604925031e-15

my_fun.m（vpasolve 函数、fsolve 函数）

function F = my_fun(x)
    F(1) =  abs(x(1)+x(2))-abs(x(3));
    F(2) =  x(1) * x(2) * x(3) + 18;
    F(3) = x(1)^2*x(2)+3*x(3);
end

Obj_fun.m（粒子群）

function f = Obj_fun(x)
    f1= abs(x(1)+x(2))-abs(x(3)) ;
    f2 = x(1) * x(2) * x(3) + 18;
    f3= x(1)^2 * x(2) + 3*x(3);
    f = abs(f1) + abs(f2) + abs(f3);
end

$SSE$是残差平方和，但为了图个方便，虽然使用这个名称，但运算的公式是abs(f1) + abs(f2) + abs(f3)，也就是 Obj_fun.m 中的最后一步，表示算法整体的误差值，当然越小越好，从最后结果来看显然是 vpasolve 函数的运算效果最好，因为它的 $SSE$ 是 0，而粒子群的 $SSE$ 也很小，在可接受范围内。因此选择哪个还得看具体模型环境，从简便和准确性之间权衡

二、多元函数的拟合

之所以不用 $MATLAB$ 的 $cftool$拟合工具箱是因为它仅有 $xyz$ 三元的函数拟合，当遇到多维拟合时，就可以用到粒子群了，但前提是知道它的拟合函数

【代码部分】

code.m

clear; clc
global x y;  % 将x和y定义为全局变量（方便在子函数中直接调用，要先声明）
load data_x_y  % 数据集内里面有x和y两个变量
narvs = 2;
% 使用粒子群算法，不需要指定初始值，只需要给定一个搜索的范围
lb = [-10 -10];  ub = [10 10];  
[k, fval] = particleswarm(@Obj_fun,narvs,lb,ub) 

% k = -0.0915    0.3169
% SSE：297.0634


% 使用粒子群算法后再利用fmincon函数混合求解
options = optimoptions('particleswarm','HybridFcn',@fmincon);
[k,fval] = particleswarm(@Obj_fun,narvs,lb,ub,options)

% k = -0.0915    0.3169
% SSE：297.0634

Obj_fun.m

function f = Obj_fun(k)
    global x y;  % 在子函数中使用全局变量前也需要声明下
    y_hat = exp(-k(1)*x(:,1)) .* sin(k(2)*x(:,2)) + x(:,3).^2;
    f = sum((y - y_hat) .^ 2);
end

从上面结果看出粒子群混合求解的方法和原方法的效果差不多

三、拟合微分方程

【求解方法】


现在的未知量：${\beta }_1、{\beta }_2、a、b$

【代码部分】

main.m

clear;clc
load mydata.mat  % 导入数据（共三列，分别是S,I,R的数量）
global true_s true_i true_r  n  % 定义全局变量
n = size(mydata,1);  % 一共有多少行数据
true_s = mydata(:,1);
true_i = mydata(:,2);
true_r = mydata(:,3);
figure(1)
plot(1:n,true_i,'b-',1:n,true_r,'g-') % 单独画出I和R的数量
legend('I','R')

% 粒子群算法来求解
% beta1,beta2,a,b
% 给定参数的搜索范围（根据题目来给），后期可缩小范围
lb = [0 0 -1 -1]; 
ub = [1 1 1 1];  
% lb = [0 0 -0.2 -0.2]; 
% ub = [0.3 0.3 0.2 0.2];  
options = optimoptions('particleswarm','Display','iter','SwarmSize',200);
[h, fval] = particleswarm(@Obj_fun, 4, lb, ub, options)  % h就是要优化的参数，fval是目标函数值, 第二个输入参数4代表我们有4个变量需要搜索

beta1 = h(1);  % t<=15时，易感染者与已感染者接触且被传染的强度
beta2 = h(2);  % t>15时，易感染者与已感染者接触且被传染的强度
a = h(3); % 康复率gamma = a+bt
b = h(4);  % 康复率gamma = a+bt
[t,p]=ode45(@(t,x) SIR_fun(t,x,beta1,beta2,a,b), [1:n],[true_s(1) true_i(1) true_r(1)]);   
p = round(p);  % 对p进行四舍五入(人数为整数)

figure(3)  % 真实的人数和拟合的人数放到一起看看
plot(1:n,true_i,'b-',1:n,true_r,'g-') 
hold on
plot(1:n,p(:,2),'b--',1:n,p(:,3),'g--') 
legend('I','R','拟合的I','拟合的R')

% 预测未来的趋势
N = 27; % 计算N天
[t,p]=ode45(@(t,x) SIR_fun(t,x,beta1,beta2,a,b), [1:N],[true_s(1) true_i(1) true_r(1)]);   
p = round(p);  % 对p进行四舍五入(人数为整数)
figure(4)
plot(1:n,true_i,'b-',1:n,true_r,'g-') 
hold on
plot(1:N,p(:,2),'b--',1:N,p(:,3),'g--') 
legend('I','R','拟合的I','拟合的R')

% beta1 = 0.2227
% beta2 = 0.1326
% a = 0.0488
% b = -0.0012
% SSE：157

SIR_fun.m

function dx=SIR_fun(t,x,beta1,beta2,a,b)  
    % beta1 : t<=15时，易感染者与已感染者接触且被传染的强度
    % beta2   : t>15时，易感染者与已感染者接触且被传染的强度
    % a  : 康复率gamma = a+bt
    % b  : 康复率gamma = a+bt
    if t <=15
        beta = beta1;
    else
        beta = beta2;
    end
    gamma = a + b * t;  
    dx = zeros(3,1);  % x(1)表示S  x(2)表示I  x(3)表示R
    C = x(1)+x(2);  % 传染病系统中的有效人群,也就是课件中的N' = S+I
    dx(1) = - beta*x(1)*x(2)/C;  
    dx(2) = beta*x(1)*x(2)/C - gamma*x(2);
    dx(3) = gamma*x(2); 
end

Obj_fun.m

function f = Obj_fun(h)  
    % 注意，h就是我们要估计的目标函数里面的参数
    % h的长度就是我们待估参数的个数
    global true_s  true_i  true_r n   % 在子函数中使用全局变量前也需要声明下
    beta1 = h(1);  % t<=15时，易感染者与已感染者接触且被传染的强度
    beta2 = h(2);  % t>15时，易感染者与已感染者接触且被传染的强度
    a = h(3); % 康复率gamma = a+bt
    b = h(4);  % 康复率gamma = a+bt
    [t,p]=ode45(@(t,x) SIR_fun(t,x,beta1,beta2,a,b) , [1:n],[true_s(1) true_i(1) true_r(1)]);   
    p = round(p);  % 对p进行四舍五入(人数为整数)
    % p的第一列是易感染者S的数量，p的第二列是患者I的数量，p的第三列是康复者R的数量
    f = sum((true_s - p(:,1)).^2  + (true_i -  p(:,2)).^2  + (true_r - p(:,3)).^2);
end

【结果图】

总结：

粒子群作为一个智能算法，在解决模型的运算方面有很强大的功能（无论是简洁性还是效率），但是仅此而已，真正建模最重要的还是模型的建立，前面也提及了，在有特定约束或方程组的条件下可以用粒子群求解，但若连公式都没有谈何解题，因此粒子群只是作为辅佐而存在的，论模型的建立还得阅读大量文献和教材