1. 引力搜索算法到底在解决什么问题第一次接触引力搜索算法Gravitational Search AlgorithmGSA是在一个多峰函数优化的项目里。当时手头有个参数寻优的活儿变量维度不高但局部极值特别密集用梯度类方法很容易一头扎进某个坑里出不来用遗传算法又觉得收敛速度不够理想。后来翻文献翻到了GSA看它的物理隐喻觉得挺有意思就花了两周时间把它从公式到代码完整复现了一遍中间踩了不少坑也积累了一些文档里不会写的经验。GSA是2009年由伊朗学者Esmat Rashedi等人提出的一种启发式群体智能优化算法。它的核心灵感来自牛顿万有引力定律在宇宙中每个有质量的粒子都会吸引其他粒子质量越大吸引力越强粒子之间距离越近吸引力也越大。把这个物理规律搬到优化问题里每个候选解就相当于一个粒子粒子的质量由它的适应度值决定——解越好质量越大对其他粒子的吸引力就越强。粒子在引力作用下不断移动最终会聚集到质量最大的位置也就是全局最优解附近。这个算法能做什么简单说凡是能转化成“在一堆参数组合里找最优”的问题它都能派上用场。比如函数极值求解、神经网络权重初始化、PID控制器参数整定、特征选择、路径规划、图像分割阈值选取等等。它适合谁适合有一定编程基础、了解基本优化概念、想尝试群体智能算法的开发者或研究者。如果你之前用过粒子群算法PSO或遗传算法GA那理解GSA会非常快因为它们的整体框架很相似区别主要在粒子更新规则上。我个人的感受是GSA最大的特点是前期探索能力强因为引力作用让粒子之间保持信息交流不容易过早收敛但后期局部搜索精度有时候不如一些改进算法。所以实际用的时候往往需要根据问题特点做一些调整。下面我就把自己从零复现GSA的完整过程拆开来讲包括公式推导、参数设置、代码实现、调优技巧和踩坑记录。2. 算法核心原理与公式拆解2.1 从万有引力到优化规则牛顿万有引力公式是F G * (M1 * M2) / R^2其中G是引力常数M1和M2是两个物体的质量R是它们之间的距离。GSA把这个公式改造成了适合优化问题的形式。假设有N个粒子每个粒子在D维空间里有一个位置。粒子i对粒子j的吸引力大小是这样算的F_ij^d G(t) * (M_pi(t) * M_aj(t)) / (R_ij(t) ε) * (x_j^d(t) - x_i^d(t))这里有几个关键点需要解释。G(t)是随时间变化的引力常数它随着迭代次数增加而减小这个设计是为了让算法前期探索范围大、后期收敛更精细。M_pi是粒子i的被动引力质量M_aj是粒子j的主动引力质量在标准GSA里通常都等于粒子的惯性质量。R_ij是两个粒子之间的欧氏距离ε是一个很小的常数防止分母为零。最后那个(x_j^d - x_i^d)是方向项保证引力方向是从i指向j。粒子i在第d维上受到的总引力是所有其他粒子对它吸引力的合力。合力算出来之后根据牛顿第二定律加速度等于力除以质量a_i^d(t) F_i^d(t) / M_ii(t)然后速度和位置更新就很简单了v_i^d(t1) rand_i * v_i^d(t) a_i^d(t) x_i^d(t1) x_i^d(t) v_i^d(t1)这里的rand_i是[0,1]之间的随机数给速度更新加一点随机扰动增加探索性。2.2 质量计算的关键细节质量怎么算这是GSA里最容易被忽略但影响很大的部分。标准做法是先用适应度值算出每个粒子的相对质量。对于最小化问题先做归一化worst(t) max fit_j(t)best(t) min fit_j(t)然后算中间变量m_i(t) (fit_i(t) - worst(t)) / (best(t) - worst(t) ε)注意这里有个容易搞反的地方适应度越小解越好fit_i - worst就越小因为是负数所以m_i反而越大。这符合直觉——好解质量大。最后归一化得到质量M_i(t) m_i(t) / Σ m_j(t)这样所有粒子的质量加起来等于1每个粒子的质量在0到1之间。注意如果所有粒子适应度相同分母会变成零所以ε不能省。我一般取1e-10到1e-16之间太小了浮点精度可能出问题太大了会影响质量区分度。2.3 引力常数G的衰减策略G(t)的衰减方式直接影响算法的探索和开发平衡。标准GSA用的是G(t) G0 * exp(-α * t / T)其中G0是初始引力常数α是衰减系数t是当前迭代次数T是总迭代次数。G0一般取100α取20这是原文推荐的默认值。但我实测下来G0取100有时候前期步长太大粒子容易飞出搜索空间。对于搜索空间范围在[-1,1]的问题G0取10到50可能更稳。α取20意味着到迭代结束时G会衰减到G0 * e^-20几乎为零这个衰减速度是合适的。还有一个细节G(t)在每一维上是一样的但有些改进版本会让每一维有不同的G值这个后面讲改进的时候再说。3. 完整实操流程与代码实现3.1 环境准备与参数初始化我用Python实现依赖只有numpy和matplotlib不需要额外的优化库。先定义问题以经典的Rastrigin函数为例这是一个多峰函数全局最优在原点最优值0搜索范围通常取[-5.12, 5.12]维度取10维或30维。import numpy as np import matplotlib.pyplot as plt def rastrigin(x): n len(x) return 10 * n np.sum(x**2 - 10 * np.cos(2 * np.pi * x)) # 参数设置 N 50 # 粒子数量 D 10 # 维度 T 500 # 最大迭代次数 G0 100 # 初始引力常数 alpha 20 # 衰减系数 lb -5.12 # 下界 ub 5.12 # 上界 epsilon 1e-10 # 防零常数粒子数量N一般取20到100之间。N太小引力信息不足容易早熟N太大计算量增加但精度提升有限。我试过N30、50、80对于10维问题50个粒子已经足够。维度高的话可以适当增加到80或100。3.2 粒子初始化与主循环初始化就是均匀随机撒点np.random.seed(42) X np.random.uniform(lb, ub, (N, D)) V np.zeros((N, D)) # 记录每代最优 best_fitness_history [] best_position None best_fitness np.inf for t in range(T): # 计算适应度 fitness np.array([rastrigin(X[i]) for i in range(N)]) # 更新全局最优 current_best_idx np.argmin(fitness) if fitness[current_best_idx] best_fitness: best_fitness fitness[current_best_idx] best_position X[current_best_idx].copy() best_fitness_history.append(best_fitness) # 计算质量 worst np.max(fitness) best np.min(fitness) m (fitness - worst) / (best - worst epsilon) M m / (np.sum(m) epsilon) # 计算引力常数 G G0 * np.exp(-alpha * t / T) # 计算加速度 # 先算所有粒子对之间的距离矩阵 # 用广播计算避免双重循环 diff X[np.newaxis, :, :] - X[:, np.newaxis, :] # (N, N, D) dist np.linalg.norm(diff, axis2) epsilon # (N, N) # 计算力的大小 # F_ij G * M_i * M_j / (dist eps) * diff_ij M_i M[:, np.newaxis] # (N, 1) M_j M[np.newaxis, :] # (1, N) force_magnitude G * M_i * M_j / (dist epsilon) # (N, N) # 合力对j求和注意排除自身 # 自身距离为0force_magnitude对角线会很大需要置零 np.fill_diagonal(force_magnitude, 0) # 加速度 合力 / 自身质量 # 合力方向由diff决定 acceleration np.sum(force_magnitude[:, :, np.newaxis] * diff, axis1) / (M[:, np.newaxis] epsilon) # 速度更新 V np.random.rand(N, D) * V acceleration # 位置更新 X X V # 边界处理越界拉回 X np.clip(X, lb, ub)这段代码里有个性能优化的点值得说。如果直接用双重循环算粒子之间的力N50时每代要算2500次500代就是125万次Python循环会非常慢。用numpy广播把距离矩阵和力矩阵一次性算出来速度能快几十倍。我实测过双重循环版本跑500代要将近两分钟广播版本只要两三秒。3.3 边界处理与收敛判断边界处理看起来简单但不同策略对结果影响不小。我试过三种第一种是直接截断clip越界就拉回边界。这是最简单的但会让粒子堆积在边界上如果最优解在边界附近还好不在的话就浪费了探索能力。第二种是反弹越界后速度反向。这个实现稍麻烦但粒子不会卡在边界。第三种是随机重置越界后重新随机撒点。这个探索性最强但可能破坏已有的搜索方向。我最后选了clip因为实现简单且对大多数问题够用。如果你发现粒子大量堆在边界可以换成反弹策略。收敛判断我一般用两个条件达到最大迭代次数或者最优值连续50代没有改善。后者可以提前结束节省时间。# 收敛判断示例 no_improve_count 0 if t 0 and abs(best_fitness_history[-1] - best_fitness_history[-2]) 1e-8: no_improve_count 1 else: no_improve_count 0 if no_improve_count 50: print(f提前收敛于第{t}代) break4. 调参经验与常见问题排查4.1 参数敏感性实测记录我拿Rastrigin函数做了参数扫描固定N50、D10、T500每次只变一个参数跑20次取平均最优值。结果如下参数取值平均最优值收敛代数备注G0102.31380前期步长偏小探索不足G0500.87320比较均衡G01001.52290前期震荡大后期收敛快G02003.76260步长过大容易跳过最优alpha101.95350衰减慢后期还在探索alpha200.87320默认值表现稳定alpha402.14280衰减太快早熟N203.42400信息不足N500.87320推荐值N1000.91300提升不明显耗时翻倍从表里能看出来G050、alpha20、N50这组参数在Rastrigin上表现最好。但这不是万能公式换一个问题可能需要重新调。我的经验是搜索空间越大G0应该适当增大问题越复杂多峰越多alpha应该适当减小让算法有更多时间探索。4.2 常见问题速查表现象可能原因排查方法解决方案早熟收敛结果稳定但不好alpha太大或N太小打印每代最优值曲线减小alpha增大N震荡不收敛G0太大观察粒子轨迹减小G0或加惯性权重后期精度不够G衰减太快看最后100代最优值变化减小alpha或加局部搜索粒子飞出搜索空间G0太大或边界处理不当检查粒子位置范围减小G0改用反弹策略计算太慢用了双重循环看代码是否有for嵌套改用numpy广播结果每次不一样随机种子未固定检查是否设了seed设np.random.seed()所有粒子质量相同适应度值太接近打印m值增大epsilon或缩放适应度4.3 三个文档里不会写的避坑技巧第一个坑质量计算时如果适应度值范围很大比如从1e-5到1e5归一化后大部分粒子的m值会非常接近0只有极少数粒子质量接近1。这会导致引力几乎全部来自那几个好粒子其他粒子被忽略算法退化成局部搜索。解决办法是对适应度做对数变换或者用排名代替原始适应度来算质量。我一般用排名法适应度排名越靠前质量越大这样质量分布更均匀。第二个坑距离计算时如果某一维的量纲特别大欧氏距离会被这一维主导其他维度的信息被淹没。比如优化一个参数范围是[0, 1000]和一个范围是[0, 1]的问题距离几乎完全由第一个参数决定。解决办法是做归一化把所有维度缩放到相同范围再算距离。这个在标准GSA里没提但实际用的时候非常关键。第三个坑速度更新里的随机数rand_i是每个粒子每个维度独立的不是整个粒子共用一个。我一开始理解错了整个粒子用一个随机数结果收敛很慢。后来改成每个维度独立随机探索能力明显提升。这个细节原文公式里写的是rand_i没有下标d但实际实现时按维度独立随机效果更好。5. 改进方向与扩展应用5.1 几个实用的改进思路标准GSA有几个明显的短板后期收敛精度不够、容易在最优附近震荡、参数敏感。针对这些问题我试过几种改进效果比较明显的有三个。第一个是加入惯性权重。速度更新公式改成v_i^d(t1) w * v_i^d(t) a_i^d(t)w从0.9线性递减到0.4。这个改动很小但能明显减少后期震荡。原理和PSO里的惯性权重一样前期大w保持探索后期小w促进收敛。第二个是精英保留策略。每代把最优的几个粒子直接复制到下一代不参与速度和位置更新。这样保证最优解不会因为随机扰动而丢失。我一般保留前10%的粒子。第三个是混合局部搜索。在GSA迭代结束后以当前最优解为中心用拟牛顿法或单纯形法做一次局部精细搜索。这个对提升最终精度非常有效代价是增加一点计算量。我试过在Rastrigin上纯GSA最优值0.87加局部搜索后能到1e-6级别。5.2 实际应用场景举例除了函数优化GSA我用过两个实际场景。一个是特征选择。有个分类任务原始特征80多个想筛到20个以内。把每个粒子编码成二进制向量1表示选该特征0表示不选适应度用分类准确率减去特征数量的惩罚项。GSA跑50代左右就能找到不错的特征子集比穷举快得多。另一个是神经网络超参数调优。学习率、批大小、隐藏层节点数这些参数组合空间很大手动调很费时间。把每个粒子编码成一组超参数适应度用验证集损失GSA自动搜索。实测下来比网格搜索快5到10倍效果还略好一点。提示GSA处理离散问题时需要做编码转换二进制问题用sigmoid函数把连续值映射到0/1整数问题用round取整。转换后适应度函数要能处理非法解一般加惩罚项。5.3 和其他算法的对比和PSO比GSA的引力机制让粒子之间信息交流更充分在多峰问题上不容易早熟。但PSO的参数更少调起来更简单收敛速度通常也更快。和GA比GSA没有交叉变异操作实现更简单但GA的种群多样性保持机制更成熟。和DE比GSA的探索能力不输DE但DE的局部搜索能力更强。我的选择逻辑是问题维度低、局部极值多优先GSA问题维度高、计算资源有限优先PSO问题需要离散编码、约束复杂优先GA。当然这只是经验之谈具体还得看问题特点。6. 完整代码封装与复现建议6.1 封装成可复用类把上面的代码整理成一个类方便换问题、换参数class GSA: def __init__(self, obj_func, dim, lb, ub, N50, T500, G0100, alpha20): self.obj_func obj_func self.dim dim self.lb lb self.ub ub self.N N self.T T self.G0 G0 self.alpha alpha self.epsilon 1e-10 def optimize(self, seedNone): if seed is not None: np.random.seed(seed) X np.random.uniform(self.lb, self.ub, (self.N, self.dim)) V np.zeros((self.N, self.dim)) best_fitness np.inf best_position None history [] for t in range(self.T): fitness np.array([self.obj_func(x) for x in X]) idx np.argmin(fitness) if fitness[idx] best_fitness: best_fitness fitness[idx] best_position X[idx].copy() history.append(best_fitness) worst np.max(fitness) best np.min(fitness) m (fitness - worst) / (best - worst self.epsilon) M m / (np.sum(m) self.epsilon) G self.G0 * np.exp(-self.alpha * t / self.T) diff X[np.newaxis, :, :] - X[:, np.newaxis, :] dist np.linalg.norm(diff, axis2) self.epsilon force G * M[:, np.newaxis] * M[np.newaxis, :] / (dist self.epsilon) np.fill_diagonal(force, 0) acc np.sum(force[:, :, np.newaxis] * diff, axis1) / (M[:, np.newaxis] self.epsilon) V np.random.rand(self.N, self.dim) * V acc X X V X np.clip(X, self.lb, self.ub) return best_position, best_fitness, history用的时候gsa GSA(rastrigin, dim10, lb-5.12, ub5.12, N50, T500, G050, alpha20) pos, fit, hist gsa.optimize(seed42) print(f最优值: {fit:.6f}) print(f最优位置: {pos})6.2 复现时的几个建议第一先用低维问题验证正确性。比如1维或2维画一下粒子轨迹看看是不是真的在往最优解移动。如果轨迹乱七八糟大概率是公式实现错了。第二固定随机种子。GSA是随机算法每次结果不一样。调试的时候固定seed方便对比不同参数的效果。最终评估的时候跑20到30次取平均单次结果没有代表性。第三画收敛曲线。每代最优值画出来能直观看到收敛速度和是否早熟。正常曲线应该是前期快速下降后期平缓。如果曲线一直震荡不下降检查G衰减是否正常如果很快平了但值不好检查是否早熟。第四注意适应度函数的计算成本。GSA每代要算N次适应度500代就是25000次。如果单次适应度计算很慢比如涉及仿真总时间会很长。这种情况要么减少N和T要么用代理模型加速。6.3 后续可以扩展的方向如果你已经把标准GSA跑通了想进一步深入有几个方向可以尝试。一是二进制GSA把位置编码成0/1向量用于特征选择或背包问题。二是多目标GSA用Pareto支配关系代替单一适应度同时优化多个目标。三是并行GSA把粒子分组每组独立进化定期交流最优解适合高维问题。四是自适应参数让G0和alpha随迭代动态调整而不是固定值。我个人在实际操作中的体会是GSA的物理隐喻虽然优雅但真正决定效果的是参数调优和问题适配。不要迷信默认参数也不要指望一个算法解决所有问题。多试、多调、多对比才是正道。
