本文已参与「新人创作礼」活动，一起开启掘金创作之路。

一、前言

数模比赛中，经常见到有同学“套用”启发式算法(数模中常称为智能优化算法）去求解一些数模问题，事实上，很大一部分问题是不需要用到启发式算法求解的，Matlab中内置的函数足够我们使用了。但是如果遇到的优化问题特别复杂的话，启发式算法就是我们求解问题的一大法宝。

今天我们就来学习第一个智能算法：粒子群算法，其全称为粒子群优化算法(Particle Swarm Optimization, PSO)。它是通过模拟鸟群觅食行为而发展起来的一种基于群体协作的搜索算法。本节我们主要侧重学习其思想，并将其用于求解函数的最值问题，下一节我们会使用粒子群算法求解几类较难处理的优化类问题。

1.1什么是启发式算法？

注：“智能算法"是指在工程实践中提出的一些比较"新颖"的算法或理论，因此智能算法的范围要比启发式算法更大一点，如果某种智能算法可以用来解决优化问题，那么这种算法也可能认为是启发式算法。

启发式算法百度百科上的定义：一个基于直观或经验构造的算法，在可接受的花费下给出待解决优化问题的一个可行解。

什么是可接受的花费？计算时间和空间能接受（求解一个问题要几万年 or 一万台电脑）

2）什么是优化问题？工程设计中优化问题 (optimization problem）指在一定约束条件下,求解一个目标函数的最大值(或最小值)问题。注：实际上和我们之前学过的规划问题的定义一样，称呼不同而已。

3）什么是可行解？得到的结果能用于工程实践（不一定非要是最优解）

常见的启发式算法：粒子群、模拟退火、遗传算法、蚁群算法、禁忌搜索算法等等 (启发式算法解决的问题大同小异，只要前三个算法学会了在数学建模中就足够了）

1.2盲目搜索和启发式搜索

从最简单的优化问题说起：

思考：怎么找到这个一元函数的最大值？（只有一个上下界约束，即函数的定义域）

1.2.1 盲目搜索

假设图中的a=1,b=10,我们要找出连续函数y=f(x)在[1，10]的最大值。

枚举法 取 $x_1$ =1, $x_2$ =1.0001, $x_3$ =1.0002，⋯， $x_{9001}$ =10 分别计算 $f(x_1)$ , $f(x_2)$ ，⋯ $f(x_{9001})$ ，看其中哪个函数值最大。（如果不考虑计算时间，我们可以划分的更小，这样可以增加计算精度）
蒙特卡罗模拟 随机在[a,b]区间上取N个样本点，即 $x_1, x_2, x_3, ..., x_N$ 都是在[a,b]取的随机数分别计算 $f(x_1)$ , $f(x_2)$ ，⋯ $f(x_N)$ ，看其中哪个函数值最大。（如果不考虑计算时间，我们可以取更多的样本点，这样可以增加计算精度）

有什么问题？

如果现在要求最值的函数是一个多变量的函数，例如是一个10个变量的函数，那么就完蛋了！（要考虑的情况随着变量数呈指数增长，计算时间肯定不够）

1.2.2 启发式搜索

爬山法

在解空间中随机生成一个初始解；
根据初始解的位置，我们在下一步有两种走法：向左邻域走一步或者向右邻域走一步（走的步长越小越好，但会加大计算量）；
比较不同走法下目标函数的大小，并决定下一步往哪个方向走。上图中向左走目标函数较大，因此我们应该往左走一小步；
不断重复这个步骤，直到找到一个极大值点(或定义域边缘处)，此时我们就结束搜索。

爬山法的缺陷：特别容易找到局部最优解

1.2.3 盲目搜索还是启发式搜索？

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

例如：蒙特卡罗模拟用来求解优化问题就是盲目搜索，还有大家熟悉的枚举法也是盲目搜索。

1.2.4 关于“启发式”，可以简单总结为下面两点：

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

1.3 为什么要学习启发式算法？

从“比赛”的角度出发：

TSP(旅行商问题)这类组合优化问题
非线性整数规划问题（例如书店买书问题）之前学过的蒙特卡罗模拟实际上是随机找解来尝试，如果解集中存在的元素特别多，那么效果就不理想。

从“学习”的角度出发：

训练最优化的思维
提高自身的编程水平
各个专业都有广泛应用

1.4 参考 Matlab2016数值计算与智能算法

二、粒子群算法

1995年，美国学者Kennedy和Eberhart共同提出了粒子群算法，其基本思想源于对乌类群体行为进行建模与仿真的研究结果的启发。

它的核心思想是：利用群体中的个体对信息的共享使整个群体的运动在问题求解空间中产生从无序到有序的演化过程，从而获得问题的可行解。

最初提出的论文：Kennedy J，Eberhart R. Particle swarm optimization[c]//Proceedings of ICNN'95 - International Conference on Neural Networks. IEEE, 1995.

2.1 粒子群算法的直观解释

设想这样一个场景：一群鸟在搜索食物

假设：

所有的乌都不知道食物在哪

它们知道自己的当前位置距离食物有多远

它们知道离食物最近的鸟的位置

那么想一下这时候会发生什么？

首先，离食物最近的鸟会对其他的鸟说：兄弟们，你们快往我这个方向来，我这离食物最近与此同时，每只鸟在搜索食物的过程中，它们的位置也在不停变化，因此每只鸟也知道自己离食物最近的位置，这也是它们的一个参考。最后，鸟在飞行中还需要考虑一个惯性。

图形的进一步解释

c1 : 个体学习因子，也称为个体加速因子
c2 :社会学习因子，也称为社会加速因子
r1,r2:[0,1]上随机数
w :称为惯性权重，也称为惯性系数

iShot_2022-10-28_13.02.59.png

这只鸟第 $d$ 步所在的位置 = 第 $d‐1$ 步所在的位置 $+$ 第 $d‐1$ 步的速度 $*$ 运动的时间 $t$ (每一步运动的时间t一般取1)

x(d) = x(d‐1) + v(d‐1) * t

这只鸟第 $d$ 步的速度 $=$ 上一步自身的速度惯性 $+$ 自我认知部分 $+$ 社会认知部分

v(d) = w*v(d‐1) + c1*r1*(pbest(d)‐x(d)) + c2*r2*(gbest(d)‐x(d))

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

粒子：优化问题的候选解
位置：候选解所在的位置
速度：候选解移动的速度
适应度：评价粒子优劣的值，一般设置为目标函数值
个体最佳位置：单个粒子迄今为止找到的最佳位置
群体最佳位置：所有粒子迄今为止找到的最佳位置

2.3 粒子群算法流程图

graph TD
Start --> 初始化参数
初始化参数 --> 初始化粒子
初始化粒子 --> 计算适应度
计算适应度 --> 更新粒子速度和位置
更新粒子速度和位置 --> 重新计算粒子的适应度
重新计算粒子的适应度 --> 找到全局最优的粒子
找到全局最优的粒子 --> 是否到达迭代次数
是否到达迭代次数 --> |No| 更新粒子速度和位置
是否到达迭代次数 --> |Yes| 结束

2.4 符号说明

符号	说明
$n$	粒子个数
$c_1$	粒子的个体学习因子，也称为个体加速因子
$c_2$	粒子的社会学习因子，也称为社会加速因子
$w$	速度的惯性权重
$v_i^d$	第d次迭代时，第i个粒子的速度
$x_i^d$	第d次迭代时，第i个粒子所在的位置
$f(x)$	在位置x时的适应度值(一般取目标函数值)
$pbest_i^d$	到第d次迭代为止，第i个粒子经过的最好的位置
$gbest_i^d$	到第d次迭代为止，所有粒子经过的最好的位置

2.5 粒子群算法的核心公式

2.5.1 速度公式

这只鸟第 $d$ 步的速度 $=$ 上一步自身的速度惯性 $+$ 自我认知部分 $+$ 社会认知部分

v_i^d = wv_i^{d-1} + c_1r_1(pbest_i^d - x_i^d) + c_2r_2(gbest_i^d - x_i^d)

其中 $r_1,r_2是[0,1]$ 的随机数。

注意，这里我们没有在变量说明的表格中放入r1和r2这两个随机数，是因为他们表示的含义不太重要，我们只需要简单的交代一下就行

2.5.2 位置公式

这只鸟第 $d$ 步所在的位置 = 第 $d‐1$ 步所在的位置 $+$ 第 $d‐1$ 步的速度 $*$ 运动的时间 $t$ (每一步运动的时间t一般取1)

x_i^{d+1} = x_i^d + v_i^d

2.5.3 预设参数：学习因子 $c_1, c_2$

v_i^d = wv_i^{d-1} + c_1r_1(pbest_i^d - x_i^d) + c_2r_2(gbest_i^d - x_i^d)

在最初提出粒子群算法的论文中指出，个体学习因子和社会(或群体)学习因子取2比较合适。（注意：最初提出粒子群算法的这篇论文没有惯性权重）

最初提出的论文：Kennedy J , Eberhart R . Particle swarm optimization[C]// Proceedings of ICNN'95 ‐ International Conference on Neural Networks. IEEE, 1995.

2.5.4 预设参数：惯性权重 $w$

论文中得到的结论：惯性权重取0.9‐1.2是比较合适的，一般取0.9就行

引入惯性权重的论文：SHI,Y. A Modified Particle Swarm Optimizer[C]// Proc. of IEEE ICEC conference, Anchorage. 1998.

2.6 例1：求一元函数的最大值

$求函数 y = 11sin(x) + 7cos(5x) 在[-3,3]内的最大值$

%% 粒子群算法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);
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]）
%  注意：这种写法只支持2017及之后的Matlab，老版本的同学请自己使用repmat函数将向量扩充为矩阵后再运算。
% 即：v = -repmat(vmax, n, 1) + 2*repmat(vmax, n, 1) .* rand(n,narvs);  
% 注意：x的初始化也可以用一行写出来：  x = x_lb + (x_ub-x_lb).*rand(n,narvs) ，原理和v的计算一样
% 老版本同学可以用x = repmat(x_lb, n, 1) + repmat((x_ub-x_lb), n, 1).*rand(n,narvs) 

%% 计算适应度
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))

