基于牛顿引力和运动定律的引力搜索算法
引力搜索算法(Gravitational Search Algorithm, GSA)是一种基于牛顿万有引力定律和运动定律的群体智能优化算法。
算法原理
GSA模拟了物理世界中的万有引力现象：
- 每个粒子(解决方案)被视为一个有质量的物体
- 质量越大表示解决方案越好
- 粒子之间通过引力相互作用，导致运动
- 随着时间的推移，质量大的粒子会吸引质量小的粒子
代码
function [bestSolution, bestFitness, convergenceCurve] = GSA(objFunc, dim, lb, ub, maxIter, popSize)
% 引力搜索算法(GSA)
% 输入:
% objFunc: 目标函数
% dim: 问题维度
% lb: 决策变量下界
% ub: 决策变量上界
% maxIter: 最大迭代次数
% popSize: 种群大小
% 输出:
% bestSolution: 最优解
% bestFitness: 最优适应度值
% convergenceCurve: 收敛曲线
% 初始化参数
G0 = 100; % 初始引力常数
alpha = 20; % 衰减系数
epsilon = 0.1; % 小常数，防止除零
% 初始化种群
positions = initialization(popSize, dim, ub, lb);
velocity = zeros(popSize, dim);
% 评估初始种群
fitness = zeros(popSize, 1);
for i = 1:popSize
fitness(i) = objFunc(positions(i, :));
end
% 记录最佳解
[bestFitness, bestIndex] = min(fitness);
bestSolution = positions(bestIndex, :);
% 初始化收敛曲线
convergenceCurve = zeros(maxIter, 1);
convergenceCurve(1) = bestFitness;
% 迭代过程
for iter = 2:maxIter
% 计算当前引力常数
G = G0 * exp(-alpha * iter / maxIter);
% 计算每个粒子的质量
[mass, M] = calculateMass(fitness, popSize);
% 计算每个粒子受到的合力
acceleration = calculateAcceleration(positions, mass, G, epsilon, iter, maxIter, popSize, dim);
% 更新速度和位置
velocity = rand(popSize, dim) .* velocity + acceleration;
positions = positions + velocity;
% 边界处理
positions = max(positions, lb);
positions = min(positions, ub);
% 评估新位置
for i = 1:popSize
fitness(i) = objFunc(positions(i, :));
end
% 更新最佳解
[currentBestFitness, currentBestIndex] = min(fitness);
if currentBestFitness < bestFitness
bestFitness = currentBestFitness;
bestSolution = positions(currentBestIndex, :);
end
% 记录收敛曲线
convergenceCurve(iter) = bestFitness;
% 显示迭代信息
if mod(iter, 100) == 0
fprintf('迭代次数: %d, 最佳适应度: %f\n', iter, bestFitness);
end
end
end
function positions = initialization(popSize, dim, ub, lb)
% 初始化种群位置
positions = zeros(popSize, dim);
for i = 1:popSize
positions(i, :) = lb + (ub - lb) .* rand(1, dim);
end
end
function [mass, M] = calculateMass(fitness, popSize)
% 计算每个粒子的质量
% 找到最差和最佳适应度值
worst = max(fitness);
best = min(fitness);
% 处理所有适应度值相同的情况
if worst == best
mass = ones(popSize, 1);
M = mass / sum(mass);
return;
end
% 计算质量
mass = (fitness - worst) ./ (best - worst);
% 归一化质量
M = mass ./ sum(mass);
end
function acceleration = calculateAcceleration(positions, mass, G, epsilon, iter, maxIter, popSize, dim)
% 计算每个粒子受到的加速度
% 确定Kbest（随着迭代逐渐减少）
Kbest = popSize * (1 - 0.9 * iter / maxIter);
Kbest = max(2, round(Kbest)); % 确保至少有两个粒子
% 按质量排序，选择最好的Kbest个粒子
[~, sortedIndices] = sort(mass, 'descend');
activeAgents = sortedIndices(1:Kbest);
% 初始化加速度矩阵
acceleration = zeros(popSize, dim);
% 计算每个粒子受到的合力
for i = 1:popSize
for j = 1:Kbest
if i ~= activeAgents(j)
% 计算粒子间的欧氏距离
R = norm(positions(i, :) - positions(activeAgents(j), :), 2);
% 计算引力
for k = 1:dim
% 引力公式: F = G * (M_i * M_j) / (R + epsilon) * (x_j - x_i)
force = G * (mass(i) * mass(activeAgents(j))) / (R + epsilon) * ...
(positions(activeAgents(j), k) - positions(i, k));
% 随机化引力方向
acceleration(i, k) = acceleration(i, k) + rand * force;
end
end
end
% 根据牛顿第二定律计算加速度: a = F / M
if mass(i) > 0
acceleration(i, :) = acceleration(i, :) / mass(i);
end
end
end
测试
使用GSA解决标准测试函数
% 定义测试函数 (Sphere函数)
sphereFunc = @(x) sum(x.^2);
% 设置参数
dim = 30; % 问题维度
lb = -100; % 下界
ub = 100; % 上界
maxIter = 1000; % 最大迭代次数
popSize = 50; % 种群大小
% 运行GSA算法
[bestSolution, bestFitness, convergenceCurve] = GSA(sphereFunc, dim, lb, ub, maxIter, popSize);
% 显示结果
fprintf('最优解: %s\n', mat2str(bestSolution, 4));
fprintf('最优适应度: %f\n', bestFitness);
% 绘制收敛曲线
figure;
plot(convergenceCurve, 'LineWidth', 2);
title('GSA收敛曲线');
xlabel('迭代次数');
ylabel('最佳适应度值');
grid on;
扩展功能：多种测试函数
% 多种测试函数定义
testFunctions = {
@(x) sum(x.^2), % Sphere函数
@(x) sum(abs(x)) + prod(abs(x)), % Schwefel 2.22函数
@(x) sum([1:length(x)] .* (x.^4)) + rand, % Quartic函数
@(x) sum(100*(x(2:end)-x(1:end-1).^2).^2 + (x(1:end-1)-1).^2) % Rosenbrock函数
};
functionNames = {'Sphere', 'Schwefel 2.22', 'Quartic', 'Rosenbrock'};
% 为每个测试函数运行GSA
results = cell(length(testFunctions), 1);
for i = 1:length(testFunctions)
fprintf('正在测试函数: %s\n', functionNames{i});
[bestSolution, bestFitness, convergenceCurve] = GSA(...
testFunctions{i}, dim, lb, ub, maxIter, popSize);
results{i}.bestSolution = bestSolution;
results{i}.bestFitness = bestFitness;
results{i}.convergenceCurve = convergenceCurve;
fprintf('%s函数的最优适应度: %f\n\n', functionNames{i}, bestFitness);
end
% 绘制所有测试函数的收敛曲线
figure;
hold on;
colors = ['r', 'g', 'b', 'm'];
for i = 1:length(testFunctions)
plot(results{i}.convergenceCurve, 'Color', colors(i), 'LineWidth', 1.5, ...
'DisplayName', functionNames{i});
end
hold off;
title('不同测试函数的GSA收敛曲线');
xlabel('迭代次数');
ylabel('最佳适应度值');
legend('show');
grid on;
参考代码 基于牛顿引力和运动定律的引力搜索算法 www.youwenfan.com/contentcng/99352.html
参数调优建议
- 种群大小(popSize): 通常设置在20-100之间，复杂问题需要更大的种群
- 最大迭代次数(maxIter): 取决于问题复杂度，一般500-5000次
- 初始引力常数(G0): 通常设置为100
- 衰减系数(alpha): 控制引力衰减速度，通常设置为20
- epsilon: 防止除零的小常数，通常设置为0.1
算法特点
- 无梯度优化: GSA不需要目标函数的梯度信息
- 全局搜索: 通过引力机制实现有效的全局搜索
- 自适应调整: 引力常数随时间衰减，平衡探索与开发
- 简单有效: 概念简单，实现容易，在许多问题上表现良好
应用领域
GSA已被成功应用于：
- 工程优化问题
- 机器学习参数优化
- 神经网络训练
- 图像处理
- 经济调度问题
这个实现提供了GSA的基本框架，您可以根据具体问题调整参数和算法细节以获得更好的性能。
浙公网安备 33010602011771号