【原理与代码讲解】粒子群优化算法 PSO（基础版）

0 前言

Gitee 代码地址： https://gitee.com/futurelqh/GA

粒子群优化算法(Particle Swarm Optimization，PSO)，其思想源于对鸟群觅食过程的模拟，把优化问题的解域类比为鸟类的飞行空间，将优化问题的候选解类比为 鸟群中的鸟，每只鸟抽象为飞行空间内一个无质量无体积的微粒，问题需搜寻的最优解 则等同于鸟群所要寻找的食物。

1. 应用研究的现状

• (1) 函数优化，可以求解复杂函数的优化问题，包括带 约束的优化问题，通常这些函数是非常复杂的，主要表现为非线性、维数高和不可微等 特点，而且不少优化函数存在多个局部极值。
• (2) 人工神经网络学习，粒子群算法可以 提高人工神经网络自组织和自学习能力，还能完成人工神经网络中的多种任务，主要包 括：网络连接权值的训练和网络拓扑结构的训练。
• (3) 电力系统中的应用，粒子群算法 在电力系统中的应用涉及到很多方面，包括：系统规划、系统运行和控制等。

2. 粒子群算法原理

2.1 粒子群算法中的基本概念

2.2 粒子群算法直观解释

进一步对整个鸟群定义符号

2.3 符号说明

2.4 核心公式

2.6 理论定义

2.7 算法执行流程

3. 案例：MATLAB 实现粒子群算法 => 优化函数

测试函数

3.1 题目

3.2 解题思路

上述算法存在的问题： 上述算法属于暴力搜素解题，当自变量增多时，时间复杂度会非常高。

3.3 引入粒子群算法求解此问题

上述算法按照预定的策略实行搜索，在搜索过程中获取的中间信息不用来改进策略，称为盲目搜索；反之，如果利用了中间信息来改进搜索策略则称为启发式搜索。

启发式，可以简单总结为下面两点：

1. 任何有助于找到问题的最优解，但不能保证找到最优解的方法均是启发式方法；
2. 有助于加速求解过程和找到较优解的方法是启发式方法

目标函数即：适应度

Obj_fun1.m

function y = Obj_fun1(x)
    y = 11*sin(x) + 7*cos(5*x);
end

粒子群算法实现

%% 粒子群算法PSO: 求解函数y = 11*sin(x) + 7*cos(5*x)在[-3,3]内的最大值
clear; clc

%% 绘制函数的图形
x = -3:0.01:3;
y = 11*sin(x) + 7*cos(5*x);
figure(1)
plot(x,y,'b-')
title('y = 11*sin(x) + 7*cos(5*x)')
hold on  % 不关闭图形，继续在上面画图

%% 粒子群算法中的预设参数（参数的设置不是固定的，可以适当修改）
n = 10; % 粒子数量
narvs = 1; % 变量个数
c1 = 2;  % 每个粒子的个体学习因子，也称为个体加速常数
c2 = 2;  % 每个粒子的社会学习因子，也称为社会加速常数
w = 0.9;  % 惯性权重
K = 50;  % 迭代的次数
vmax = 1.2; % 粒子的最大速度
x_lb = -3; % x的下界
x_ub = 3; % x的上界

%% 初始化粒子的位置和速度
x = zeros(n,narvs); % 初始化位置为全 0
for i = 1: narvs
    % 随机初始化粒子所在的位置在定义域内，方法为下线 + (上线 - 下线) * [0,1] 的随机数
    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_fun1(x(i,:));   % 调用Obj_fun1函数来计算适应度（这里写成x(i,:)主要是为了和以后遇到的多元函数互通）
end
pbest = x;   % 初始化这n个粒子迄今为止找到的最佳位置（是一个n*narvs的向量）
ind = find(fit == max(fit), 1);  % 找到适应度最大的那个粒子的下标
gbest = x(ind,:);  % 定义所有粒子迄今为止找到的最佳位置（是一个1*narvs的向量）

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

