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

一、前言

数模竞赛中，常常见到有同学“套用”启发式算法(数模中常称为智能优化算法）去求解一些数模问题，事实上，很大一部分问题是不需求用到启发式算法求解的，Matlab中内置的函数满足咱们运用了。可是假如遇到的优化问题特别杂乱的话，启发式算法便是咱们求解问题的一大法宝。

今日咱们就来学习第一个智能算法：粒子群算法，其全称为粒子群优化算法(Particle Swarm Optimization, PSO)。它是通过模仿鸟群寻食行为而发展起来的一种根据集体协作的查找算法。本节咱们首要侧重学习其思想，并将其用于求解函数的最值问题，下一节咱们会运用粒子群算法求解几类较难处理的优化类问题。

1.1什么是启发式算法？

注：“智能算法”是指在工程实践中提出的一些比较”新颖”的算法或理论，因而智能算法的规模要比启发式算法更大一点，假如某种智能算法能够用来处理优化问题，那么这种算法也或许认为是启发式算法。

启发式算法百度百科上的界说：一个根据直观或经历构造的算法，在可承受的花费下给出待处理优化问题的一个可行解。

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

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

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

常见的启发式算法： 粒子群、模仿退火、遗传算法、蚁群算法、禁忌查找算法等等 (启发式算法处理的问题迥然不同，只需前三个算法学会了在数学建模中就满足了）

1.2盲目查找和启发式查找

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

思考：怎样找到这个一元函数的最大值？（只要一个上下界束缚，即函数的界说域）

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

1.2.1 盲目查找

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

  • 枚举法 取x1x_1=1,x2x_2=1.0001,x3x_3=1.0002，⋯，x9001x_{9001} =10 别离核算f(x1)f(x_1),f(x2)f(x_2)，⋯ f(x9001)f(x_{9001})，看其间哪个函数值最大。 （假如不考虑核算时刻，咱们能够区分的更小，这样能够增加核算精度）

  • 蒙特卡罗模仿 随机在[a,b]区间上取N个样本点，即x1,x2,x3,…,xNx_1, x_2, x_3, …, x_N都是在[a,b]取的随机数别离核算f(x1)f(x_1),f(x2)f(x_2)，⋯ f(xN)f(x_N)，看其间哪个函数值最大。 （假如不考虑核算时刻，咱们能够取更多的样本点，这样能够增加核算精度）

有什么问题？

假如现在要求最值的函数是一个多变量的函数，例如是一个10个变量的函数，那么就完蛋了！（要考虑的状况跟着变量数呈指数增长，核算时刻必定不行）

1.2.2 启发式查找

  • 爬山法
  1. 在解空间中随机生成一个初始解；
  2. 根据初始解的方位，咱们鄙人一步有两种走法：向左邻域走一步或者向右邻域走一步（走的步长越小越好，但会加大核算量）；
  3. 比较不同走法下方针函数的巨细，并决议下一步往哪个方向走。上图中向左走方针函数较大，因而咱们应该往左走一小步；
  4. 不断重复这个进程，直到找到一个极大值点(或界说域边际处)，此刻咱们就完毕查找。

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

爬山法的缺陷：特别简单找到部分最优解

1.2.3 盲目查找仍是启发式查找？

依照预订的战略实施查找，在查找进程中获取的中心信息不用来改善战略，称为盲目查找； 反之，假如运用了中心信息来改善查找战略则称为启发式查找。

例如：蒙特卡罗模仿用来求解优化问题便是盲目查找，还有咱们熟悉的枚举法也是盲目查找。

1.2.4 关于“启发式”，能够简略总结为下面两点：

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

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

  • 从“竞赛”的角度动身：
  1. TSP(旅行商问题)这类组合优化问题
  2. 非线性整数规划问题（例如书店买书问题）之前学过的蒙特卡罗模仿实际上是随机找解来尝试，假如解会集存在的元素特别多，那么效果就不理想。
  • 从“学习”的角度动身：
  1. 练习最优化的思想
  2. 进步自身的编程水平
  3. 各个专业都有广泛运用

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 粒子群算法的直观解说

想象这样一个场景：一群鸟在查找食物

假定：

  1. 一切的乌都不知道食物在哪
  2. 它们知道自己的当时方位距离食物有多远
  3. 它们知道离食物最近的鸟的方位

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

首先，离食物最近的鸟会对其他的鸟说：兄弟们，你们快往我这个方历来，我这离食物最近与此一起，每只鸟在查找食物的进程中，它们的方位也在不停改动，因而每只鸟也知道自己离食物最近的方位，这也是它们的一个参阅。最终，鸟在飞行中还需求考虑一个惯性。

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

图形的进一步解说

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

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

这只鸟第dd步地点的方位 = 第d‐1d‐1步地点的方位 ++ 第d‐1d‐1步的速度 ∗* 运动的时刻tt (每一步运动的时刻t一般取1)

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

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

v(d)=w∗v(d‐1)+c1∗r1∗(pbest(d)‐x(d))+c2∗r2∗(gbest(d)‐x(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 符号阐明

符号 阐明
nn 粒子个数
c1c_1 粒子的个别学习因子，也称为个别加快因子
c2c_2 粒子的社会学习因子，也称为社会加快因子
ww 速度的惯性权重
vidv_i^d 第d次迭代时，第i个粒子的速度
xidx_i^d 第d次迭代时，第i个粒子地点的方位
f(x)f(x) 在方位x时的习惯度值(一般取方针函数值)
pbestidpbest_i^d 到第d次迭代停止，第i个粒子通过的最好的方位
gbestidgbest_i^d 到第d次迭代停止，一切粒子通过的最好的方位

2.5 粒子群算法的中心公式

2.5.1 速度公式

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

vid=wvid−1+c1r1(pbestid−xid)+c2r2(gbestid−xid)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)

其间r1,r2是[0,1]r_1,r_2是[0,1]的随机数。

留意，这儿咱们没有在变量阐明的表格中放入r1和r2这两个随机数，是由于他 们表明的含义不太重要，咱们只需求简略的交代一下就行

2.5.2 方位公式

这只鸟第dd步地点的方位 = 第d‐1d‐1步地点的方位 ++ 第d‐1d‐1步的速度 ∗* 运动的时刻tt (每一步运动的时刻t一般取1)

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

2.5.3 预设参数：学习因子 c1,c2c_1, c_2

vid=wvid−1+c1r1(pbestid−xid)+c2r2(gbestid−xid)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 预设参数：惯性权重 ww

论文中得到的定论：惯性权重取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]内的最大值求函数 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))
  1. 增加粒子数量会增加咱们找到更好成果的或许性，但会增加运算的时刻。
  2. c1,c2,w这三个量有许多改善空间，未来我再来专门讲怎样去调整。
  3. 迭代的次数K越大越好吗？不必定哦，假如现在现已找到最优解了，再增加迭代次数是没有意义的。
  4. 这儿出现了粒子的最大速度，是为了保证下一步的方位离当时方位别太远，一般取自变量可行域的10%至20%(不同文献看法不同)。
  5. 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 改善：线性递减惯性权重

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

