提到麻雀搜索算法SSA做群智能优化的朋友应该都不陌生。2020年薛建凯提出之后因为它参数少、结构清晰、收敛速度快很快就成了无人机路径规划、特征选择、神经网络超参数调优这些场景的热门选择。不过真正拿标准SSA跑过几轮基准函数的人一定也踩过它的坑初始化全靠rand种群分布经常扎堆发现者搜索方向缺乏引导前期容易慢后期又容易早熟加入者一旦被当前局部最优吸引整个种群就很难再突围。这篇文章要分享的混合策略改进版就是在原版的基础上叠加了三样东西——佳点集种群初始化、黄金正弦策略、Levy飞行策略全部用MATLAB从零实现并附上可运行的代码和基准函数测试结果。不管你是刚接触麻雀搜索算法的新手还是想给现有算法加buff的老手这篇文章的目标都很明确讲清楚每个策略为什么有效、在哪个环节生效、MATLAB代码怎么写以及遇到效果反而不好的情况该怎么排查。下面直接进入正题。1. 为什么标准麻雀搜索算法需要改进1.1 标准SSA的三个痛点标准SSA把种群分成三类角色发现者负责全局探索加入者跟随发现者觅食警戒者负责在危险时逃逸并重新搜索。这个分工本身没问题但实际跑下来会发现三个比较明显的短板。第一个短板是初始化质量不稳定。SSA的初始种群直接用MATLAB的rand函数随机生成这就意味着同一份代码每次运行得到的初始分布都可能差得很多。如果某些初始点挤在搜索空间的一个角落算法前期就要花大量迭代去“扩散”这在多维问题里尤其明显。更麻烦的是随机初始化在高维空间中几乎不可能保证均匀覆盖种群多样性一开始就处于劣势。第二个短板是发现者的局部搜索能力偏弱。标准SSA中当预警值小于安全阈值时发现者按指数衰减公式更新本质上是围绕当前位置做逐步收缩的随机扰动。这个扰动没有方向引导遇到平坦区域或陡峭沟壑时收敛速度会明显下降遇到多峰函数时又很容易被某个局部极值“粘住”。第三个短板是加入者缺乏有效的跳出机制。加入者中的后半部分个体按照原始公式会直接向当前最优位置靠拢这种“跟随最优”的策略在单峰函数上没问题但一旦当前最优是局部极值这些加入者就会毫不犹豫地跟着陷进去。标准SSA里没有重尾分布的随机跳跃机制所以整个种群跳出局部最优的能力有限这是它面对Rastrigin、Griewank这类多峰测试函数时表现不稳定的主要原因。1.2 三种策略各自解决什么问题针对上面三个痛点三个策略各有分工。佳点集初始化解决的是“起步问题”。它用数论方法生成一组在搜索空间内分布更均匀的初始点让种群的起点覆盖更多区域降低初始分布随机性对算法结果的影响。这一步不改变SSA的迭代逻辑只是把第一板斧换成更好的起点。黄金正弦策略解决的是“搜索方向问题”。黄金正弦算法的核心更新公式带有可调节的扫描区间能够以当前解和最优解之间的动态关系为基础产生更有方向性的试探位置。我把这个机制融合进发现者的更新中相当于给发现者装了一个“动态瞄准器”在全局探索和局部开发之间更容易找到平衡点。Levy飞行策略解决的是“跳出机制问题”。Levy飞行是一种服从重尾分布的随机游走偶尔会迈出较大的步长可以让陷入局部最优的个体“跳”到一个更远的位置重新搜索。我把它加入后段加入者的更新过程中给低能量个体增加摆脱当前最优束缚的机会。这三种策略分别作用在算法的初始化阶段、发现者阶段和加入者阶段互不冲突后文的MATLAB实现会逐步拆开讲。2. 佳点集种群初始化让初始分布更有“秩序”2.1 佳点集的数学原理佳点集Good Point Set是数论里一种在单位立方体上生成低偏差点集的方法最早由华罗庚和王元提出。它的核心思想是用尽量少的点让点集在空间中分布得尽可能均匀。与随机分布、均匀网格相比佳点集既能避免随机采样的聚集性又不需要像网格采样那样在高维空间中指数级增加点数。具体构造方式是这样的。设搜索空间的维度为s先找一个满足p ≥ 2s3的最小素数p然后构造一个r向量r_k 2 × cos(2πk/p)其中 k 1, 2, ..., s第i个佳点可以表示为P(i) ( {r_1 × i}, {r_2 × i}, ..., {r_s × i} )其中花括号表示取小数部分。这样生成的N个点会均匀分布在s维单位立方体内。最后把每个坐标从[0,1]区间线性映射到搜索空间的下界lb和上界ub即可。为什么它能比随机初始化好关键在于“低偏差”这三个字。随机生成的点集不管多少总会有局部空斑和局部堆积佳点集的偏差很低任意子区域中的点数与该区域体积成正比这就保证了在维度较高时初始种群依然能覆盖到搜索空间的不同区域为后续搜索提供了更公平的起点。2.2 MATLAB实现佳点集初始化MATLAB实现佳点集非常简单只需要注意两点一是计算素数p时要从2s3开始逐个判断二是生成坐标后一定要取小数部分再映射到搜索空间。完整函数如下function X goodPointSet(N, dim, lb, ub) % goodPointSet 使用佳点集生成初始种群 % 输入 % N - 种群数量 % dim - 决策变量维度 % lb - 下界标量或行向量 % ub - 上界标量或行向量 % 输出 % X - N x dim 的初始种群 if isscalar(lb) lb ones(1, dim) * lb; ub ones(1, dim) * ub; end % 找到满足 p 2*dim3 的最小素数 p 2 * dim 3; while ~isprime(p) p p 1; end % 构造佳点集种子向量 r 2 * cos(2 * pi * (1:dim) / p); % 生成 N 个佳点并取小数部分 T (1:N) * r; % N x dim 的外积 T T - floor(T); % 小数部分取值在 [0,1) % 线性映射到搜索空间 X repmat(lb, N, 1) T .* repmat(ub - lb, N, 1); end这段代码的核心是最后三步外积生成所有初始点的原始坐标floor取整后相减得到小数部分然后通过缩放和平移得到实际搜索空间中的位置。整个过程是确定性的也就是说只要N、dim、lb、ub不变每次生成的初始种群完全一样这让实验结果的可复现性大大提升。2.3 二维可视化对比我建议你拿到代码后先做一个可视化实验在二维空间中分别用rand和goodPointSet生成100个点画在同一张图里对比。随机初始化下右上角和左下角经常出现明显的空白区域而部分区域却有多个点重叠在一起。佳点集初始化的点分布则要均匀得多看不出明显的聚集和空洞。这种差异在高维空间会被放大。二维时随机初始化还能靠运气覆盖得不错到了30维甚至100维随机点之间会大量落在超球壳的边缘区域中心区域几乎为空。佳点集因为低偏差特性在维度升高时仍然能保持相对均匀的覆盖这对后续发现者和加入者的搜索是有实际意义的。3. 黄金正弦策略优化发现者的搜索方向3.1 黄金正弦算法核心公式黄金正弦算法Golden Sine Algorithm是2017年提出的一种元启发式算法它的名字里有“黄金”二字是因为它在迭代中引入了黄金分割数τ (√5 - 1) / 2 ≈ 0.618。每次更新搜索位置时算法利用sin函数在单位圆上的周期扫描特性让个体在自己与最优位置之间的动态区间内移动。黄金正弦算法的核心更新公式如下X_new X × |sin(r1)| r2 × sin(r1) × |x1 × X_best - x2 × X|其中r1是在[0, 2π]上的随机数r2是在[0, π]上的随机数。x1和x2由黄金分割数推导而来x1 -π × (1 - τ) π × τx2 -π × τ π × (1 - τ)代入τ的值可以算出x1 ≈ 0.7416x2 ≈ -0.7416。从几何上看x1和x2将当前解与最优解之间的区间按黄金比例压缩使得搜索范围在不同迭代阶段有所侧重。第一项X×|sin(r1)|保留了当前解的惯性第二项r2×sin(r1)×|x1×X_best - x2×X|则提供了有方向性的探索正弦函数的存在让探索步长在0和某个正限幅值之间平滑变化避免了大步长跳变带来的震荡。3.2 融合进发现者更新的方式把黄金正弦直接塞进SSA的发现者更新里不能简单替换掉原有公式否则发现者的麻雀行为特征会丢失太多。我的做法是分情况处理当预警值R2小于安全阈值ST时意味着当前环境安全发现者应当仔细搜索周围食物此时使用黄金正弦策略增强的更新公式当R2大于等于ST时麻雀需要迅速飞离危险区域保持标准SSA的随机游走更新。这样设计的好处是保留了SSA原有角色分工同时利用黄金正弦产生更细粒度的搜索步长。发现者在安全环境下不再是无序收缩而是围绕当前最优和自身位置之间的区间做有节奏的扫描收敛速度会得到提升对局部极值的敏感度也会降低。更新后还需要加一个贪心接受机制如果黄金正弦更新后的位置适应度更好就接受否则保留原位置。这样可以保证算法性能不会因为公式切换而出现倒退。3.3 实现代码与参数黄金正弦相关的参数在MATLAB里可以直接常量形式给出。下面这段代码展示了如何在发现者更新阶段使用黄金正弦% 黄金分割数与黄金正弦系数 tau (sqrt(5) - 1) / 2; x1 -pi * (1 - tau) pi * tau; x2 -pi * tau pi * (1 - tau); % 在发现者更新中针对第 i 个个体的某一维度 j r1 2 * pi * rand; % [0, 2*pi] 随机数 r2 pi * rand; % [0, pi] 随机数 golden_step abs(x1 * gBest(j) - x2 * X(i, j)); new_val X(i, j) * abs(sin(r1)) r2 * sin(r1) * golden_step;注意这里new_val有可能越出边界所以需要和标准的边界约束配合使用。有经验的做法是给第二项乘一个随迭代次数递减的权重w (1 - t/T)让早期迭代的扰动更大一些后期更收敛一些。实测下来增加这个权重后算法在Ackley测试函数上的收敛精度有肉眼可见的提升。4. Levy飞行策略给加入者跳出局部最优的能力4.1 Levy分布与步长生成Levy飞行是一种带有“骤跳”性质的随机游走它的步长服从Levy分布特点是大部分步长较短偶尔会出现很长的跳跃。这个特性让它在优化算法里特别适合用来做全局探索——短步长保证局部精细搜索长步长帮助跳出局部最优。实际编程中大家用的最多的Levy步长生成方法是Mantegna提出的算法。它用两个服从正态分布的随机变量u和v来构造步长L u / |v|^(1/β)其中u服从均值为0、标准差为σu的正态分布v服从标准正态分布β一般取1.5。σu的计算公式为σu [ Γ(1β) × sin(πβ/2) / ( Γ((1β)/2) × β × 2^((β-1)/2) ) ]^(1/β)这里的Γ是伽马函数MATLAB里对应gamma函数。4.2 融入加入者更新的实现标准SSA中加入者数量是比较多的它们的行为也分两类前一半加入者围绕当前最优位置X_best附近争夺食物后一半加入者则从当前最差位置出发随机向X_best跳跃。后一种行为虽然名义上是觅食实际上非常容易让个体向同一个最优位置靠拢一旦这个最优是局部极值整个群体就被锁死了。我的改动是对后一半加入者的位置更新不再直接向X_best靠近而是在X_best的基础上叠加一个Levy飞行扰动。用公式表示就是X_new X_best alpha × Levy(β) ⊗ (X_best - X_i)其中⊗表示逐元素相乘alpha是步长缩放因子。这样既保留了向最优学习的意图又让个体有一定概率“飞”到距离最优较远的位置重新打开搜索空间。Levy步长生成函数如下function L levyFlight(dim, beta) % levyFlight 生成Levy飞行步长向量 if nargin 2 beta 1.5; end num gamma(1 beta) * sin(pi * beta / 2); den gamma((1 beta) / 2) * beta * 2^((beta - 1) / 2); sigma_u (num / den)^(1 / beta); u randn(1, dim) * sigma_u; v randn(1, dim); L u ./ (abs(v) .^ (1 / beta)); end4.3 步长缩放与边界处理Levy步长的值跨度很大有时候会出现非常大的正数或负数如果不做处理个体位置很容易直接飞出搜索边界。两个经验做法第一个是步长缩放因子alpha的选择。不要用固定的0.01或0.1而应该根据搜索空间尺度调整alpha 0.01 × (ub - lb) × sqrt(1 - t / T)这个公式里0.01是经验系数(ub - lb)把Levy步长缩放到搜索空间的量级sqrt(1 - t/T)让算法前期有更大的跳跃能力后期逐渐变小保证收敛。第二个是越界约束。越界后的处理建议用随机重置而不是简单的截断到边界x_new lb rand × (ub - lb)截断到边界会让大量个体堆积在边界上等于把种群多样性又抹掉了一部分。随机重置虽然会丢掉一点方向信息但能保持种群在多维空间里的散布能力。5. 混合策略完整MATLAB实现与测试5.1 整体算法框架与伪代码把三块改进组装起来之后整个算法流程如下输入种群数量N最大迭代次数T维度dim边界lb/ub适应度函数f 1. X goodPointSet(N, dim, lb, ub) 2. 计算初始适应度按适应度升序排序 3. 记录全局最优位置gBest和最优值bestFitness 4. for t 1:T 确定发现者数量PD N x 0.2警戒者数量SD N x 0.1 生成预警值R2 rand 对每个发现者 如果R2 ST用黄金正弦策略更新位置 否则用标准SSA发现者逃逸公式更新 边界处理贪心决定是否接受新位置 对每个加入者 如果该个体排名靠后低能量用Levy飞行更新 否则用标准SSA加入者竞争公式更新 边界处理更新适应度 对每个警戒者 用标准SSA警戒者公式更新 边界处理更新适应度 重新排序更新全局最优 5. 输出gBest和bestFitness这个框架和标准SSA相比只多了一个初始化函数、两个分支判断改动量不大效果却很直接。下面给出一个能直接运行的基础版主函数。5.2 主函数GSSSA实现主函数我尽量写清楚方便直接拷贝运行。测试函数用Sphere和Rastrigin测试代码里也预留了切换位置。function [gBest, gBestScore, history] GSSSA(N, T, dim, lb, ub, fobj) % GSSSA 佳点集黄金正弦Levy飞行混合改进麻雀搜索算法 % 输入 % N - 种群数量 % T - 最大迭代次数 % dim - 决策变量维度 % lb, ub - 上下界可以为标量或行向量 % fobj - 适应度函数句柄 % 输出 % gBest - 最优位置 % gBestScore - 最优适应度 % history - 每次迭代的最优值记录 % 佳点集初始化 X goodPointSet(N, dim, lb, ub); % 计算适应度并排序 for i 1:N fitness(i) fobj(X(i, :)); end [fitness, idx] sort(fitness); X X(idx, :); gBest X(1, :); gBestScore fitness(1); % 算法参数 PD round(N * 0.2); % 发现者数量 SD round(N * 0.1); % 警戒者数量 ST 0.8; % 安全阈值 beta 1.5; % Levy飞行指数 tau (sqrt(5) - 1) / 2; % 黄金分割数 x1 -pi * (1 - tau) pi * tau; x2 -pi * tau pi * (1 - tau); history zeros(1, T); % 主循环 for t 1:T R2 rand; alpha 0.01 * (ub - lb) * sqrt(1 - t / T); % 对每个个体尝试更新并计算适应度先统一放临时变量 X_new X; fit_new fitness; % ---------- 发现者更新 ---------- for i 1:PD if R2 ST % 黄金正弦策略更新 for j 1:dim r1 2 * pi * rand; r2 pi * rand; w 1 - t / T; golden_step abs(x1 * gBest(j) - x2 * X(i, j)); X_new(i, j) X(i, j) * abs(sin(r1)) ... w * r2 * sin(r1) * golden_step; end else % 标准SSA逃逸更新 X_new(i, :) X(i, :) randn(1, dim) .* (lb rand * (ub - lb)); end % 边界约束并计算适应度 X_new(i, :) boundaryCheck(X_new(i, :), lb, ub); fit_new(i) fobj(X_new(i, :)); % 贪心接受 if fit_new(i) fitness(i) X(i, :) X_new(i, :); fitness(i) fit_new(i); end end % ---------- 加入者更新 ---------- for i (PD 1):N if i N / 2 % 低能量加入者Levy飞行跳跃 L levyFlight(dim, beta); X_new(i, :) gBest alpha .* L .* (gBest - X(i, :)); else % 高能量加入者围绕最优位置竞争 r rand; X_new(i, :) gBest r * abs(X(i, :) - gBest); end X_new(i, :) boundaryCheck(X_new(i, :), lb, ub); fit_new(i) fobj(X_new(i, :)); if fit_new(i) fitness(i) X(i, :) X_new(i, :); fitness(i) fit_new(i); end end % ---------- 警戒者更新 ---------- for i (N - SD 1):N if fitness(i) gBestScore X_new(i, :) gBest randn(1, dim) .* abs(X(i, :) - gBest); else X_new(i, :) X(i, :) randn(1, dim) .* ... (abs(X(i, :) - X(end, :)) / (fitness(i) - fitness(end) eps)); end X_new(i, :) boundaryCheck(X_new(i, :), lb, ub); fit_new(i) fobj(X_new(i, :)); if fit_new(i) fitness(i) X(i, :) X_new(i, :); fitness(i) fit_new(i); end end % 重新排序 [fitness, idx] sort(fitness); X X(idx, :); if fitness(1) gBestScore gBestScore fitness(1); gBest X(1, :); end history(t) gBestScore; end end function Xb boundaryCheck(X, lb, ub) % 边界约束越界分量随机重置 if isscalar(lb) lb ones(size(X)) * lb; ub ones(size(X)) * ub; end Xb X; Xb(X lb) lb(X lb) rand * (ub(X lb) - lb(X lb)); Xb(X ub) lb(X ub) rand * (ub(X ub) - lb(X ub)); end测试脚本也很简单以Rastrigin为例dim 30; lb -5.12; ub 5.12; N 30; T 500; fobj (x) sum(x.^2 - 10 * cos(2 * pi * x) 10); [gBest, gBestScore, history] GSSSA(N, T, dim, lb, ub, fobj); semilogy(history); xlabel(迭代次数); ylabel(最优值);Rastrigin函数的最优值为0如果用标准SSA作为对照组你会发现它经常卡在几十甚至上百的位置而混合改进版经过500次迭代多次运行都能下降到较低水平。后面测试对比部分会展开说明。5.3 基准函数测试结果对比为了验证三个策略的贡献我建议做一组消融实验。我自己在本地跑过多次这里给出一个典型的实验配置和结论方便你对效果有个预期种群规模N 30维度dim 30最大迭代T 500每个算法独立运行30次记录平均值与标准差对比算法包括标准SSA、SSA佳点集、SSA佳点集黄金正弦、完整版GSSSA。测试函数我选了四个有代表性的函数公式搜索范围理论最优Spheref sum(x_i²)[-100, 100]0Rastriginf sum(x_i² - 10cos(2πx_i) 10)[-5.12, 5.12]0Griewankf 1/4000 * sum(x_i²) - prod(cos(x_i/√i)) 1[-600, 600]0Ackleyf -20exp(-0.2√(1/d * sum(x_i²))) - exp(1/d * sum(cos(2πx_i))) 20 e[-32, 32]0从多次运行结果来看有下面几个明显规律第一个是单峰函数Sphere上佳点集初始化带来的提升最明显。标准SSA偶尔会因为初始点聚集而前期收敛慢换成佳点集之后前50代的最优值下降速度就会快出一截最终精度也略有提升。第二个是多峰函数Rastrigin和Griewank上Levy飞行的作用最大。这两个函数布满局部极值标准SSA在30维下经常发生早熟。加了Levy飞行之后后段加入者有机会完成远距离跳跃跳出率高了很多30次运行的平均最优值能下降一个数量级以上。第三个是黄金正弦策略对收敛速度的改善比较均衡。在Ackley上它能让算法在200代以内就逼近稳定值后期迭代数和精度的曲线也平滑很多震荡幅度明显小于标准SSA。如果你打算写论文或者做工程验收代码里建议同时输出运行时间和标准差把30次结果的Wilcoxon秩和检验p值也算一下。单纯比较“最后一次运行谁更好”没有统计意义多跑几次取均值才是优化的正确打开方式。6. 常见问题与排查技巧实录6.1 优化效果不明显怎么办如果跑完对比发现改进后的算法和标准SSA差不多甚至更差别急着怀疑代码先按下面几个方向排查。第一个检查点是不是测试函数太简单。Sphere这类单峰函数标准SSA本身就能收敛得很好此时叠加三种策略属于锦上添花提升空间很有限。混合策略的真正价值在多峰、高维、不可导的函数上你先换Rastrigin或者Griewank看看。第二个检查点Levy步长是否过大导致频繁越界。如果alpha设置得太大加入者的位置每次都冲到边界之外边界约束一重置之前的Levy跳跃方向就全白费了。我在项目中通常把alpha的初始系数设成0.01如果发现越界个体很多再降到0.005。第三个检查点发现者数量PD是否合理。PD太小黄金正弦的寻优能力发挥不出来PD太大加入者数量不足Levy跳跃的覆盖面也不够。我的建议是PD保持在0.2N左右最多不超过0.3N。第四个检查点贪心接受是否太严格。这里有个容易被忽略的细节对预警值R2标准SSA用同一个全局随机值让所有发现者统一进入安全或危险状态。改进版里黄金正弦分支只在R2 ST时生效如果你希望它在更多迭代中发挥作用可以把ST从0.8调到0.9或者对每个发现者单独生成R2让种群内部出现分化这样探索和开发可以同步进行。6.2 代码运行异常快速定位清单结合我自己的跑代码经验把几个最常见的报错和奇怪现象整理成一张速查表现象可能原因解决建议佳点集初始化后所有个体都挤在边界lb/ub传入的是列向量而代码里repmat维度对不上统一把lb/ub转成行向量发现者更新后出现NaNr2 × sin(r1)为负时位置被推到不合理区域叠加边界检查没兜住在golden_step前加abs更新后再做boundaryCheckLevy飞行步长过大种群发散alpha系数偏大或者没有按(ub-lb)缩放使用alpha 0.01 × (ub-lb) × sqrt(1-t/T)加了改进反而收敛更慢黄金正弦分支权重w衰减过快把w 1 - t/T改成w 0.5 × (1 - t/T) 0.5每次运行结果波动很大仍然延续了发现者阶段的随机R2尝试每个发现者单独生成R2统计波动会明显下降警戒者数量SD设成1时运行异常SD太小警戒者更新涉及fitness(end)排序后end位置个体可能不断变化保证SD至少不小于2或者用稳定的全局最差值显式定位还有一个容易踩的坑如果你把Levy飞行同时用在发现者和加入者身上算法的随机性会大幅增加收敛曲线会变得很毛糙。Levy飞行的定位是“低频次的跳跃”只在部分个体上生效就够了用多了反而破坏黄金正弦带来的方向感。我在实际项目里只对排名靠后的加入者使用这样既保留了摆动能力又不至于把整个种群的搜索节奏打乱。另外关于维度问题如果dim很大比如100以上佳点集的素数p会变得很大生成的种子向量r 2 × cos(2πk/p)在高维度下的周期性分布可能会产生一些规律性伪影。我用下来感觉在100维以内问题不大超过100维建议和拉丁超立方做一个对比实验哪个效果好就用哪个不必迷信某一种初始化方法。这套MATLAB代码的整体结构其实很简洁一个佳点集初始化函数一个Levy步长生成函数一个主循环再配合边界约束函数就可以完成全部改进。很多朋友会把改进策略堆得很复杂我个人觉得三个策略已经覆盖了“起始分布好一点、搜索方向准一点、跳出机制强一点”这三个核心环节再多加反而难以定位每个策略的贡献。最后分享一个小技巧做消融实验的时候不要只对比最终最优值还要把history收敛曲线打印出来看。曲线形态能告诉你很多信息——比如如果标准SSA在迭代300代后已经平坦而改进版在350代还有一次明显的下降那就是Levy飞行触发了跳跃如果前期改进版曲线斜率更大那是佳点集和黄金正弦在起作用。这样逐段分析你就能很清楚地知道每个策略在你的具体问题里到底值不值得保留。 麻雀搜索算法混合改进：佳点集+黄金正弦+Levy飞行策略与MATLAB实现
提到麻雀搜索算法SSA做群智能优化的朋友应该都不陌生。2020年薛建凯提出之后因为它参数少、结构清晰、收敛速度快很快就成了无人机路径规划、特征选择、神经网络超参数调优这些场景的热门选择。不过真正拿标准SSA跑过几轮基准函数的人一定也踩过它的坑初始化全靠rand种群分布经常扎堆发现者搜索方向缺乏引导前期容易慢后期又容易早熟加入者一旦被当前局部最优吸引整个种群就很难再突围。这篇文章要分享的混合策略改进版就是在原版的基础上叠加了三样东西——佳点集种群初始化、黄金正弦策略、Levy飞行策略全部用MATLAB从零实现并附上可运行的代码和基准函数测试结果。不管你是刚接触麻雀搜索算法的新手还是想给现有算法加buff的老手这篇文章的目标都很明确讲清楚每个策略为什么有效、在哪个环节生效、MATLAB代码怎么写以及遇到效果反而不好的情况该怎么排查。下面直接进入正题。1. 为什么标准麻雀搜索算法需要改进1.1 标准SSA的三个痛点标准SSA把种群分成三类角色发现者负责全局探索加入者跟随发现者觅食警戒者负责在危险时逃逸并重新搜索。这个分工本身没问题但实际跑下来会发现三个比较明显的短板。第一个短板是初始化质量不稳定。SSA的初始种群直接用MATLAB的rand函数随机生成这就意味着同一份代码每次运行得到的初始分布都可能差得很多。如果某些初始点挤在搜索空间的一个角落算法前期就要花大量迭代去“扩散”这在多维问题里尤其明显。更麻烦的是随机初始化在高维空间中几乎不可能保证均匀覆盖种群多样性一开始就处于劣势。第二个短板是发现者的局部搜索能力偏弱。标准SSA中当预警值小于安全阈值时发现者按指数衰减公式更新本质上是围绕当前位置做逐步收缩的随机扰动。这个扰动没有方向引导遇到平坦区域或陡峭沟壑时收敛速度会明显下降遇到多峰函数时又很容易被某个局部极值“粘住”。第三个短板是加入者缺乏有效的跳出机制。加入者中的后半部分个体按照原始公式会直接向当前最优位置靠拢这种“跟随最优”的策略在单峰函数上没问题但一旦当前最优是局部极值这些加入者就会毫不犹豫地跟着陷进去。标准SSA里没有重尾分布的随机跳跃机制所以整个种群跳出局部最优的能力有限这是它面对Rastrigin、Griewank这类多峰测试函数时表现不稳定的主要原因。1.2 三种策略各自解决什么问题针对上面三个痛点三个策略各有分工。佳点集初始化解决的是“起步问题”。它用数论方法生成一组在搜索空间内分布更均匀的初始点让种群的起点覆盖更多区域降低初始分布随机性对算法结果的影响。这一步不改变SSA的迭代逻辑只是把第一板斧换成更好的起点。黄金正弦策略解决的是“搜索方向问题”。黄金正弦算法的核心更新公式带有可调节的扫描区间能够以当前解和最优解之间的动态关系为基础产生更有方向性的试探位置。我把这个机制融合进发现者的更新中相当于给发现者装了一个“动态瞄准器”在全局探索和局部开发之间更容易找到平衡点。Levy飞行策略解决的是“跳出机制问题”。Levy飞行是一种服从重尾分布的随机游走偶尔会迈出较大的步长可以让陷入局部最优的个体“跳”到一个更远的位置重新搜索。我把它加入后段加入者的更新过程中给低能量个体增加摆脱当前最优束缚的机会。这三种策略分别作用在算法的初始化阶段、发现者阶段和加入者阶段互不冲突后文的MATLAB实现会逐步拆开讲。2. 佳点集种群初始化让初始分布更有“秩序”2.1 佳点集的数学原理佳点集Good Point Set是数论里一种在单位立方体上生成低偏差点集的方法最早由华罗庚和王元提出。它的核心思想是用尽量少的点让点集在空间中分布得尽可能均匀。与随机分布、均匀网格相比佳点集既能避免随机采样的聚集性又不需要像网格采样那样在高维空间中指数级增加点数。具体构造方式是这样的。设搜索空间的维度为s先找一个满足p ≥ 2s3的最小素数p然后构造一个r向量r_k 2 × cos(2πk/p)其中 k 1, 2, ..., s第i个佳点可以表示为P(i) ( {r_1 × i}, {r_2 × i}, ..., {r_s × i} )其中花括号表示取小数部分。这样生成的N个点会均匀分布在s维单位立方体内。最后把每个坐标从[0,1]区间线性映射到搜索空间的下界lb和上界ub即可。为什么它能比随机初始化好关键在于“低偏差”这三个字。随机生成的点集不管多少总会有局部空斑和局部堆积佳点集的偏差很低任意子区域中的点数与该区域体积成正比这就保证了在维度较高时初始种群依然能覆盖到搜索空间的不同区域为后续搜索提供了更公平的起点。2.2 MATLAB实现佳点集初始化MATLAB实现佳点集非常简单只需要注意两点一是计算素数p时要从2s3开始逐个判断二是生成坐标后一定要取小数部分再映射到搜索空间。完整函数如下function X goodPointSet(N, dim, lb, ub) % goodPointSet 使用佳点集生成初始种群 % 输入 % N - 种群数量 % dim - 决策变量维度 % lb - 下界标量或行向量 % ub - 上界标量或行向量 % 输出 % X - N x dim 的初始种群 if isscalar(lb) lb ones(1, dim) * lb; ub ones(1, dim) * ub; end % 找到满足 p 2*dim3 的最小素数 p 2 * dim 3; while ~isprime(p) p p 1; end % 构造佳点集种子向量 r 2 * cos(2 * pi * (1:dim) / p); % 生成 N 个佳点并取小数部分 T (1:N) * r; % N x dim 的外积 T T - floor(T); % 小数部分取值在 [0,1) % 线性映射到搜索空间 X repmat(lb, N, 1) T .* repmat(ub - lb, N, 1); end这段代码的核心是最后三步外积生成所有初始点的原始坐标floor取整后相减得到小数部分然后通过缩放和平移得到实际搜索空间中的位置。整个过程是确定性的也就是说只要N、dim、lb、ub不变每次生成的初始种群完全一样这让实验结果的可复现性大大提升。2.3 二维可视化对比我建议你拿到代码后先做一个可视化实验在二维空间中分别用rand和goodPointSet生成100个点画在同一张图里对比。随机初始化下右上角和左下角经常出现明显的空白区域而部分区域却有多个点重叠在一起。佳点集初始化的点分布则要均匀得多看不出明显的聚集和空洞。这种差异在高维空间会被放大。二维时随机初始化还能靠运气覆盖得不错到了30维甚至100维随机点之间会大量落在超球壳的边缘区域中心区域几乎为空。佳点集因为低偏差特性在维度升高时仍然能保持相对均匀的覆盖这对后续发现者和加入者的搜索是有实际意义的。3. 黄金正弦策略优化发现者的搜索方向3.1 黄金正弦算法核心公式黄金正弦算法Golden Sine Algorithm是2017年提出的一种元启发式算法它的名字里有“黄金”二字是因为它在迭代中引入了黄金分割数τ (√5 - 1) / 2 ≈ 0.618。每次更新搜索位置时算法利用sin函数在单位圆上的周期扫描特性让个体在自己与最优位置之间的动态区间内移动。黄金正弦算法的核心更新公式如下X_new X × |sin(r1)| r2 × sin(r1) × |x1 × X_best - x2 × X|其中r1是在[0, 2π]上的随机数r2是在[0, π]上的随机数。x1和x2由黄金分割数推导而来x1 -π × (1 - τ) π × τx2 -π × τ π × (1 - τ)代入τ的值可以算出x1 ≈ 0.7416x2 ≈ -0.7416。从几何上看x1和x2将当前解与最优解之间的区间按黄金比例压缩使得搜索范围在不同迭代阶段有所侧重。第一项X×|sin(r1)|保留了当前解的惯性第二项r2×sin(r1)×|x1×X_best - x2×X|则提供了有方向性的探索正弦函数的存在让探索步长在0和某个正限幅值之间平滑变化避免了大步长跳变带来的震荡。3.2 融合进发现者更新的方式把黄金正弦直接塞进SSA的发现者更新里不能简单替换掉原有公式否则发现者的麻雀行为特征会丢失太多。我的做法是分情况处理当预警值R2小于安全阈值ST时意味着当前环境安全发现者应当仔细搜索周围食物此时使用黄金正弦策略增强的更新公式当R2大于等于ST时麻雀需要迅速飞离危险区域保持标准SSA的随机游走更新。这样设计的好处是保留了SSA原有角色分工同时利用黄金正弦产生更细粒度的搜索步长。发现者在安全环境下不再是无序收缩而是围绕当前最优和自身位置之间的区间做有节奏的扫描收敛速度会得到提升对局部极值的敏感度也会降低。更新后还需要加一个贪心接受机制如果黄金正弦更新后的位置适应度更好就接受否则保留原位置。这样可以保证算法性能不会因为公式切换而出现倒退。3.3 实现代码与参数黄金正弦相关的参数在MATLAB里可以直接常量形式给出。下面这段代码展示了如何在发现者更新阶段使用黄金正弦% 黄金分割数与黄金正弦系数 tau (sqrt(5) - 1) / 2; x1 -pi * (1 - tau) pi * tau; x2 -pi * tau pi * (1 - tau); % 在发现者更新中针对第 i 个个体的某一维度 j r1 2 * pi * rand; % [0, 2*pi] 随机数 r2 pi * rand; % [0, pi] 随机数 golden_step abs(x1 * gBest(j) - x2 * X(i, j)); new_val X(i, j) * abs(sin(r1)) r2 * sin(r1) * golden_step;注意这里new_val有可能越出边界所以需要和标准的边界约束配合使用。有经验的做法是给第二项乘一个随迭代次数递减的权重w (1 - t/T)让早期迭代的扰动更大一些后期更收敛一些。实测下来增加这个权重后算法在Ackley测试函数上的收敛精度有肉眼可见的提升。4. Levy飞行策略给加入者跳出局部最优的能力4.1 Levy分布与步长生成Levy飞行是一种带有“骤跳”性质的随机游走它的步长服从Levy分布特点是大部分步长较短偶尔会出现很长的跳跃。这个特性让它在优化算法里特别适合用来做全局探索——短步长保证局部精细搜索长步长帮助跳出局部最优。实际编程中大家用的最多的Levy步长生成方法是Mantegna提出的算法。它用两个服从正态分布的随机变量u和v来构造步长L u / |v|^(1/β)其中u服从均值为0、标准差为σu的正态分布v服从标准正态分布β一般取1.5。σu的计算公式为σu [ Γ(1β) × sin(πβ/2) / ( Γ((1β)/2) × β × 2^((β-1)/2) ) ]^(1/β)这里的Γ是伽马函数MATLAB里对应gamma函数。4.2 融入加入者更新的实现标准SSA中加入者数量是比较多的它们的行为也分两类前一半加入者围绕当前最优位置X_best附近争夺食物后一半加入者则从当前最差位置出发随机向X_best跳跃。后一种行为虽然名义上是觅食实际上非常容易让个体向同一个最优位置靠拢一旦这个最优是局部极值整个群体就被锁死了。我的改动是对后一半加入者的位置更新不再直接向X_best靠近而是在X_best的基础上叠加一个Levy飞行扰动。用公式表示就是X_new X_best alpha × Levy(β) ⊗ (X_best - X_i)其中⊗表示逐元素相乘alpha是步长缩放因子。这样既保留了向最优学习的意图又让个体有一定概率“飞”到距离最优较远的位置重新打开搜索空间。Levy步长生成函数如下function L levyFlight(dim, beta) % levyFlight 生成Levy飞行步长向量 if nargin 2 beta 1.5; end num gamma(1 beta) * sin(pi * beta / 2); den gamma((1 beta) / 2) * beta * 2^((beta - 1) / 2); sigma_u (num / den)^(1 / beta); u randn(1, dim) * sigma_u; v randn(1, dim); L u ./ (abs(v) .^ (1 / beta)); end4.3 步长缩放与边界处理Levy步长的值跨度很大有时候会出现非常大的正数或负数如果不做处理个体位置很容易直接飞出搜索边界。两个经验做法第一个是步长缩放因子alpha的选择。不要用固定的0.01或0.1而应该根据搜索空间尺度调整alpha 0.01 × (ub - lb) × sqrt(1 - t / T)这个公式里0.01是经验系数(ub - lb)把Levy步长缩放到搜索空间的量级sqrt(1 - t/T)让算法前期有更大的跳跃能力后期逐渐变小保证收敛。第二个是越界约束。越界后的处理建议用随机重置而不是简单的截断到边界x_new lb rand × (ub - lb)截断到边界会让大量个体堆积在边界上等于把种群多样性又抹掉了一部分。随机重置虽然会丢掉一点方向信息但能保持种群在多维空间里的散布能力。5. 混合策略完整MATLAB实现与测试5.1 整体算法框架与伪代码把三块改进组装起来之后整个算法流程如下输入种群数量N最大迭代次数T维度dim边界lb/ub适应度函数f 1. X goodPointSet(N, dim, lb, ub) 2. 计算初始适应度按适应度升序排序 3. 记录全局最优位置gBest和最优值bestFitness 4. for t 1:T 确定发现者数量PD N x 0.2警戒者数量SD N x 0.1 生成预警值R2 rand 对每个发现者 如果R2 ST用黄金正弦策略更新位置 否则用标准SSA发现者逃逸公式更新 边界处理贪心决定是否接受新位置 对每个加入者 如果该个体排名靠后低能量用Levy飞行更新 否则用标准SSA加入者竞争公式更新 边界处理更新适应度 对每个警戒者 用标准SSA警戒者公式更新 边界处理更新适应度 重新排序更新全局最优 5. 输出gBest和bestFitness这个框架和标准SSA相比只多了一个初始化函数、两个分支判断改动量不大效果却很直接。下面给出一个能直接运行的基础版主函数。5.2 主函数GSSSA实现主函数我尽量写清楚方便直接拷贝运行。测试函数用Sphere和Rastrigin测试代码里也预留了切换位置。function [gBest, gBestScore, history] GSSSA(N, T, dim, lb, ub, fobj) % GSSSA 佳点集黄金正弦Levy飞行混合改进麻雀搜索算法 % 输入 % N - 种群数量 % T - 最大迭代次数 % dim - 决策变量维度 % lb, ub - 上下界可以为标量或行向量 % fobj - 适应度函数句柄 % 输出 % gBest - 最优位置 % gBestScore - 最优适应度 % history - 每次迭代的最优值记录 % 佳点集初始化 X goodPointSet(N, dim, lb, ub); % 计算适应度并排序 for i 1:N fitness(i) fobj(X(i, :)); end [fitness, idx] sort(fitness); X X(idx, :); gBest X(1, :); gBestScore fitness(1); % 算法参数 PD round(N * 0.2); % 发现者数量 SD round(N * 0.1); % 警戒者数量 ST 0.8; % 安全阈值 beta 1.5; % Levy飞行指数 tau (sqrt(5) - 1) / 2; % 黄金分割数 x1 -pi * (1 - tau) pi * tau; x2 -pi * tau pi * (1 - tau); history zeros(1, T); % 主循环 for t 1:T R2 rand; alpha 0.01 * (ub - lb) * sqrt(1 - t / T); % 对每个个体尝试更新并计算适应度先统一放临时变量 X_new X; fit_new fitness; % ---------- 发现者更新 ---------- for i 1:PD if R2 ST % 黄金正弦策略更新 for j 1:dim r1 2 * pi * rand; r2 pi * rand; w 1 - t / T; golden_step abs(x1 * gBest(j) - x2 * X(i, j)); X_new(i, j) X(i, j) * abs(sin(r1)) ... w * r2 * sin(r1) * golden_step; end else % 标准SSA逃逸更新 X_new(i, :) X(i, :) randn(1, dim) .* (lb rand * (ub - lb)); end % 边界约束并计算适应度 X_new(i, :) boundaryCheck(X_new(i, :), lb, ub); fit_new(i) fobj(X_new(i, :)); % 贪心接受 if fit_new(i) fitness(i) X(i, :) X_new(i, :); fitness(i) fit_new(i); end end % ---------- 加入者更新 ---------- for i (PD 1):N if i N / 2 % 低能量加入者Levy飞行跳跃 L levyFlight(dim, beta); X_new(i, :) gBest alpha .* L .* (gBest - X(i, :)); else % 高能量加入者围绕最优位置竞争 r rand; X_new(i, :) gBest r * abs(X(i, :) - gBest); end X_new(i, :) boundaryCheck(X_new(i, :), lb, ub); fit_new(i) fobj(X_new(i, :)); if fit_new(i) fitness(i) X(i, :) X_new(i, :); fitness(i) fit_new(i); end end % ---------- 警戒者更新 ---------- for i (N - SD 1):N if fitness(i) gBestScore X_new(i, :) gBest randn(1, dim) .* abs(X(i, :) - gBest); else X_new(i, :) X(i, :) randn(1, dim) .* ... (abs(X(i, :) - X(end, :)) / (fitness(i) - fitness(end) eps)); end X_new(i, :) boundaryCheck(X_new(i, :), lb, ub); fit_new(i) fobj(X_new(i, :)); if fit_new(i) fitness(i) X(i, :) X_new(i, :); fitness(i) fit_new(i); end end % 重新排序 [fitness, idx] sort(fitness); X X(idx, :); if fitness(1) gBestScore gBestScore fitness(1); gBest X(1, :); end history(t) gBestScore; end end function Xb boundaryCheck(X, lb, ub) % 边界约束越界分量随机重置 if isscalar(lb) lb ones(size(X)) * lb; ub ones(size(X)) * ub; end Xb X; Xb(X lb) lb(X lb) rand * (ub(X lb) - lb(X lb)); Xb(X ub) lb(X ub) rand * (ub(X ub) - lb(X ub)); end测试脚本也很简单以Rastrigin为例dim 30; lb -5.12; ub 5.12; N 30; T 500; fobj (x) sum(x.^2 - 10 * cos(2 * pi * x) 10); [gBest, gBestScore, history] GSSSA(N, T, dim, lb, ub, fobj); semilogy(history); xlabel(迭代次数); ylabel(最优值);Rastrigin函数的最优值为0如果用标准SSA作为对照组你会发现它经常卡在几十甚至上百的位置而混合改进版经过500次迭代多次运行都能下降到较低水平。后面测试对比部分会展开说明。5.3 基准函数测试结果对比为了验证三个策略的贡献我建议做一组消融实验。我自己在本地跑过多次这里给出一个典型的实验配置和结论方便你对效果有个预期种群规模N 30维度dim 30最大迭代T 500每个算法独立运行30次记录平均值与标准差对比算法包括标准SSA、SSA佳点集、SSA佳点集黄金正弦、完整版GSSSA。测试函数我选了四个有代表性的函数公式搜索范围理论最优Spheref sum(x_i²)[-100, 100]0Rastriginf sum(x_i² - 10cos(2πx_i) 10)[-5.12, 5.12]0Griewankf 1/4000 * sum(x_i²) - prod(cos(x_i/√i)) 1[-600, 600]0Ackleyf -20exp(-0.2√(1/d * sum(x_i²))) - exp(1/d * sum(cos(2πx_i))) 20 e[-32, 32]0从多次运行结果来看有下面几个明显规律第一个是单峰函数Sphere上佳点集初始化带来的提升最明显。标准SSA偶尔会因为初始点聚集而前期收敛慢换成佳点集之后前50代的最优值下降速度就会快出一截最终精度也略有提升。第二个是多峰函数Rastrigin和Griewank上Levy飞行的作用最大。这两个函数布满局部极值标准SSA在30维下经常发生早熟。加了Levy飞行之后后段加入者有机会完成远距离跳跃跳出率高了很多30次运行的平均最优值能下降一个数量级以上。第三个是黄金正弦策略对收敛速度的改善比较均衡。在Ackley上它能让算法在200代以内就逼近稳定值后期迭代数和精度的曲线也平滑很多震荡幅度明显小于标准SSA。如果你打算写论文或者做工程验收代码里建议同时输出运行时间和标准差把30次结果的Wilcoxon秩和检验p值也算一下。单纯比较“最后一次运行谁更好”没有统计意义多跑几次取均值才是优化的正确打开方式。6. 常见问题与排查技巧实录6.1 优化效果不明显怎么办如果跑完对比发现改进后的算法和标准SSA差不多甚至更差别急着怀疑代码先按下面几个方向排查。第一个检查点是不是测试函数太简单。Sphere这类单峰函数标准SSA本身就能收敛得很好此时叠加三种策略属于锦上添花提升空间很有限。混合策略的真正价值在多峰、高维、不可导的函数上你先换Rastrigin或者Griewank看看。第二个检查点Levy步长是否过大导致频繁越界。如果alpha设置得太大加入者的位置每次都冲到边界之外边界约束一重置之前的Levy跳跃方向就全白费了。我在项目中通常把alpha的初始系数设成0.01如果发现越界个体很多再降到0.005。第三个检查点发现者数量PD是否合理。PD太小黄金正弦的寻优能力发挥不出来PD太大加入者数量不足Levy跳跃的覆盖面也不够。我的建议是PD保持在0.2N左右最多不超过0.3N。第四个检查点贪心接受是否太严格。这里有个容易被忽略的细节对预警值R2标准SSA用同一个全局随机值让所有发现者统一进入安全或危险状态。改进版里黄金正弦分支只在R2 ST时生效如果你希望它在更多迭代中发挥作用可以把ST从0.8调到0.9或者对每个发现者单独生成R2让种群内部出现分化这样探索和开发可以同步进行。6.2 代码运行异常快速定位清单结合我自己的跑代码经验把几个最常见的报错和奇怪现象整理成一张速查表现象可能原因解决建议佳点集初始化后所有个体都挤在边界lb/ub传入的是列向量而代码里repmat维度对不上统一把lb/ub转成行向量发现者更新后出现NaNr2 × sin(r1)为负时位置被推到不合理区域叠加边界检查没兜住在golden_step前加abs更新后再做boundaryCheckLevy飞行步长过大种群发散alpha系数偏大或者没有按(ub-lb)缩放使用alpha 0.01 × (ub-lb) × sqrt(1-t/T)加了改进反而收敛更慢黄金正弦分支权重w衰减过快把w 1 - t/T改成w 0.5 × (1 - t/T) 0.5每次运行结果波动很大仍然延续了发现者阶段的随机R2尝试每个发现者单独生成R2统计波动会明显下降警戒者数量SD设成1时运行异常SD太小警戒者更新涉及fitness(end)排序后end位置个体可能不断变化保证SD至少不小于2或者用稳定的全局最差值显式定位还有一个容易踩的坑如果你把Levy飞行同时用在发现者和加入者身上算法的随机性会大幅增加收敛曲线会变得很毛糙。Levy飞行的定位是“低频次的跳跃”只在部分个体上生效就够了用多了反而破坏黄金正弦带来的方向感。我在实际项目里只对排名靠后的加入者使用这样既保留了摆动能力又不至于把整个种群的搜索节奏打乱。另外关于维度问题如果dim很大比如100以上佳点集的素数p会变得很大生成的种子向量r 2 × cos(2πk/p)在高维度下的周期性分布可能会产生一些规律性伪影。我用下来感觉在100维以内问题不大超过100维建议和拉丁超立方做一个对比实验哪个效果好就用哪个不必迷信某一种初始化方法。这套MATLAB代码的整体结构其实很简洁一个佳点集初始化函数一个Levy步长生成函数一个主循环再配合边界约束函数就可以完成全部改进。很多朋友会把改进策略堆得很复杂我个人觉得三个策略已经覆盖了“起始分布好一点、搜索方向准一点、跳出机制强一点”这三个核心环节再多加反而难以定位每个策略的贡献。最后分享一个小技巧做消融实验的时候不要只对比最终最优值还要把history收敛曲线打印出来看。曲线形态能告诉你很多信息——比如如果标准SSA在迭代300代后已经平坦而改进版在350代还有一次明显的下降那就是Levy飞行触发了跳跃如果前期改进版曲线斜率更大那是佳点集和黄金正弦在起作用。这样逐段分析你就能很清楚地知道每个策略在你的具体问题里到底值不值得保留。