增加粒子数量会增加我们找到更好结果的可能性，但会增加运算的时间。
c1,c2,w这三个量有很多改进空间，未来我再来专门讲怎么去调整。
迭代的次数K越大越好吗？不一定哦，如果现在已经找到最优解了，再增加迭代次数是没有意义的。
这里出现了粒子的最大速度，是为了保证下一步的位置离当前位置别太远，一般取自变量可行域的10%至20%(不同文献看法不同)。
X的下界和上界是保证粒子不会飞出定义域。

注意：假设目标函数是二元函数，且x1和x2都位于-3至3之间，则:

narvs = 2
x_lb=[-3 -3]; x_ub=[3 3]; vmax=[1.2 1.2]

2.7 例2：求二元函数的最小值

注意: 用粒子群算法求解函数最小值时，粒子适应度的计算我们仍设置为目标函数值，但是此时我们希望找到适应度最小的解。因此希望大家不要用我们中文的内涵去理解这里的 “适应度” （中文的内涵就是越适应越好），为了避免混淆你可以就直接用目标函数值来代替适应度的说法。

%% 粒子群算法PSO: 求解函数y = x1^2+x2^2-x1*x2-10*x1-4*x2+60在[-15,15]内的最小值(动画演示）
clear; clc

%% 绘制函数的图形
x1 = -15:1:15;
x2 = -15:1:15;
[x1,x2] = meshgrid(x1,x2);
y = x1.^2 + x2.^2 - x1.*x2 - 10*x1 - 4*x2 + 60;
mesh(x1,x2,y)
xlabel('x1');  ylabel('x2');  zlabel('y');  % 加上坐标轴的标签
axis vis3d % 冻结屏幕高宽比，使得一个三维对象的旋转不会改变坐标轴的刻度显示
hold on  % 不关闭图形，继续在上面画图

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

%% 初始化粒子的位置和速度
x = zeros(n,narvs);
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]）
%  注意：这种写法只支持2017及之后的Matlab，老版本的同学请自己使用repmat函数将向量扩充为矩阵后再运算。
% 即：v = -repmat(vmax, n, 1) + 2*repmat(vmax, n, 1) .* rand(n,narvs);  
% 注意：x的初始化也可以用一行写出来：  x = x_lb + (x_ub-x_lb).*rand(n,narvs) ，原理和v的计算一样
% 老版本同学可以用x = repmat(x_lb, n, 1) + repmat((x_ub-x_lb), n, 1).*rand(n,narvs) 


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

%% 在图上标上这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个粒子的速度与位置
        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_fun2(x(i,:));  % 重新计算第i个粒子的适应度
        if fit(i) < Obj_fun2(pbest(i,:))   % 如果第i个粒子的适应度小于这个粒子迄今为止找到的最佳位置对应的适应度
           pbest(i,:) = x(i,:);   % 那就更新第i个粒子迄今为止找到的最佳位置
        end
        if  fit(i) < Obj_fun2(gbest)  % 如果第i个粒子的适应度小于所有的粒子迄今为止找到的最佳位置对应的适应度
            gbest = pbest(i,:);   % 那就更新所有粒子迄今为止找到的最佳位置
        end
    end
    fitnessbest(d) = Obj_fun2(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_fun2(gbest))

三、改进粒子群算法

3.1 改进：线性递减惯性权重

惯性权重 $w$ 体现的是粒子继承先前的速度的能力，Shi,Y最先将惯性权重 $w$ 引入到粒子群算法中，并分析指出一个较大的惯性权值有利于全局搜索，而一个较小的权值则更利于局部搜索。为了更好地平衡算法的全局搜索以及局部搜索能力，Shi,Y提出了线性递减惯性权重LDIW(Linear Decreasing Inertia Weight),公式如下：

