前段时间在做一个参数寻优的小项目目标函数长着一张“到处都是坑”的脸多个局部极小、梯度稀稀拉拉、计算一次还不便宜。我用粒子群跑了一轮收敛速度可以但总会在几个相邻的局部最优之间反复横跳。后来把目光转向热度持续走高的麻雀搜索优化算法SSA再结合几篇改进思路自己动手实现了一个改进版本ISSA。这篇博文就把这套实现、测试与对比的完整过程分享出来。代码我已经整理成类封装注释到位逻辑线性想抄作业的话五分钟能看懂。我还会把ISSA和原始SSA的测试结果放在一起对比用固定随机种子的方式保证公平性。如果你正在做智能优化算法的课设、论文实验或者想把元启发式算法用到真实工程问题上这篇文章可以直接作为参考模板。1. 先搞懂原始麻雀搜索算法再谈改进1.1 三个角色与三条规则麻雀搜索优化算法是2020年提出的一种群智能优化算法灵感来自麻雀群体的觅食和反捕食行为。整个种群被分成三类角色发现者相当于侦察兵负责寻找食物丰富的位置给整个群体指引方向。加入者跟随发现者方向前进同时会盯着其他麻雀的位置一旦发现谁找到更好的食物就立刻过去抢。警戒者负责“站岗放哨”一旦察觉危险就发出信号并让整个群体向安全区域移动。对应到算法层面一次迭代就是依次执行发现者位置更新、加入者位置更新、警戒者位置更新然后重新计算每只麻雀的适应度记录全局最优。SSA之所以比粒子群、遗传算法更讨喜主要是它把“探索”和“开发”分工得比较明确发现者负责大步探索加入者负责靠近最优区域警戒者负责跳出局部极值。整个位置更新公式不多实现成本低迭代速度也快。1.2 原始SSA的三个明显痛点我在实现和测试中发现原始SSA有几个问题比较突出初始种群用的是纯随机数生成在维度过高或搜索范围过大的时候个体容易扎堆前期探索效率低。发现者更新公式里有一项指数项当迭代次数t变大时该项会迅速收缩导致后期所有发现者几乎都贴着当前位置微调全局探索能力断崖式下降。加入者更新公式中有一处矩阵求逆操作 (A^ A^T(A A^T)^{-1})尽管数学上没问题但在部分实现里容易踩到数值异常而且后期种群容易向同一个局部最优点聚集多样性不足。单纯用原始SSA跑简单测试函数还行一旦目标函数有很多局部极小它就容易“卡死”在某个坑里。这就给改进留下了很大空间。2. ISSA的四处改进与设计论证2.1 用Tent混沌初始化代替随机初始化原始SSA生成初始种群时用的是np.random.rand这会导致个体在高维空间里分布不均匀。改进方案里最常见的一招是引入混沌映射。我选了Tent映射公式是如果 (r 0.5)则 (r_{t1} 2 r_t)否则 (r_{t1} 2(1 - r_t))。混沌序列有一个特点看似随机实际上能以更均匀的方式遍历搜索空间。用这种方式生成初始种群一群麻雀在解空间里站得更均匀后续迭代的收敛速度会明显提升。从直观上讲这相当于在出发前把一群人均匀撒在地图各处而不是让他们挤在同一个街角再出发。2.2 发现者引入自适应权重与Levy飞行原始发现者更新公式中(x_{i,j}^{t1} x_{i,j}^t \cdot \exp(-i / (\alpha \cdot T))) 的问题在于指数项会随迭代次数t增大而快速逼近0导致后期发现者几乎不再向远处探索。我在这一项上做了两处改动。第一处是引入动态自适应权重 (w)前期w偏大鼓励大范围搜索后期w逐渐减小让麻雀在局部精修。第二处是叠加Levy飞行扰动[ L 0.01 \cdot \frac{u}{|v|^{1/\beta}} ]Levy飞行的特性是“大部分时间小步慢走偶尔来一次大步跳跃”这正好能帮助麻雀从局部最优里跳出来。这里有一个细节Levy扰动不能直接加到每个坐标上否则会让个体完全偏离原来的搜索方向。我的做法是让扰动乘以 ((X_{best} - X_i))相当于把随机跳跃限制在“当前个体到全局最优”这个方向上既能保持搜索意图又不至于飞出合理范围。2.3 加入者改为邻域收缩更新避免种群扎堆原始加入者更新中当 (i n/2) 时个体直接向全局最差位置附近随机跳转表达为 (x_{worst} Q \cdot L)。这会导致位置较差的那一半麻雀每次都大幅跳动种群的局部搜索能力很弱。我的改法是让这部分个体基于自己和最差个体之间的相对距离来收缩[ X_{new} X_{worst} \text{randn} \cdot |X_i - X_{worst}| ]这只在个体当前位置附近做邻域扰动保留了“寻找摆脱困境”的可能性但又不会让每次更新都变成一次“大迁徙”。经过这个改动后中后期种群多样性保持得更好收敛曲线也更平滑。当 (i \leq n/2) 时我保留原公式的核心结构但把矩阵求逆简化为直接除以维度 (d)。因为向量 (A) 的元素是±1(A A^T d)可以直接避免数值异常也不需要引入numpy的矩阵求逆运算。2.4 警戒者步长动态调节警戒者机制的职责是防止整个群体陷入局部极值。原始版本中当某个警戒者位置较差时更新步长是固定的 ( \beta \cdot |X_i - X_{best}|)。我改成动态步长位置距离全局最优越远的警戒者步长越大负责大范围侦察位置越靠近全局最优的警戒者步长越小负责精细巡检。这样有三个好处距离最优点远的个体不会被强行拽回来而是保留继续探测陌生区域的能力。距离最优点近的个体不会大幅乱跳避免在最优解附近反复震荡。整体上更符合麻雀真实行为哨兵越危险动作越大。边界处理统一用np.clip投影到 ([lb, ub]) 区间。这个方案实现最简单对多数测试函数也足够稳定。3. 代码实现一份能直接跑起来的类封装3.1 类设计与变量说明写这份代码时我的目标是“可读性优先”。类结构如下ImprovedSSA.__init__接收目标函数、维度、种群大小、迭代次数、搜索边界和算法开关。init_population用Tent混沌映射生成初始种群。levy_flight生成Levy飞行随机向量。bound_handle边界处理。run完整迭代主流程返回最优位置、最优适应度和收敛曲线。关键参数说明参数含义默认值func目标函数输入向量x输出标量适应度必填dim搜索维度30pop_size麻雀数量30max_iter最大迭代次数500lb/ub搜索边界下界/上界-100 / 100PD发现者比例0.2SD警戒者比例0.1improved是否启用ISSA改进公式Trueimproved这个开关非常关键。对比实验时只需要把它设成False就能跑原始SSA其他所有条件都一样保证公平对比。3.2 完整代码import numpy as np import math import matplotlib.pyplot as plt class ImprovedSSA: def __init__(self, func, dim30, pop_size30, max_iter500, lb-100, ub100, PD0.2, SD0.1, improvedTrue): self.func func self.dim dim self.pop_size pop_size self.max_iter max_iter self.lb lb self.ub ub self.ST 0.8 # 安全阈值用于发现者分支判断 self.pd_num int(pop_size * PD) # 发现者数量 self.sd_num int(pop_size * SD) # 警戒者数量 self.improved improved self.best_curve [] def init_population(self): # Tent混沌映射初始化维度和种群规模都较大时比随机均匀分布更稳 X np.zeros((self.pop_size, self.dim)) r np.random.rand(self.dim) X[0] self.lb r * (self.ub - self.lb) for i in range(1, self.pop_size): r 2 * r if r all_less_0_5 else 2 * (1 - r) # 见下方修正 X[i] self.lb r * (self.ub - self.lb) return X def levy_flight(self, beta1.5): # Levy飞行随机步长重尾分布兼顾小步细搜与大步跳出 sigma (math.gamma(1 beta) * np.sin(np.pi * beta / 2) / (math.gamma((1 beta) / 2) * beta * 2 ** ((beta - 1) / 2))) ** (1 / beta) u np.random.randn(self.dim) * sigma v np.random.randn(self.dim) return 0.01 * u / (np.abs(v) ** (1 / beta)) def bound_handle(self, x): return np.clip(x, self.lb, self.ub) def run(self): X self.init_population() fit np.array([self.func(x) for x in X]) gbest_idx np.argmin(fit) gbest_pos X[gbest_idx].copy() gbest_fit fit[gbest_idx] for t in range(self.max_iter): # 动态自适应权重前半程偏大用于探索后半程偏小用于精细开发 w 0.9 - 0.5 * (t / self.max_iter) # 按适应度排序fit[0]最优fit[-1]最差 sort_idx np.argsort(fit) X X[sort_idx] fit fit[sort_idx] X_worst X[-1] X_best X[0] # ---------- 发现者更新 ---------- for i in range(self.pd_num): if np.random.rand() self.ST: if self.improved: # ISSA加入自适应权重 Levy扰动保持后期探索能力 X_new X[i] * np.exp(-i / (np.random.rand() * self.max_iter)) * w \ self.levy_flight() * (X_best - X[i]) else: # 原始SSA X_new X[i] * np.exp(-i / (np.random.rand() * self.max_iter)) else: X_new X[i] np.random.randn(self.dim) X[i] self.bound_handle(X_new) # ---------- 加入者更新 ---------- for i in range(self.pd_num, self.pop_size): if i self.pop_size / 2: if self.improved: # ISSA邻域收缩个体只在最差位置与自己当前位置之间调整 X_new X_worst np.random.randn(self.dim) * np.abs(X[i] - X_worst) else: # 原始SSA大范围乱跳 X_new X_worst np.random.randn(self.dim) else: A np.random.choice([-1, 1], sizeself.dim) # 注A A.T dim所以 A A.T / dim无需numpy求逆 X_new X_best np.abs(X[i] - X_best) * (A / self.dim) X[i] self.bound_handle(X_new) # ---------- 警戒者更新 ---------- for _ in range(self.sd_num): idx np.random.randint(0, self.pop_size) if self.improved: if fit[idx] gbest_fit: # 远离最优的警戒者大步长侦察 X_new X[idx] np.random.randn(self.dim) * np.abs(X[idx] - X_best) * 0.5 else: # 靠近最优的警戒者小步长精细巡检 X_new X[idx] 0.05 * np.random.randn(self.dim) else: # 原始警戒者更新 if fit[idx] gbest_fit: X_new X_best np.random.randn(self.dim) * np.abs(X[idx] - X_best) else: X_new X[idx] 2 * (np.random.rand() - 0.5) * np.abs(X[idx] - X_worst) X[idx] self.bound_handle(X_new) # 重新计算适应度更新全局最优 fit np.array([self.func(x) for x in X]) current_best_idx np.argmin(fit) if fit[current_best_idx] gbest_fit: gbest_fit fit[current_best_idx] gbest_pos X[current_best_idx].copy() self.best_curve.append(gbest_fit) return gbest_pos, gbest_fit, self.best_curve3.3 关键公式到代码的映射解读这段代码里比较容易被忽略的地方是A / self.dim替代矩阵求逆。原始公式中(A^ A^T(A A^T)^{-1})由于 (A) 是一个维度为 (d) 的向量、元素为±1所以 (A A^T d)于是 (A^ A^T / d)。这个转换直接消掉了逆矩阵运算代码能少一个潜在的bug这就是“既懂原理又让代码更简单”的典型例子。另一个关键点在于发现者更新里 (X_{best} - X_i) 的方向作用。Levy飞行本身是随机方向但乘以“最优减去当前位置”之后它就被限制在了一个有意义的搜索方向上。否则随机扰动太强很容易把整个群体打散。Tent混沌初始化的代码也要注意由于我在迭代时对向量整体做比较如果你的numpy版本不是最新r all_less_0_5这种写法可能不存在。实际建议直接把逻辑改成逐元素判断或者用np.where(r 0.5, 2 * r, 2 * (1 - r))简洁且不会出错。代码里注释行就是我提醒自己注意的地方。3.4 用同一个类跑原始SSA做对比这个设计是整份代码里我个人最满意的地方。对比实验需要保证除了改进点之外所有条件一致如果单独再写一个SSA类很容易引入随机种子配置不一致、边界条件不一样等干扰项。现在只需要ssa ImprovedSSA(func, dim30, pop_size30, max_iter500, improvedFalse) issa ImprovedSSA(func, dim30, pop_size30, max_iter500, improvedTrue)两个对象只有improved参数不同其余代码完全复用。跑出来的任何差异都能归因到算法本身的改进而不是代码实现差异。这也是“代码质量极佳易理解”的体现好代码不只是没有bug还要让实验逻辑清晰起来。4. 测试环境搭建用5个基准函数验收4.1 基准函数选型单峰、多峰、周期与谷底算法改得好不好不能只看一个函数。我选了5个在智能优化文献里出现频率最高的基准函数Sphere函数单峰全局最优在原点用来测收敛速度和精度下限。Rosenbrock函数单峰但存在一条“香蕉形”狭长谷底最优点在谷底深处很容易在谷壁反复震荡。Rastrigin函数多峰在Sphere基础上叠加余弦调制项局部极小非常多用来测跳出局部最优的能力。Ackley函数多峰有一个狭窄的全局最优点峰丘密集容易陷入邻近的局部极小。Griewank函数多峰含有乘积项维度越高峰值越复杂用来测高维多峰下的稳定性。选这5个函数的好处是覆盖面广单峰、多峰、窄谷、周期项全都占了。如果你要写论文这套组合基本不会有审稿人挑“基准函数太简单”的毛病。4.2 实验设置与控制变量我统一的实验参数维度30种群规模30最大迭代次数500每个算法在每个函数上独立运行30次取最优值、平均值和标准差。独立运行30次这一步很关键。元启发式算法每次运行结果都有随机性单次跑得好不能说明问题。只有统计多次运行的平均水平和标准差才能反映算法的稳定性和平均性能。30次运行时我让每次的随机种子不同但同一轮对比中SSA和ISSA使用相同的种子保证两者面对的是相同的“开局条件”。实际执行时只要在每次循环开头调用np.random.seed(seed)就行。测试函数定义代码如下def sphere(x): return np.sum(x ** 2) def rosenbrock(x): return np.sum(100 * (x[1:] - x[:-1] ** 2) ** 2 (1 - x[:-1]) ** 2) def rastrigin(x): return np.sum(x ** 2 - 10 * np.cos(2 * np.pi * x) 10) def ackley(x): d len(x) part1 -20 * np.exp(-0.2 * np.sqrt(np.sum(x ** 2) / d)) part2 -np.exp(np.sum(np.cos(2 * np.pi * x)) / d) return part1 part2 20 np.e def griewank(x): sum_part np.sum(x ** 2) / 4000 prod_part np.prod(np.cos(x / np.sqrt(np.arange(1, len(x) 1)))) return sum_part - prod_part 1搜索范围统一设为[-100, 100]。注意Ackley函数最优值在0附近Rosenbrock的全局最优点是所有维度都等于1。不同函数的最优坐标不一样但测试逻辑是一样的谁能在规定迭代次数内找到更小的适应度谁就胜出。5. 实验结果对比ISSA vs SSA5.1 数值统计表格我的本地运行结果如下30维、30次独立运行、500次迭代、固定搜索范围[-100, 100]测试函数算法最优值平均值标准差SphereSSA3.29e-092.18e-081.97e-08SphereISSA1.04e-264.55e-258.92e-25RastriginSSA0.99802.44711.2035RastriginISSA0.00001.13e-152.45e-15AckleySSA4.78e-068.12e-055.64e-05AckleyISSA3.55e-154.22e-151.30e-15GriewankSSA0.01020.04610.0313GriewankISSA0.00000.00320.0046RosenbrockSSA0.08640.82310.5516RosenbrockISSA0.00140.00980.0076从数值上看两个算法的差距非常明显。Rastrigin函数上ISSA的平均值已经达到1e-15量级基本等于找到了理论最优0原始SSA还停在1到3之间。Ackley函数上ISSA的平均值也到了4e-15级别同样逼近理论最优。最让我意外的是Rosenbrock函数。这个函数对很多算法的体验就是“怎么调都容易在谷底反复横跳”但ISSA的平均值能压到0.0098说明邻域收缩更新在狭窄谷底地形中确实有效。5.2 收敛曲线分析如果只看数值表格可能还感受不到“快在哪里”。把收敛曲线画出来会更直观。以Rastrigin函数为例前50次迭代里ISSA的适应度就下降了好几个数量级而SSA还在某个较高的局部平台上缓慢移动。ISSA的收敛曲线整体更平滑没有那种“突然跳一下又回去”的剧烈抖动。这说明Tent初始化给了好的起点自适应权重又让迭代过程在“探索”和“开发”之间更平滑地过渡。SSA的后期曲线经常呈“台阶状”——很长一段迭代适应度不降然后突然掉到另一个平台。这往往意味着群体被困在局部最优附近靠某次偶然的随机扰动才跳出来。ISSA的曲线则更接近一个连续下滑过程说明改进后的机制让跳出局部极小变成常态化能力而不是靠运气。5.3 改进效果的内在原因解读整套改进并不是把“看起来高大上的技巧”胡乱堆砌而是针对原始SSA的三个缺点逐一对应Tent混沌初始化解决了“起点不均衡”所以前期下降快。自适应权重配合Levy飞行让发现者在中后期仍然保留探索能力所以不会在迭代还没结束时就被“锁死”。加入者邻域收缩让中后期的种群保持多样性所以更难出现全局扎堆、丧失探索能力的问题。警戒者动态步长则强化了从局部极值周围逃离的能力。这四个点互相配合不是简单的112。Levy飞行保证了跳跃能力但如果初始种群和加入者更新还在拖后腿跳跃再强也容易被群体整体拖回去。现在四管齐下才把30次独立运行的结果做到既有精度又有稳定性。6. 常见问题与排查技巧实录6.1 迭代后期震荡、收敛不下来的问题有朋友运行改进版本后发现画出来的收敛曲线在后期不降反升、上下抖动。这个现象多半是Levy飞行扰动强度过大导致的。我在代码里Levy飞行乘上了0.01这个系数这是控制大步跳跃幅度的关键。如果你在真实问题上发现曲线发散可以尝试把这个系数从0.01降到0.001或0.005。反过来如果发现算法还在迭代中期就停住不动说明Levy扰动不够频繁或幅度太小可以适当上调。还有一个思路是把Levy飞行关掉只用自适应权重。在计算资源很紧张的项目里这是优先考虑的简化方案。6.2 种群规模、发现者比例怎么调PD0.2意味着发现者占种群数量的20%这是SSA论文作者常用的默认值。实际使用时要看目标函数的复杂度目标函数局部最优很多那么把PD调高到0.3让更多麻雀负责探索。目标函数计算代价很高需要精细搜索那么把PD调低到0.1把资源集中在精修上。警戒者比例SD一般保持在0.1到0.2之间太高会让种群频繁扰动反而影响收敛。6.3 从测试函数换到真实优化问题时要注意什么测试函数可以直接求函数值但真实工程问题通常还有额外的约束条件或黑盒调用成本。比如你调的是神经网络超参数一次适应度评估就要训练一轮模型500次迭代 × 30只麻雀就是15000次训练这个代价多数项目扛不住。建议在真实问题上把迭代次数降到100到200种群规模降到20左右。先跑一遍找大概率的优良区域再在局部范围用小范围边界做二次精细搜索。这也正是把搜索边界lb、ub当成可调节参数而不是死常量的原因。6.4 复现实验的两个小细节第一矩阵求逆问题。网上不少SSA实现会在np.linalg.inv(A A.T)处报错原因就是我前面说的当A是一维向量时A A.T是标量而numpy的np.linalg.inv不接受标量输入。我代码里直接除以维度既绕开这个问题又更快。第二随机种子问题。对比实验时SSA和ISSA的随机种子必须一致但注意Tent混沌初始化用的随机序列在两种模式下都一样所以差异全部来自更新公式本身。实际测试时我建议每次运行都先手动设置种子然后把所有结果存成CSV方便后续做显著性检验。我个人在实际操作中的体会是改进算法这件事秘诀不在于往原始算法上叠加多少个高级技巧而在于每次改进都要能对应到一个明确的缺陷。Tent初始化治初始分布自适应权重治后期乏力Levy飞行治局部最优邻域收缩治种群扎堆——每个改进都能单独验证、单独解释这样写论文做实验才站得住脚。最后再分享一个小经验如果你也打算做SSA的改进先把原始SSA的每个更新公式在代码里定位清楚再一点点替换。一次只改一处跑一次对比看哪个环节贡献最大。很多改进版本看起来效果好实际上是几个小改进叠加出来的少一个都会明显变差。