%% 迭代K次来更新速度与位置
fitnessbest = ones(K,1);  % 初始化每次迭代得到的最佳的适应度
for d = 1:K  % 开始迭代，一共迭代K次
    for i = 1:n   % 依次更新第i个粒子的速度与位置
        v(i,:) = w*v(i,:) + c1*rand(1)*(pbest(i,:) - x(i,:)) + c2*rand(1)*(gbest - x(i,:));  % 更新第i个粒子的速度
        % 如果粒子的速度超过了最大速度限制，就对其进行调整
        for j = 1: narvs
            if v(i,j) < -vmax(j)
                v(i,j) = -vmax(j);
            elseif v(i,j) > vmax(j)
                v(i,j) = vmax(j);
            end
        end
        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_fun1(x(i,:));  % 重新计算第i个粒子的适应度
        if fit(i) > Obj_fun1(pbest(i,:))   % 如果第i个粒子的适应度大于这个粒子迄今为止找到的最佳位置对应的适应度
            pbest(i,:) = x(i,:);   % 那就更新第i个粒子迄今为止找到的最佳位置
        end
        if  fit(i) > Obj_fun1(gbest)  % 如果第i个粒子的适应度大于所有的粒子迄今为止找到的最佳位置对应的适应度
            gbest = pbest(i,:);   % 那就更新所有粒子迄今为止找到的最佳位置
        end
    end
    fitnessbest(d) = Obj_fun1(gbest);  % 更新第d次迭代得到的最佳的适应度
    pause(0.1)  % 暂停0.1s
    h.XData = x;  % 更新散点图句柄的x轴的数据（此时粒子的位置在图上发生了变化）
    h.YData = fit; % 更新散点图句柄的y轴的数据（此时粒子的位置在图上发生了变化）
end

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

执行结果

Fitness 曲线

3.4 案例二： 求解具有多变量的 De Jong’s function 1 函数

Obj_fun1.m

function y = Obj_fun1(x)
% x 是一个矩阵
     y = sum(x.^2, 1);  % x.^2 表示 x 中每一个元素分别计算平方
end

PSO.m

%% 粒子群算法PSO: 求解函数 f1(x)=sum(x(i)^2) 在[-5.12,5.12]内的最小值
clear; clc


%% 粒子群算法中的预设参数（参数的设置不是固定的，可以适当修改）
n = 1000; % 粒子数量
narvs = 30; % 变量个数
c1 = 2;  % 每个粒子的个体学习因子，也称为个体加速常数
c2 = 2;  % 每个粒子的社会学习因子，也称为社会加速常数  
w = 0.9;  % 惯性权重
K = 500;  % 迭代的次数
vmax = ones(1,30); % 粒子的最大速度
x_lb = -5.12*ones(1,30); % x的下界
x_ub = 5.12*ones(1,30); % x的上界

%% 初始化粒子的位置和速度
x = zeros(n,narvs); % 初始化位置为全 0
x = x_lb + (x_ub - x_lb) .* rand(n, narvs); % 随机初始化粒子所在的位置在定义域内，方法为下线 + (上线 - 下线) * [0,1] 的随机数
v = -vmax + 2*vmax .* rand(n,narvs);  % 随机初始化粒子的速度（这里我们设置为[-vmax,vmax]）

%% 计算适应度
fit = zeros(n,1);  % 初始化这 n 个粒子的适应度全为0
fit = Obj_fun1(x);

pbest = x;   % 初始化这n个粒子迄今为止找到的最佳位置（是一个n*narvs的向量）
ind = find(fit == min(fit), 2);  % 找到适应度最小的那个粒子的下标
gbest = x(ind,:);  % 定义所有粒子迄今为止找到的最佳位置（是一个1*narvs的向量）


%% 迭代K次来更新速度与位置
fitnessbest = ones(K,1);  % 初始化每次迭代得到的最佳的适应度
for d = 1:K  % 开始迭代，一共迭代K次
    
    v = w * v + c1 * rand(1) * (pbest - x) + c2 * rand(1) * (gbest - x);
    v(v > 1) = 1;
    v(v < -1) = -1 ;

    x = x + v; % 计算新位置

    x(x < -5.12) = -5.12;
    x(x > 5.12) = 5.12;

    fit = Obj_fun1(x); %计算新位置的适应度

    judge = Obj_fun1(pbest) > fit; % 替换每个个体的最佳位置

    pbest(judge,:) = x(judge,:);
    
    % 更新最佳位置
    ind = find(Obj_fun1(pbest) == min(Obj_fun1(pbest)), 2); 

    gbest = pbest(ind,:);

    fitnessbest(d) = Obj_fun1(gbest);  % 更新第d次迭代得到的最佳的适应度
end

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