参考文献：万有引力搜索算法的分析与改进_范炜锋

一、算法原理

引力搜索算法将所有粒子当作有质量的物体，能够作无阻力运动。每个粒子会受到解空间中其它粒子的万有引力的影响，并产生加速度向质量更大的粒子运动。由于粒子的质量与粒子的适度值相关，适度值大的粒子其质量也会更大，因此，质量小的粒子在朝质量大趋近的过程中逐渐逼近优化问题中的最优解。万有引力搜索算法的特点在于粒子不需要通过环境因素来感知环境中的情况，而是通过个体之间的万有引力的相互作用来实现优化信息的共享，因此，在没有环境因素的影响下，粒子也感知全局的情况，从而对环境展开搜索。

假设在一个D维搜索空间中包含N个粒子，第i个粒子的位置为：

其中 表示第i个粒子在第k个维度上的位置。

1.惯性质量计算

每个粒子的惯性质量是由粒子所在位置所求得的适应值有关，在时刻t，粒子Xi的质量用Mi(t)来表示。由于惯性质量根据其相应的适应度值的大小来计算，因此，惯性质量越大的粒子表明越接近于解空间中的最优解，对其它粒子的吸引力越大。质量Mi(t)的计算公式如下：

其中 表示粒子Xi的适应度函数。best(t)表示时刻 t中的最佳解，worst(t)表示时刻t中的最差解，如果是求最小值的问题，best对应粒子适应度的最小值，worst对应最大值；如果是求最大值问题，那么正好相反。

2.引力计算

在时刻t，物体j在第k维上受到物体i的引力如下：

其中: 表示一个非常小的常量; Maj(t)表示粒子j的惯性质量，而Mpi(t)表示粒子i的惯性质量。G(t)表示 t时刻的万有引力常数，会随着时间增加而会变小，计算公式如下：

Rij表示粒子i和粒子j的欧式距离，计算公式如下：

综上，在t时刻，第k个维度上作用于粒子i的合力的计算公式为：

另外，为了增强算法的收敛速度，还可以在计算粒子i受到的合力时，只考虑适应度值排名靠前粒子的影响，但这样有可能会导致算法陷入局部最优。

3.位置更新

当粒子受到其它粒子的引力作用后就会产生加速度，则粒子i在第k个维维度上获得的加速度为其作用力与惯性质量的比值:

在每一次迭代过程中，粒子根据计算得到的加速度来更新速度和位置，更新方式如下：

4.算法步骤

标准引力搜索算法的具体步骤如下:

1.初始化算法中的所有粒子的位置与加速度，并设置迭代次数与算法中的参.

数。

2.对每个粒子计算该粒子的适应值，并更新重力常数。

3.由计算得到的适应值计算每个粒子的质量，进一步计算每个粒子的加速度。

4.根据公式更新每个粒子的速度，然后更新粒子的位置。

5.如果未满足终止条件，返回步骤2；否则，输出此次算法的最优解。

二、测试函数

1.Sphere函数

其中x的取值范围为[-5.12，5.12]，最优解在[0 0...0]处取得，最优值为0。

function fitness=Sphere(pop)
    %取值范围[-5.12，5.12]，最优解在[0 0...0]处取得，最优值为0
    fitness=sum(pop.^2);
end

2.Griewank函数

其中，x的取值范围在[-600，600]，最优解在[0 0...0]处取得，最优值为0。

function fitness=Griewank(pop)
    %取值范围[-600，600]，最优解在[0 0...0]处取得，最优值为0
    fitness=sum(pop.^2)/4000-prod(cos(pop)./sqrt(1:length(pop)))+1;
end

3.Rosenbrock函数

其中x的取值范围[-5，10]，最优解在[1 1...1]处取得，最优值为0。

function fitness=Rosenbrock(pop)
    %取值范围[-5，10]，最优解在[1 1...1]处取得，最优值为0
    fitness=sum((pop(1:end-1).^2-pop(2:end)).^2+(pop(1:end-1)-1).^2);
end

4.Ackley函数

其中x的取值范围在[-15，30]，最优解在[0 0...0]处取得，最优值为0。原文中有个笔误，没有列出c的取值，c=20。

function fitness=Ackley(pop)
    %取值范围[-15，30]，最优解在[0 0...0]处取得，最优值为0
    a=-0.2*sqrt(mean(pop.^2));
    b=mean(cos(20*pop));
    fitness=20+exp(1)-20*exp(a)-exp(b);
end

5.Rastrign函数

其中x的取值范围在[-5.12，5.12]，最优解在[0 0...0]处取得，最优值为0。

function fitness=Rastrign(pop)
    %取值范围[-5.12，5.12]，最优解在[0 0...0]处取得，最优值为0
    fitness=sum(pop.^2-10*cos(2*pi*pop)+10);
end

三、matlab部分代码

%% 清理内存空间
clc
clear
close all
%% 粒子参数的设定
index=input('请选择测试函数：1-Sphere，2-Griewank，3-Rosenbrock，4-Ackley，5-Rastrign');
pop_num=100;%粒子数量
dim=2;%问题的维度/决策变量的个数
[x_min,x_max,v_min,v_max]=set_pop(index,dim);%设置粒子位置，速度的上下限
iteration_num=200;%最大迭代次数
fitness=zeros(1,pop_num);%各粒子的适应度函数值
worst=zeros(1,iteration_num);%粒子的最差适应度
best=zeros(1,iteration_num);%粒子的最佳适应度
%% 初始化种群
pop=(x_max-x_min).*rand(pop_num,dim)+x_min;%随机初始化粒子位置
pop_v=zeros(pop_num,dim);%初始化粒子速度
Dij=zeros(pop_num,pop_num);%粒子i和j之间的距离
Fij=zeros(pop_num,pop_num);%粒子i和j之间的引力
F=zeros(pop_num,dim);%粒子受到的合力
Alpha=zeros(pop_num,dim);%粒子的加速度
%% 算法参数的设定
% G0:引力常量
% r:引力公式中的常量
% K:只取适应度前20%的粒子的引力
% dt:每次迭代的时间
[G0,r,K,dt]=deal(100,1,0.2,2);
iteration=1;
%% 迭代求最优解
while iteration<=iteration_num
    省略......
    iteration=iteration+1;
end
figure(2)
plot(1:iteration_num,best)
title('算法收敛情况')
xlabel('迭次次数')
ylabel('适应度函数')

四、测试结果

测试图像中始终有一个粒子不会收敛，这是算法中设置的一个全局搜索粒子，目的是避免算法陷入局部最优。为了便于迭代过程的可视化，设定粒子的维度为2，五个测试函数的迭代结果如下：

1.Sphere函数

2.Griewank函数

3.Rosenbrock函数

4.Ackley函数

5.Rastrign函数

结果表明，针对上面五个测试函数，采用基本的引力搜索算法就可以迭代求得最优解，而且最终所有粒子都会因引力的作用固定到全局最优的位置上。