w^d = w_{start} - (w_{start} - w_{end}) * d/K

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

与原来的相比，现在惯性权重和迭代次数有关

参考论文：Shi, Y. and Eberhart, R.C. (1999) Empirical Study of Particle Swarm Optimization. Proceedings of the 1999 Congress on Evolutionary Computation, Washington DC, 6‐9 July 1999, 1945‐1950.

%% 线性递减惯性权重的粒子群算法PSO: 求解函数y = x1^2+x2^2-x1*x2-10*x1-4*x2+60在[-15,15]内的最小值(动画演示）
clear; clc

%% 绘制函数的图形
x1 = -15:1:15;
x2 = -15:1:15;
[x1,x2] = meshgrid(x1,x2);
y = x1.^2 + x2.^2 - x1.*x2 - 10*x1 - 4*x2 + 60;
mesh(x1,x2,y)
xlabel('x1');  ylabel('x2');  zlabel('y');  % 加上坐标轴的标签
axis vis3d % 冻结屏幕高宽比，使得一个三维对象的旋转不会改变坐标轴的刻度显示
hold on  % 不关闭图形，继续在上面画图

%% 粒子群算法中的预设参数（参数的设置不是固定的，可以适当修改）
n = 30; % 粒子数量
narvs = 2; % 变量个数
c1 = 2;  % 每个粒子的个体学习因子，也称为个体加速常数
c2 = 2;  % 每个粒子的社会学习因子，也称为社会加速常数
w_start = 0.9;  % 初始惯性权重，通常取0.9
w_end = 0.4; % 粒子群最大迭代次数时的惯性权重，通常取0.4
K = 100;  % 迭代的次数
vmax = [6 6]; % 粒子的最大速度
x_lb = [-15 -15]; % x的下界
x_ub = [15 15]; % x的上界

%% 初始化粒子的位置和速度
x = zeros(n,narvs);
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_fun2(x(i,:));   % 调用Obj_fun2函数来计算适应度
end 
pbest = x;   % 初始化这n个粒子迄今为止找到的最佳位置（是一个n*narvs的向量）
ind = find(fit == min(fit), 1);  % 找到适应度最小的那个粒子的下标
gbest = x(ind,:);  % 定义所有粒子迄今为止找到的最佳位置（是一个1*narvs的向量）

%% 在图上标上这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个粒子的速度与位置       
        w = w_start - (w_start - w_end) * d / K;  
        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_fun2(x(i,:));  % 重新计算第i个粒子的适应度
        if fit(i) < Obj_fun2(pbest(i,:))   % 如果第i个粒子的适应度小于这个粒子迄今为止找到的最佳位置对应的适应度
           pbest(i,:) = x(i,:);   % 那就更新第i个粒子迄今为止找到的最佳位置
        end
        if  fit(i) < Obj_fun2(gbest)  % 如果第i个粒子的适应度小于所有的粒子迄今为止找到的最佳位置对应的适应度
            gbest = pbest(i,:);   % 那就更新所有粒子迄今为止找到的最佳位置
        end
    end
    fitnessbest(d) = Obj_fun2(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_fun2(gbest))


%% 三种递减惯性权重的对比
w_start = 0.9;  w_end = 0.4; 
K = 30;  d = 1:K; 
w1 = w_start-(w_start-w_end)*d/K;
w2 = w_start-(w_start-w_end)*(d/K).^2;
w3 = w_start-(w_start-w_end)*(2*d/K-(d/K).^2);
figure(3)
plot(d,w1,'r',d,w2,'b',d,w3,'g'); 
legend('w1','w2','w3')

3.1.1 拓展：非线性递减惯性权重

w^d = w_{start} - (w_{start} - w_{end}) * (d/K)^2

w^d = w_{start} - (w_{start} - w_{end}) * [2d/K - (d/K)^2]

3.1.2 三种递减惯性权重的图形

w^d = w_{start} - (w_{start} - w_{end}) * d/K

w^d = w_{start} - (w_{start} - w_{end}) * (d/K)^2

w^d = w_{start} - (w_{start} - w_{end}) * [2d/K - (d/K)^2]

3.2 改进：自适应惯性权重

3.2.1 假设现在求最小值问题，那么：

w_i^d=\left\{\begin{array}{cc} 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) \leq f_{\text {average }}^d \\ w_{\max } & , f\left(x_i^d\right)>f_{\text {average }}^d \end{array}\right.

其中：

$w_{min}, w_{max}$ 是我们预先给定的最小惯性系数和最大惯性系数，一般取0.4和0.9
$f_{average}^d = \sum_{i=1}^n f(x_i^d)/n$ 即第d次迭代时所有粒子的平均适应度
$f_{min}^d = min\{f(x_1^d),f(x_2^d),...,f(x_n)^d\}$ 即第d次迭代时所有粒子的最小适应度

与原来的相比，现在惯性权重和迭代次数以及每个粒子适应度有关

一个较大的惯性权值有利于全局搜索

而一个较小的权值则更利于局部搜索

假设现在一共五个粒子ABCDE，此时它们的适应度分别为1, 2, 3, 4, 5 取最大惯性权重为0.9，最小惯性权重为0.4 那么，这五个粒子的惯性权重应该为：0.4, 0.65, 0.9, 0.9, 0.9适应度越小，说明距离最优解越近，此时更需要局部搜索适应度越大，说明距离最优解越远，此时更需要全局搜索

3.2.2 假设现在求最大值问题，那么：

w_i^d=\left\{\begin{array}{cc} 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) \geq f_{\text {average }}^d \\ w_{\max } & , f\left(x_i^d\right)<f_{\text {average }}^d \end{array}\right.

与原来的相比，现在惯性权重和迭代次数以及每个粒子适应度有关

一个较大的惯性权值有利于全局搜索

而一个较小的权值则更利于局部搜索

其中：

$w_{min}, w_{max}$ 是我们预先给定的最小惯性系数和最大惯性系数，一般取0.4和0.9
$f_{average}^d = \sum_{i=1}^n f(x_i^d)/n$ 即第d次迭代时所有粒子的平均适应度
$f_{min}^d = min\{f(x_1^d),f(x_2^d),...,f(x_n)^d\}$ 即第d次迭代时所有粒子的最小适应度

假设现在一共五个粒子ABCDE，此时它们的适应度分别为1, 2, 3, 4, 5 取最大惯性权重为0.9，最小惯性权重为0.4 那么，这五个粒子的惯性权重应该为：0.9, 0.9, 0.9,0.65, 0.4适应度越小，说明距离最优解越远，此时更需要全局搜索适应度越大，说明距离最优解越近，此时更需要局部搜索

%% 自适应权重的粒子群算法PSO: 求解函数y = x1^2+x2^2-x1*x2-10*x1-4*x2+60在[-15,15]内的最小值(动画演示）
clear; clc

%% 绘制函数的图形
x1 = -15:1:15;
x2 = -15:1:15;
[x1,x2] = meshgrid(x1,x2);
y = x1.^2 + x2.^2 - x1.*x2 - 10*x1 - 4*x2 + 60;
mesh(x1,x2,y)
xlabel('x1');  ylabel('x2');  zlabel('y');  % 加上坐标轴的标签
axis vis3d % 冻结屏幕高宽比，使得一个三维对象的旋转不会改变坐标轴的刻度显示
hold on  % 不关闭图形，继续在上面画图

%% 粒子群算法中的预设参数（参数的设置不是固定的，可以适当修改）
n = 30; % 粒子数量
narvs = 2; % 变量个数
c1 = 2;  % 每个粒子的个体学习因子，也称为个体加速常数
c2 = 2;  % 每个粒子的社会学习因子，也称为社会加速常数
w_max = 0.9;  % 最大惯性权重，通常取0.9
w_min = 0.4; % 最小惯性权重，通常取0.4
K = 100;  % 迭代的次数
vmax = [6 6]; % 粒子的最大速度
x_lb = [-15 -15]; % x的下界
x_ub = [15 15]; % x的上界

%% 初始化粒子的位置和速度
x = zeros(n,narvs);
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_fun2(x(i,:));   % 调用Obj_fun2函数来计算适应度
end 
pbest = x;   % 初始化这n个粒子迄今为止找到的最佳位置（是一个n*narvs的向量）
ind = find(fit == min(fit), 1);  % 找到适应度最小的那个粒子的下标
gbest = x(ind,:);  % 定义所有粒子迄今为止找到的最佳位置（是一个1*narvs的向量）

%% 在图上标上这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_min = min(fit); % 计算此时适应度的最小值
        if f_i <= f_avg  
            if f_avg ~= f_min  % 如果分母为0，我们就干脆令w=w_min
                w = w_min + (w_max - w_min)*(f_i - f_min)/(f_avg - f_min);
            else
                w = w_min;
            end
        else
            w = w_max;
        end
        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_fun2(x(i,:));  % 重新计算第i个粒子的适应度
        if fit(i) < Obj_fun2(pbest(i,:))   % 如果第i个粒子的适应度小于这个粒子迄今为止找到的最佳位置对应的适应度
           pbest(i,:) = x(i,:);   % 那就更新第i个粒子迄今为止找到的最佳位置
        end
        if  fit(i) < Obj_fun2(gbest)  % 如果第i个粒子的适应度小于所有的粒子迄今为止找到的最佳位置对应的适应度
            gbest = pbest(i,:);   % 那就更新所有粒子迄今为止找到的最佳位置
        end
    end
    fitnessbest(d) = Obj_fun2(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_fun2(gbest))

3.3 改进：随机惯性权重

使用随机的惯性权重，可以避免在迭代前期局部搜索能力的不足；也可以避免在迭代后期全局搜索能力的不足。

%% 随机权重的粒子群算法PSO: 求解函数y = x1^2+x2^2-x1*x2-10*x1-4*x2+60在[-15,15]内的最小值(动画演示）
clear; clc

%% 绘制函数的图形
x1 = -15:1:15;
x2 = -15:1:15;
[x1,x2] = meshgrid(x1,x2);
y = x1.^2 + x2.^2 - x1.*x2 - 10*x1 - 4*x2 + 60;
mesh(x1,x2,y)
xlabel('x1');  ylabel('x2');  zlabel('y');  % 加上坐标轴的标签
axis vis3d % 冻结屏幕高宽比，使得一个三维对象的旋转不会改变坐标轴的刻度显示
hold on  % 不关闭图形，继续在上面画图

%% 粒子群算法中的预设参数（参数的设置不是固定的，可以适当修改）
n = 30; % 粒子数量
narvs = 2; % 变量个数
c1 = 2;  % 每个粒子的个体学习因子，也称为个体加速常数
c2 = 2;  % 每个粒子的社会学习因子，也称为社会加速常数
w_max = 0.9;  % 最大惯性权重，通常取0.9
w_min = 0.4; % 最小惯性权重，通常取0.4
sigma = 0.3; % 正态分布的随机扰动项的标准差（一般取0.2-0.5之间）
K = 100;  % 迭代的次数
vmax = [6 6]; % 粒子的最大速度
x_lb = [-15 -15]; % x的下界
x_ub = [15 15]; % x的上界

%% 初始化粒子的位置和速度
x = zeros(n,narvs);
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_fun2(x(i,:));   % 调用Obj_fun2函数来计算适应度
end 
pbest = x;   % 初始化这n个粒子迄今为止找到的最佳位置（是一个n*narvs的向量）
ind = find(fit == min(fit), 1);  % 找到适应度最小的那个粒子的下标
gbest = x(ind,:);  % 定义所有粒子迄今为止找到的最佳位置（是一个1*narvs的向量）

%% 在图上标上这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个粒子的速度与位置
        w = w_min + (w_max - w_min)*rand(1) + sigma*randn(1);
        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_fun2(x(i,:));  % 重新计算第i个粒子的适应度
        if fit(i) < Obj_fun2(pbest(i,:))   % 如果第i个粒子的适应度小于这个粒子迄今为止找到的最佳位置对应的适应度
           pbest(i,:) = x(i,:);   % 那就更新第i个粒子迄今为止找到的最佳位置
        end
        if  fit(i) < Obj_fun2(gbest)  % 如果第i个粒子的适应度小于所有的粒子迄今为止找到的最佳位置对应的适应度
            gbest = pbest(i,:);   % 那就更新所有粒子迄今为止找到的最佳位置
        end
    end
    fitnessbest(d) = Obj_fun2(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_fun2(gbest))

最开始提出随机惯性权重的论文：Zhang L , Yu H , Hu S . A New Approach to Improve Particle Swarm Optimization[J]. lecture notes in computer science, 2003, 2723:134‐139.

3.4 改进：压缩(收缩)因子法

个体学习因子c1和社会(群体)学习因子c2决定了粒子本身经验信息和其他粒子的经验信息对粒子运行轨迹的影响，其反映了粒子群之间的信息交流。设置c1较大的值，会使粒子过多地在自身的局部范围内搜索，而较大的c2的值，则又会促使粒子过早收敛到局部最优值。为了有效地控制粒子的飞行速度，使算法达到全局搜索与局部搜索两者间的有效平衡，Clerc构造了引入收缩因子的 PSO模型，采用了压缩因子，这种调整方法通过合适选取参数，可确保PSO算法的收敛性，并可取消对速度的边界限制。

%% 带有压缩因子(收缩因子)的粒子群算法PSO: 求解函数y = x1^2+x2^2-x1*x2-10*x1-4*x2+60在[-15,15]内的最小值(动画演示）
clear; clc

%% 绘制函数的图形
x1 = -15:1:15;
x2 = -15:1:15;
[x1,x2] = meshgrid(x1,x2);
y = x1.^2 + x2.^2 - x1.*x2 - 10*x1 - 4*x2 + 60;
mesh(x1,x2,y)
xlabel('x1');  ylabel('x2');  zlabel('y');  % 加上坐标轴的标签
axis vis3d % 冻结屏幕高宽比，使得一个三维对象的旋转不会改变坐标轴的刻度显示
hold on  % 不关闭图形，继续在上面画图

%% 粒子群算法中的预设参数（参数的设置不是固定的，可以适当修改）
n = 30; % 粒子数量
narvs = 2; % 变量个数
c1 = 2.05;  % 每个粒子的个体学习因子，也称为个体加速常数
c2 = 2.05;  % 每个粒子的社会学习因子，也称为社会加速常数
C = c1+c2; 
fai = 2/abs((2-C-sqrt(C^2-4*C))); % 收缩因子
w = 0.9;  % 惯性权重 
K = 100;  % 迭代的次数
vmax = [6 6]; % 粒子的最大速度
x_lb = [-15 -15]; % x的下界
x_ub = [15 15]; % x的上界

%% 初始化粒子的位置和速度
x = zeros(n,narvs);
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_fun2(x(i,:));   % 调用Obj_fun2函数来计算适应度
end 
pbest = x;   % 初始化这n个粒子迄今为止找到的最佳位置（是一个n*narvs的向量）
ind = find(fit == min(fit), 1);  % 找到适应度最小的那个粒子的下标
gbest = x(ind,:);  % 定义所有粒子迄今为止找到的最佳位置（是一个1*narvs的向量）

%% 在图上标上这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个粒子的速度与位置
        v(i,:) = fai * (w*v(i,:) + c1*rand(1)*(pbest(i,:) - x(i,:)) + c2*rand(1)*(gbest - x(i,:)));  % 更新第i个粒子的速度
        %  有时候会看到直接写成：0.657*v(i,:) +  1.496*rand(1)*(pbest(i,:) - x(i,:)) + 1.496*rand(1)*(gbest - x(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_fun2(x(i,:));  % 重新计算第i个粒子的适应度
        if fit(i) < Obj_fun2(pbest(i,:))   % 如果第i个粒子的适应度小于这个粒子迄今为止找到的最佳位置对应的适应度
           pbest(i,:) = x(i,:);   % 那就更新第i个粒子迄今为止找到的最佳位置
        end
        if  fit(i) < Obj_fun2(gbest)  % 如果第i个粒子的适应度小于所有的粒子迄今为止找到的最佳位置对应的适应度
            gbest = pbest(i,:);   % 那就更新所有粒子迄今为止找到的最佳位置
        end
    end
    fitnessbest(d) = Obj_fun2(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_fun2(gbest))

参考文献：M. Clerc. The swarm and queen: towards a deterministic and adaptive particle swarm op‐timization. Proc. Congress on Evolutionary Computation, Washington, DC,. Piscataway,NJ:IEEE Service Center (1999) 1951–1957

3.5 改进：非对称学习因子

%% 非对称学习因子的粒子群算法PSO: 求解函数y = x1^2+x2^2-x1*x2-10*x1-4*x2+60在[-15,15]内的最小值(动画演示）
clear; clc
% 毛开富, 包广清, 徐驰. 基于非对称学习因子调节的粒子群优化算法[J]. 计算机工程, 2010(19):188-190.（公式10中的减号应该改为加号）


%% 绘制函数的图形
x1 = -15:1:15;
x2 = -15:1:15;
[x1,x2] = meshgrid(x1,x2);
y = x1.^2 + x2.^2 - x1.*x2 - 10*x1 - 4*x2 + 60;
mesh(x1,x2,y)
xlabel('x1');  ylabel('x2');  zlabel('y');  % 加上坐标轴的标签
axis vis3d % 冻结屏幕高宽比，使得一个三维对象的旋转不会改变坐标轴的刻度显示
hold on  % 不关闭图形，继续在上面画图

%% 粒子群算法中的预设参数（参数的设置不是固定的，可以适当修改）
n = 30; % 粒子数量
narvs = 2; % 变量个数
c1_ini = 2.5;  % 个体学习因子的初始值
c1_fin = 0.5;  % 个体学习因子的最终值
c2_ini = 1;  % 社会学习因子的初始值
c2_fin = 2.25;  % 社会学习因子的最终值
w = 0.9;  % 惯性权重
K = 100;  % 迭代的次数
vmax = [6 6]; % 粒子的最大速度
x_lb = [-15 -15]; % x的下界
x_ub = [15 15]; % x的上界

%% 初始化粒子的位置和速度
x = zeros(n,narvs);
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_fun2(x(i,:));   % 调用Obj_fun2函数来计算适应度
end 
pbest = x;   % 初始化这n个粒子迄今为止找到的最佳位置（是一个n*narvs的向量）
ind = find(fit == min(fit), 1);  % 找到适应度最小的那个粒子的下标
gbest = x(ind,:);  % 定义所有粒子迄今为止找到的最佳位置（是一个1*narvs的向量）

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

%% 迭代K次来更新速度与位置
fitnessbest = ones(K,1);  % 初始化每次迭代得到的最佳的适应度
for d = 1:K  % 开始迭代，一共迭代K次
    %    w = 0.9 - 0.5*d/K;  % 效果不是很好的话可以同时考虑线性递减权重
        c1 = c1_ini + (c1_fin - c1_ini)*d/K;
        c2 = c2_ini + (c2_fin - c2_ini)*d/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_fun2(x(i,:));  % 重新计算第i个粒子的适应度
        if fit(i) < Obj_fun2(pbest(i,:))   % 如果第i个粒子的适应度小于这个粒子迄今为止找到的最佳位置对应的适应度
           pbest(i,:) = x(i,:);   % 那就更新第i个粒子迄今为止找到的最佳位置
        end
        if  fit(i) < Obj_fun2(gbest)  % 如果第i个粒子的适应度小于所有的粒子迄今为止找到的最佳位置对应的适应度
            gbest = pbest(i,:);   % 那就更新所有粒子迄今为止找到的最佳位置
        end
    end
    fitnessbest(d) = Obj_fun2(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_fun2(gbest))

参考文献：毛开富, 包广清, 徐驰. 基于非对称学习因子调节的粒子群优化算法[J]. 计算机工程, 2010(19):188‐190.

3.6 改进：自动退出迭代循环

当粒子已经找到最佳位置后，再增加迭代次数只会浪费计算时间，那么我们能否设计一个策略，能够自动退出迭代呢?

初始化最大迭代次数、计数器以及最大计数值（例如分别取100, 0, 20）
定义“函数变化量容忍度”，一般取非常小的正数，例如10^(‐6)；
在迭代的过程中，每次计算出来最佳适应度后，都计算该适应度和上一次迭代时最佳适应度的变化量(取绝对值)；
判断这个变化量和“函数变化量容忍度”的相对大小，如果前者小，则计数器加1；否则计数器清0；
不断重复这个过程，有以下两种可能:

① 此时还没有超过最大迭代次数，计数器的值超过了最大计数值，那么我们就跳出迭代循环，搜索结束。
② 此时已经达到了最大迭代次数，那么直接跳出循环，搜索结束。

%% 自适应权重且带有收缩因子的粒子群算法PSO: 求解四种不同测试函数的最小值(动画演示）
clear; clc

%% 粒子群算法中的预设参数（参数的设置不是固定的，可以适当修改）
tic % 开始计时
n = 1000; % 粒子数量
narvs = 30; % 变量个数
c1 = 2.05;  % 每个粒子的个体学习因子，也称为个体加速常数
c2 = 2.05;  % 每个粒子的社会学习因子，也称为社会加速常数
C = c1+c2; 
fai = 2/abs((2-C-sqrt(C^2-4*C))); % 收缩因子
w_max = 0.9;  % 最大惯性权重，通常取0.9
w_min = 0.4; % 最小惯性权重，通常取0.4
K = 1000;  % 迭代的次数
% Sphere函数
% vmax = 30*ones(1,30); % 粒子的最大速度
% x_lb = -100*ones(1,30); % x的下界
% x_ub = 100*ones(1,30); % x的上界
% Rosenbrock函数
vmax = 10*ones(1,30); % 粒子的最大速度
x_lb = -30*ones(1,30); % x的下界
x_ub = 30*ones(1,30); % x的上界
% Rastrigin函数
% vmax = 1.5*ones(1,30); % 粒子的最大速度
% x_lb = -5.12*ones(1,30); % x的下界
% x_ub = 5.12*ones(1,30); % x的上界
% Griewank函数
% vmax = 150*ones(1,30); % 粒子的最大速度
% x_lb = -600*ones(1,30); % x的下界
% x_ub = 600*ones(1,30); % x的上界

%% 改进：自动退出迭代循环
Count = 0; % 计数器初始化为0
max_Count = 30;  % 最大计数值初始化为30
tolerance = 1e-6;  % 函数变化量容忍度，取10^(-6)

%% 初始化粒子的位置和速度
x = zeros(n,narvs);
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_fun3(x(i,:));   % 调用Obj_fun3函数来计算适应度
end 
pbest = x;   % 初始化这n个粒子迄今为止找到的最佳位置（是一个n*narvs的向量）
ind = find(fit == min(fit), 1);  % 找到适应度最小的那个粒子的下标
gbest = x(ind,:);  % 定义所有粒子迄今为止找到的最佳位置（是一个1*narvs的向量）


%% 迭代K次来更新速度与位置
% fitnessbest = ones(K,1);  % 初始化每次迭代得到的最佳的适应度
% 注意我把上面这行注释掉了，因为我们不知道最后一共要迭代多少次，可能后面会自动跳出循环
% 初始化的目的一般是节省运行的时间，在这里这个初始化的步骤我们可以跳过
for d = 1:K  % 开始迭代，一共迭代K次
    tem = gbest; % 将上一步找到的最佳位置保存为临时变量
    for i = 1:n   % 依次更新第i个粒子的速度与位置
        f_i = fit(i);  % 取出第i个粒子的适应度
        f_avg = sum(fit)/n;  % 计算此时适应度的平均值
        f_min = min(fit); % 计算此时适应度的最小值
        if f_i <= f_avg  
            w = w_min + (w_max - w_min)*(f_i - f_min)/(f_avg - f_min);
        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_fun3(x(i,:));  % 重新计算第i个粒子的适应度
        if fit(i) < Obj_fun3(pbest(i,:))   % 如果第i个粒子的适应度小于这个粒子迄今为止找到的最佳位置对应的适应度
           pbest(i,:) = x(i,:);   % 那就更新第i个粒子迄今为止找到的最佳位置
        end
        if  fit(i) < Obj_fun3(gbest)  % 如果第i个粒子的适应度小于所有的粒子迄今为止找到的最佳位置对应的适应度
            gbest = pbest(i,:);   % 那就更新所有粒子迄今为止找到的最佳位置
        end
    end
    fitnessbest(d) = Obj_fun3(gbest);  % 更新第d次迭代得到的最佳的适应度
    delta_fit = abs(Obj_fun3(gbest) - Obj_fun3(tem));   % 计算相邻两次迭代适应度的变化量
    if delta_fit < tolerance  % 判断这个变化量和“函数变化量容忍度”的相对大小，如果前者小，则计数器加1
        Count = Count + 1;
    else
        Count = 0;  % 否则计数器清0
    end   
    if Count > max_Count  % 如果计数器的值达到了最大计数值
        break;  % 跳出循环
    end
end

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

toc % 结束计时

四、优化问题的测试函数

当你提出了一种新的优化算法后，你需要和别人之前提出的算法来进行PK，看你的算法有没有提高，下表给出了四种常见的测试函数：

%% 自适应权重且带有收缩因子的粒子群算法PSO: 求解四种不同测试函数的最小值(动画演示）
clear; clc

%% 粒子群算法中的预设参数（参数的设置不是固定的，可以适当修改）
tic % 开始计时
n = 1000; % 粒子数量
narvs = 30; % 变量个数
c1 = 2.05;  % 每个粒子的个体学习因子，也称为个体加速常数
c2 = 2.05;  % 每个粒子的社会学习因子，也称为社会加速常数
C = c1+c2; 
fai = 2/abs((2-C-sqrt(C^2-4*C))); % 收缩因子
w_max = 0.9;  % 最大惯性权重，通常取0.9
w_min = 0.4; % 最小惯性权重，通常取0.4
K = 1000;  % 迭代的次数
% Sphere函数
vmax = 30*ones(1,30); % 粒子的最大速度
x_lb = -100*ones(1,30); % x的下界
x_ub = 100*ones(1,30); % x的上界
% Rosenbrock函数
% vmax = 10*ones(1,30); % 粒子的最大速度
% x_lb = -30*ones(1,30); % x的下界
% x_ub = 30*ones(1,30); % x的上界
% Rastrigin函数
% vmax = 1.5*ones(1,30); % 粒子的最大速度
% x_lb = -5.12*ones(1,30); % x的下界
% x_ub = 5.12*ones(1,30); % x的上界
% Griewank函数
% vmax = 150*ones(1,30); % 粒子的最大速度
% x_lb = -600*ones(1,30); % x的下界
% x_ub = 600*ones(1,30); % x的上界

%% 初始化粒子的位置和速度
x = zeros(n,narvs);
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_fun3(x(i,:));   % 调用Obj_fun3函数来计算适应度
end 
pbest = x;   % 初始化这n个粒子迄今为止找到的最佳位置（是一个n*narvs的向量）
ind = find(fit == min(fit), 1);  % 找到适应度最小的那个粒子的下标
gbest = x(ind,:);  % 定义所有粒子迄今为止找到的最佳位置（是一个1*narvs的向量）


%% 迭代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_min = min(fit); % 计算此时适应度的最小值
        if f_i <= f_avg  
            w = w_min + (w_max - w_min)*(f_i - f_min)/(f_avg - f_min);
        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_fun3(x(i,:));  % 重新计算第i个粒子的适应度
        if fit(i) < Obj_fun3(pbest(i,:))   % 如果第i个粒子的适应度小于这个粒子迄今为止找到的最佳位置对应的适应度
           pbest(i,:) = x(i,:);   % 那就更新第i个粒子迄今为止找到的最佳位置
        end
        if  fit(i) < Obj_fun3(gbest)  % 如果第i个粒子的适应度小于所有的粒子迄今为止找到的最佳位置对应的适应度
            gbest = pbest(i,:);   % 那就更新所有粒子迄今为止找到的最佳位置
        end
    end
    fitnessbest(d) = Obj_fun3(gbest);  % 更新第d次迭代得到的最佳的适应度
end

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

toc % 结束计时

维数：自变量x的个数,也是上面表达式中n的大小
取值范围：每个x对应的变化范围
理论极值：这个函数理论上的最小值
误差目标：只要我们求出来的最小值小于这个目标值就能被接受

下面我们就来用粒子群算法求解这四个测试函数

张玮, 王华奎. 粒子群算法稳定性的参数选择策略分析[J]. 系统仿真学报, 2009, 21(014):4339‐4344.

五、粒子群算法进阶

纪震等《粒子群算法及应用》科学出版社，2009，P13‐14

六、Matlab自带的粒子群函数

6.1 示例

%% Matlab自带的粒子群函数 particleswarm
% particleswarm函数是求最小值的
% 如果目标函数是求最大值则需要添加负号从而转换为求最小值。
clear;clc
%  Matlab中粒子群算法函数的参考文献
%       Mezura-Montes, E., and C. A. Coello Coello. "Constraint-handling in nature-inspired numerical optimization: Past, present and future." Swarm and Evolutionary Computation. 2011, pp. 173–194.
%       Pedersen, M. E. "Good Parameters for Particle Swarm Optimization." Luxembourg: Hvass Laboratories, 2010.
%       Iadevaia, S., Lu, Y., Morales, F. C., Mills, G. B. & Ram, P. T. Identification of optimal drug combinations targeting cellular networks: integrating phospho-proteomics and computational network analysis. Cancer Res. 70, 6704-6714 (2010).
%       Liu, Mingshou , D. Shin , and H. I. Kang . "Parameter estimation in dynamic biochemical systems based on adaptive Particle Swarm Optimization." Information, Communications and Signal Processing, 2009. ICICS 2009. 7th International Conference on IEEE Press, 2010.

%% 求解函数y = x1^2+x2^2-x1*x2-10*x1-4*x2+60在[-15,15]内的最小值（最小值为8）
narvs = 2; % 变量个数
x_lb = [-15 -15]; % x的下界(长度等于变量的个数，每个变量对应一个下界约束)
x_ub = [15 15]; % x的上界
[x,fval,exitflag,output] = particleswarm(@Obj_fun2, narvs, x_lb, x_ub)   

%% 直接调用particleswarm函数进行求解测试函数
narvs = 30; % 变量个数
% Sphere函数
% x_lb = -100*ones(1,30); % x的下界
% x_ub = 100*ones(1,30); % x的上界
% Rosenbrock函数
x_lb = -30*ones(1,30); % x的下界
x_ub = 30*ones(1,30); % x的上界
% Rastrigin函数
% x_lb = -5.12*ones(1,30); % x的下界
% x_ub = 5.12*ones(1,30); % x的上界
% Griewank函数
% x_lb = -600*ones(1,30); % x的下界
% x_ub = 600*ones(1,30); % x的上界
[x,fval,exitflag,output] = particleswarm(@Obj_fun3,narvs,x_lb,x_ub)  


%% 绘制最佳的函数值随迭代次数的变化图
options = optimoptions('particleswarm','PlotFcn','pswplotbestf')   
[x,fval] = particleswarm(@Obj_fun3,narvs,x_lb,x_ub,options)

%% 展示函数的迭代过程
options = optimoptions('particleswarm','Display','iter');
[x,fval] = particleswarm(@Obj_fun3,narvs,x_lb,x_ub,options)

%% 修改粒子数量，默认的是：min(100,10*nvars)
options = optimoptions('particleswarm','SwarmSize',50);
[x,fval] = particleswarm(@Obj_fun3,narvs,x_lb,x_ub,options)

%% 在粒子群算法结束后继续调用其他函数进行混合求解（hybrid  n.混合物合成物; adj.混合的; 杂种的;） 
options = optimoptions('particleswarm','HybridFcn',@fmincon);
[x,fval] = particleswarm(@Obj_fun3,narvs,x_lb,x_ub,options)

%% 惯性权重的变化范围，默认的是0.1-1.1
options = optimoptions('particleswarm','InertiaRange',[0.2 1.2]);
[x,fval] = particleswarm(@Obj_fun3,narvs,x_lb,x_ub,options)

%% 个体学习因子，默认的是1.49（压缩因子）
options = optimoptions('particleswarm','SelfAdjustmentWeight',2);
[x,fval] = particleswarm(@Obj_fun3,narvs,x_lb,x_ub,options)

%% 社会学习因子，默认的是1.49（压缩因子）
options = optimoptions('particleswarm','SocialAdjustmentWeight',2);
[x,fval] = particleswarm(@Obj_fun3,narvs,x_lb,x_ub,options)

%% 最大的迭代次数，默认的是200*nvars
options = optimoptions('particleswarm','MaxIterations',10000);
[x,fval] = particleswarm(@Obj_fun3,narvs,x_lb,x_ub,options)

%% 领域内粒子的比例 MinNeighborsFraction，默认是0.25 
options = optimoptions('particleswarm','MinNeighborsFraction',0.2);
[x,fval] = particleswarm(@Obj_fun3,narvs,x_lb,x_ub,options)

%% 函数容忍度FunctionTolerance, 默认1e-6, 用于控制自动退出迭代的参数
options = optimoptions('particleswarm','FunctionTolerance',1e-8);
[x,fval] = particleswarm(@Obj_fun3,narvs,x_lb,x_ub,options)

%% 最大停滞迭代数MaxStallIterations, 默认20, 用于控制自动退出迭代的参数
options = optimoptions('particleswarm','MaxStallIterations',50);
[x,fval] = particleswarm(@Obj_fun3,narvs,x_lb,x_ub,options)

%% 不考虑计算时间，同时修改三个控制迭代退出的参数
tic
options = optimoptions('particleswarm','FunctionTolerance',1e-12,'MaxStallIterations',100,'MaxIterations',100000);
[x,fval] = particleswarm(@Obj_fun3,narvs,x_lb,x_ub,options)
toc

%% 在粒子群结束后调用其他函数进行混合求解
tic
options = optimoptions('particleswarm','FunctionTolerance',1e-12,'MaxStallIterations',50,'MaxIterations',20000,'HybridFcn',@fmincon);
[x,fval] = particleswarm(@Obj_fun3,narvs,x_lb,x_ub,options)
toc

6.2 Matlab粒子群函数的细节讲解

Matlab自带的这个函数优化的特别好，运行速度快且找到的解也比较精准，但是内部的实现比较复杂，后面我们会简单介绍这个函数的实现思路。（代码如何跑得快：少写循环和判断语句，多基于矩阵运算来进行操作）

注意：

这个函数在R2014b版本后推出，之前的老版本Matlab中不存在。
这个函数是求最小值的，如果目标函数是求最大值则需要添加负号从而转换为求最小值。

Matlab中 particleswarm 函数采用的是自适应的邻域模式

该函数主要参考的两篇文章：（该信息可在Matlab官网中查询） Mezura‐Montes, E., and C. A. Coello Coello. "Constraint‐handling in nature‐inspired numerical optimization: Past, present and future." Swarm and Evolutionary Computation. 2011, pp. 173–194. Pedersen, M. E. "Good Parameters for Particle Swarm Optimization." Luxembourg: Hvass Laboratories, 2010.

6.2 预设参数的选取

粒子个数 SwarmSize 默认设置为： min{100,10*nvars}, nvars是变量个数

惯性权重 InertiaRange 默认设置的范围为：[0.1,1.1]，注意，在迭代过程中惯性权重会采取自适应措施，随着迭代过程不断调整。

个体学习因子 SelfAdjustmentWeight 默认设置为：1.49 （和压缩因子的系数几乎相同）

社会学习因子 SocialAdjustmentWeight 默认设置为：1.49 （和压缩因子的系数几乎相同）

邻域内粒子的比例 MinNeighborsFraction 默认设置为：0.25，由于采取的是邻域模式，因此定义了一个“邻域最少粒子数目”： minNeighborhoodSize = max{2,(粒子数目*邻域内粒子的比例)的整数部分}，在迭代开始后，每个粒子会有一个邻域，初始时邻域内的粒子个数(记为Q)就等于“邻域最少粒子数目”，后续邻域内的粒子个数Q会自适应调整。

6.3 变量初始化和适应度的计算

速度初始化: 和我们之前的类似，只不过最大速度就是上界和下界的差额 vmax = ub – lb； v = ‐vmax + 2vmax . rand(n,narvs);

位置初始化: 和我们之前的类似会将每个粒子的位置均匀分布在上下界约束内

计算每个粒子的适应度适应度仍设置为我们要优化的目标函数，由于该函数求解的是最小值问题，因此，最优解应为适应度最小即目标函数越小的解。

初始化个体最优位置和我们自己写的代码一样，因为还没有开始循环，因此这里的个体最优位置就是我们初始化得到的位置。

初始化所有粒子的最优位置因为每个粒子的适应度我们都已经求出来了，所以我们只需要找到适应度最低的那个粒子，并记其位置为所有粒子的最优位置。

6.4 更新粒子的速度和位置

在每次迭代中，我们要分别更新每一个粒子的信息。例如: 对于现在要更新的粒子i ,我们要进行以下几个操作：

随机生成粒子i的邻域，邻域内一共Q个粒子(包含粒子i本身)，并找到这Q个粒子中位置最佳的那个粒子，此时它的目标函数值最小，记其位置为lbest；
更新粒子i的速度，公式为：v = wv + c1u1*(pbest‐x) + c2u2(lbest‐x)； (注：这里省略了一些上下标，大家应该能理解)，这个速度公式和基本粒子群算法最大的不同在于：这里的群体信息交流部分使用的是邻域内的最优位置，而不是整个粒子群的最优位置；
更新粒子i的位置，公式为： x=x+ v，这个和基本粒子群算法一样；
修正位置和速度：如果粒子i的位置超过了约束，就将其位置修改到边界处；另外如果这个粒子的位置在边界处，我们还需要查看其速度是否超过了最大速度，如果超过的话将这个速度变为0。(注意，如果是多元函数的话可能只有某个分量超过了约束，我们的修改只需要针对这个分量即可)
计算粒子i的适应度，如果小于其历史最佳的适应度，就更新粒子i的历史最佳位置为现在的位置；另外还需要判断粒子i的适应度是否要小于所有粒子迄今为止找到的最小适应度，如果小的话需要更新所有的粒子的最佳位置为粒子i的位置。

6.5 自适应调整参数

假设在第d次迭代过程中，所有的粒子的信息都已经更新好了，那么在开始下一次的迭代之前，需要更新模型中的参数，这里就体现了自适应过程。规则如下：记此时所有粒子的最小适应度为a, 上一次迭代完成后所有粒子的最小适应度为b。比较a和b的相对大小，如果a<b，则记flag = 1; 否则记flag = 0. 如果：flag = 0，那么做下面的操作：更新c=c+1 ；这里的c表示“停滞次数计数器”，在开始迭代前就初始为0 更新 Q = min{Q+ minNeighborhoodSize, SwarmSize} ；Q: 邻域内的粒子个数; minNeighborhoodSize: 邻域最少粒子数目； SwarmSize: 粒子的总数如果：flag = 1，那么做下面的操作：更新Q = minNeighborhoodSize 更新c = max{c‐1, 0} 判断c的大小，如果c< 2，则更新w= 2*w；如果c> 5，则更新w= w/2；这里的w是惯性权重，如果计算的结果超过了惯性权重的上界或低于下界都需要将其进行调整到边界处。自适应体现在: 如果适应度开始停滞时，粒子群搜索会从邻域模式向全局模式转换，一旦适应度开始下降，则又恢复到邻域模式，以免陷入局部最优。当适应度的停滞次数足够大时，惯性系数开始逐渐变小，从而利于局部搜索。

附录

function y = Obj_fun1(x)
    y = 11*sin(x) + 7*cos(5*x);
    %    y = -(11*sin(x) + 7*cos(5*x));  % 如果调用fmincon函数，则需要添加负号改为求最小值
end

function y = Obj_fun2(x)  
% y = x1^2+x2^2-x1*x2-10*x1-4*x2+60
% x是一个向量哦！
    y = x(1)^2 + x(2)^2 - x(1)*x(2) - 10*x(1) - 4*x(2) + 60; % 千万别写成了x1
end

function y = Obj_fun3(x)
% 注意，x是一个向量
% Sphere函数
%     y = sum(x.^2);  % x.^2 表示x中每一个元素分别计算平方

% Rosenbrock函数
    tem1 = x(1:end-1);  
    tem2 = x(2:end);
    y = sum(100 * (tem2-tem1.^2).^2 + (tem1-1).^2);

% Rastrigin函数
%     y = sum(x.^2-10*cos(2*pi*x)+10);

% Griewank函数
%     y = 1/4000*sum(x.*x)-prod(cos(x./sqrt(1:30)))+1;
end

粒子群算法从入门到高阶【全面详尽】