智能优化算法：粒子群算法

1. 粒子群算法简介

假设在一个$D$维的目标搜索空间中，有$N$个例子组成的一个群落，其中第$i$个粒子表示为一个$D$维度的向量：$${\mathbf{X}}_i=[x_{i1},x_{i1},\cdots ,x_{iD}],i=1,2,\cdots ,N\text{\hspace{1em}(1)}$$
第$i$个粒子的“飞行”速度也是一个$D$维向量，记为
$${\mathbf{V}}_i=[V_{i1},V_{i2,},\cdots ,V_{iD}],i=1,2,\cdots ,N\text{\hspace{1em}(2)}$$
第$i$个粒子迄今为止搜索到的最优位置称为个体极值，记作
$${\mathbf{p}}_{\text{best}}=[p_{i1},p_{i2},\cdots ,p_{iD}],i=1,2,\cdots ,N\text{\hspace{1em}(3)}$$
整个粒子群搜索到的最优位置为全局极值，记为
$${\mathbf{g}}_{\text{best}}=[g_1,g_2,\cdots ,g_D]$$

在找到这两个最优值时，粒子根据下式来更新自己的速度和位置：
$$\{\begin{array}{l}{v}_{ij}(t+1)=v_{ij}(t)+c_1{r}_1(t)[p_{ij}(t)-x_{ij}(t)]+c_2{r}_2(t)[g_j(t)-x_{ij}(t)]\\ x_{ij}(t+1)=x_{ij}(t)+v_{ij}(t+1)\end{array}$$式中：$i=1,2,\cdots ,N$，$j=1,2,\cdots ,D$，$c_1$和$c_2$为学习因子，也称加速常数；$r_1$和$r_2$为$[0,1]$范围内的均匀随机数，$v_{ij}$是第$i$个粒子第$j$维度的速度分量，$v_{ij}\in [-v_{\max },v_{\max }]$，v_{\max}是常数，由用户设定来限制粒子的速度。

速度更新方程右三个部分组成：

1. “惯性”或“动量”部分：分赢了粒子的运动习惯，代表粒子有维持自己先前速度的趋势；
2. “认知”部分：反映了粒子对自身历史经验的记忆或回忆，代表粒子有向自身历史最佳位置逼近的趋势；
3. “社会”部分：反映了粒子间协同合作与知识共享的群体历史经验，代表粒子有向群体或邻域历史最佳位置逼近的趋势。

2. 粒子群算法内容

2.1 粒子群算法流程

2.2 参数说明

3. 粒子群算法代码实现

3.1 Matlab

3.1.1 PSO求解函数极值

求$f(x)=x-10\sin (5x)+7\cos (4x)$的最大值，其中$x$的范围为$x\in [0,10]$。

% 粒子群优化算法求解函数最值

% 系统初始化
close all
clear
clc

% 参数设置
N = 100;            % 粒子数
D = 1;              % 维度
T = 100;            % 最大迭代次数
c1 = 1.5;           % 学习因子1
c2 = 1.5;           % 学习因子2
w = 0.8;            % 惯性权重

% 变量范围
Xmax = 0;           % 自变量范围
Xmin = 10;          % 自变量范围
Vmax = 2;           % 速度范围
Vmin = -2;          % 速度范围


% 初始化粒子位置与速度
x = rand(N,D) * (Xmax-Xmin) + Xmin;
v = rand(N,D) * (Vmax-Vmin) + Vmin;

% 初始化个体最优位置和最优值
p = x;
pbest = ones(N,1);
for i = 1:N
    pbest(i) = f(x(i,:));
end

% 初始化全局最优位置和最优值
g = ones(1,D);
gbest = -inf;
for i = 1:N
    if(pbest(i) > gbest)
        g = p(i,:);
        gbest = pbest(i);
    end
end

% 记录迭代过程
gb = ones(1,T);

% 开始迭代
for i = 1:T
    for j = 1:N
        % 更新个体最优位置和最优值
        if (f(x(j,:)) > pbest(j))
            p(j,:) = x(j,:);
            pbest(j) = f(x(j,:));
        end
        % 更新全局最优位置和最优值
        if pbest(j) > gbest
            g = p(j,:);
            gbest = pbest(j);
        end
        
        % 更新位置和速度
        v(j,:) = w * v(j,:) + c1 * rand * (p(j,:) - x(j,:)) + c2 * rand * (g - x(j,:));
        x(j,:) = x(j,:) + v(j,:);
        
        % 范围检测
        for ii = 1:D
            if (v(j,ii) > Vmax) || (v(j,ii) < Vmin)
                v(j,ii) = rand * (Vmax-Vmin) + Vmin;
            end
            if (x(j,ii) > Xmax) || (x(j,ii) < Xmin)
                x(j,ii) = rand * (Xmax-Xmin) + Xmin;
            end
        end
    end
    gb(i) = gbest;
end % 结束迭代

% 绘图
figure('name','迭代曲线','color','w')
plot(gb,'b','linewidth',1);
grid on
set(gca,'FontName','Times New Roman','FontSize',12);
xlabel('迭代次数','FontName','宋体');
ylabel('最优解','FontName','宋体');
xlim([1,T])

X = linspace(Xmin,Xmax,1000);
Y = f(X);
figure('name','函数图像','color','w');
plot(X,Y,'b','linewidth',1);
grid on
hold on
set(gca,'FontName','Times New Roman','FontSize',12);
xlabel('$x$','interpreter','Latex');
ylabel('$f(x)$','interpreter','Latex');
title('$f(x)=x - 10\sin(5x) + 7\cos(4x)$','interpreter','Latex');
ylim([-20,30])
[ymax,xindex] = max(Y);
line1 = plot(X(xindex),ymax,'go');
line2 = plot(g,f(g),'r*');
hold off
legend([line1,line2],{'Max','PSO Result'},'Location','southeast');

% 函数
function y = f(x)
y = x + 10*sin(5*x) + 7 * cos(4*x);
end