wd=wstart−(wstart−wend)∗d/Kw^d = w_{start} – (w_{start} – w_{end}) * d/K

其间dd是当时迭代的次数，KK是迭代总次数，wstart一般取0.9，wend一般取0.4w_{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 拓宽：非线性递减惯性权重

wd=wstart−(wstart−wend)∗(d/K)2w^d = w_{start} – (w_{start} – w_{end}) * (d/K)^2
wd=wstart−(wstart−wend)∗[2d/K−(d/K)2]w^d = w_{start} – (w_{start} – w_{end}) * [2d/K – (d/K)^2]

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

wd=wstart−(wstart−wend)∗d/Kw^d = w_{start} – (w_{start} – w_{end}) * d/K
wd=wstart−(wstart−wend)∗(d/K)2w^d = w_{start} – (w_{start} – w_{end}) * (d/K)^2
wd=wstart−(wstart−wend)∗[2d/K−(d/K)2]w^d = w_{start} – (w_{start} – w_{end}) * [2d/K – (d/K)^2]

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

3.2 改善：自习惯惯性权重

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

wid={wmin⁡+(wmax⁡−wmin⁡)f(xid)−fmin⁡dfaveraged−fmin⁡d,f(xid)≤faveragedwmax⁡,f(xid)>faveragedw_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.

其间：

  1. wmin,wmaxw_{min}, w_{max}是咱们预先给定的最小惯性系数和最大惯性系数，一般取0.4和0.9
  2. faveraged=∑i=1nf(xid)/nf_{average}^d = \sum_{i=1}^n f(x_i^d)/n即第d次迭代时一切粒子的平均习惯度
  3. fmind=min{f(x1d),f(x2d),…,f(xn)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 假定现在求最大值问题，那么：

wid={wmin⁡+(wmax⁡−wmin⁡)f(xid)−fmin⁡dfaveraged−fmin⁡d,f(xid)≥faveragedwmax⁡,f(xid)<faveragedw_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.

与本来的比较，现在惯性权重和迭代次数以及每个粒子习惯度有关

一个较大的惯性权值有利于大局查找

而一个较小的权值则更利于部分查找

其间：

  1. wmin,wmaxw_{min}, w_{max}是咱们预先给定的最小惯性系数和最大惯性系数，一般取0.4和0.9
  2. faveraged=∑i=1nf(xid)/nf_{average}^d = \sum_{i=1}^n f(x_i^d)/n即第d次迭代时一切粒子的平均习惯度
  3. fmind=min{f(x1d),f(x2d),…,f(xn)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 改善：主动退出迭代循环

当粒子现已找到最佳方位后，再增加迭代次数只会浪费核算时刻，那么咱们能否设计一个战略，能够主动退出迭代呢?

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

  1. 初始化最大迭代次数、计数器以及最大计数值（例如别离取100, 0, 20）
  2. 界说“函数改动量容忍度”，一般取十分小的正数，例如10^(‐6)；
  3. 在迭代的进程中，每次核算出来最佳习惯度后，都核算该习惯度和上一次迭代时最佳习惯度的改动量(取绝对值)；
  4. 判别这个改动量和“函数改动量容忍度”的相对巨细，假如前者小，则计数器加1；否则计数器清0；
  5. 不断重复这个进程，有以下两种或许:
  • ① 此刻还没有超过最大迭代次数，计数器的值超过了最大计数值，那么咱们就跳出迭代循环，查找完毕。
  • ② 此刻现已到达了最大迭代次数，那么直接跳出循环，查找完毕。
%% 自习惯权重且带有缩短因子的粒子群算法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 ,咱们要进行以下几个操作：

  1. 随机生成粒子i的邻域，邻域内总共Q个粒子(包含粒子i自身)，并找到这Q个粒子中方位最佳的那个粒子，此刻它的方针函数值最小，记其方位为lbest；
  2. 更新粒子i的速度，公式为：v = wv + c1u1*(pbest‐x) + c2u2(lbest‐x)； (注：这儿省略了一些上下标，咱们应该能了解)，这个速度公式和基本粒子群算法最大的不同在于：这儿的集体信息沟通部分运用的是邻域内的最优方位，而不是整个粒子群的最优方位；
  3. 更新粒子i的方位，公式为： x=x+ v，这个和基本粒子群算法相同；
  4. 修正方位和速度：假如粒子i的方位超过了束缚，就将其方位修正到鸿沟处；别的假如这个粒子的方位在鸿沟处，咱们还需求查看其速度是否超过了最大速度，如 果超过的话将这个速度变为0。(留意，假如是多元函数的话或许只要某个重量超过了 束缚，咱们的修正只需求针对这个重量即可)
  5. 核算粒子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