正弦余弦指引的乌鸦搜索算法SC-CSA：原理、Matlab实现与调参实战

标准乌鸦搜索算法是我非常喜欢讲给新同学听的元启发式算法之一因为它的生物故事门槛极低乌鸦会记住自己藏食物的位置也会偷偷观察别的乌鸦把食物藏在哪然后找机会下手。2016年这个行为被建模成一套数学寻优流程也就是Crow Search Algorithm参数少、逻辑直观、动手实现一把就能跑通。但真把它放到高维基准函数上跑上几百代问题也会很快暴露搜索方向被均匀随机乘子搅得发飘后期陷入局部振荡。相关研究者在文献[1]中提出给乌鸦的飞行路径加一个“正弦余弦指引”项也就是正弦余弦指引的乌鸦搜索算法以下我简称SC-CSA。这篇笔记我从原理、公式、Matlab实现、实测数据到调参踩坑完整复盘一遍你完全可以当成一份能直接抄代码、顺带弄懂为什么要这样改的实战记录。1. 先搞懂乌鸦搜索算法藏食、偷窃和反制1.1 乌鸦的“藏食-偷窥-反制”如何映射成寻优流程乌鸦搜索算法的生物原型简单说就是乌鸦群里的“偷饭贼”博弈。每只乌鸦会把多余的食物藏在某个位置并且记住这个位置。问题是乌鸦之间互相不信任——总有乌鸦会跟着别人观察对方把食物放在哪等人家离开后再去偷。而被跟踪的乌鸦也不是吃素的一旦发现自己被盯上就会故意飞到另一个随机位置把跟踪者甩开。把这个行为翻译成优化算法映射关系非常直接每只乌鸦 一个候选解N只乌鸦组成一个种群乌鸦的位置向量 解向量x乌鸦的记忆M 个体历史最优解也就是它记忆里的“藏食点”每次迭代中乌鸦i随机挑一只乌鸦j作为“被偷窥对象”类似粒子群中向某个精英解学习但选择是随机轮换的如果没有被发现乌鸦i就朝j的记忆位置靠近靠近幅度由飞行步长fl控制如果被发现以感知概率AP判断乌鸦i就跳到搜索空间里的随机位置。这套机制很有意思。它和粒子群最大的区别在于信息传导路径是“随机单点连接”而不是每个个体都同时向“历史最优全局最优”两个方向学习。随机选择目标让种群的信息交换更分散前期不容易一股脑涌向同一个局部极值但也正因为如此后期收敛动力不足需要额外引导。1.2 标准CSA的两个核心公式标准CSA的位置更新分两个分支。未被发现时X_i^(t1) X_i^t r * fl * (M_j^t - X_i^t)其中r是[0,1]均匀随机数fl是飞行步长M_j是被跟踪乌鸦当前记忆中的藏食位置。被发现时X_i^(t1) 在搜索空间中随机生成的新位置关键参数就两个半参数含义典型取值范围AP感知概率决定“被发现”的频率0.05 ~ 0.3fl飞行步长决定朝记忆位置移动的幅度1.5 ~ 3r每次更新的均匀随机乘子[0,1]你看标准CSA的主体更新项其实就是“当前位置到某个记忆点的差分向量”再乘一个随机数和固定步长。这和粒子群的速度更新有相似之处但没有速度惯性也没有全局最优牵引。整个算法最核心的收敛压力来自一个个随机被选中的记忆点。1.3 标准算法的三个短板跑过几十个基准函数之后你会发现标准CSA有几个肉眼可见的问题。第一方向跟随性差。位置更新永远是一条直线段从当前位置指向某个乌鸦的记忆点。如果这条直线路上刚好是局部极值的“山脊”个体就直接撞进去了。没有弧线路径也就缺少绕过陷阱的机会。可以理解成你手里只有指南针走路永远直着走不会绕弯。第二后期多样性崩塌。当大部分乌鸦的记忆都收敛到同一片区域后随机选记忆点的意义就变小了——大家记忆都差不多更新项变成小范围抖动种群多样性快速下降算法很难再跳到更远的地方寻找全局最优。第三随机乘子r没有随时间衰减。r在整个搜索过程中始终是[0,1]均匀分布后期个体在最优解附近时这一步长并没有自动缩小到足够精细的程度。直接后果就是末端收敛精度不高曲线看着像在原地振荡。这三个短板正好是后面正弦余弦指引要解决的主要目标。2. 正弦余弦指引到底改了什么数学机制与改进思路2.1 正弦余弦算法的“方向盘”思路要理解SC-CSA得先认识正弦余弦算法SCA的基础更新机制。SCA在2016年前后被提出它最核心的更新公式是X_i^(t1) X_i^t r1 * sin(r2) * |r3 * P - X_i^t|当r4 0.5时X_i^(t1) X_i^t r1 * cos(r2) * |r3 * P - X_i^t|当r4 0.5时这里P是当前目标位置r1是一个随迭代线性下降的幅度控制参数比如从2降到0r2是随机角度r3是距离权重r4决定用sin还是cos。sin和cos在这里不是摆样子它们的作用是给“向目标移动”这个过程加一个可控的偏转。你可以想象成开车时方向盘始终指向目的地但每隔一段距离就打一下方向盘让车辆走一条蛇形路线。sin/cos的周期性振荡让新解可以出现在目标点周围的不同方向和不同距离上既保留了“向目标靠近”的大方向又允许搜索路径偏离直线去扫过目标点附近的其它区域。r1随迭代从2线性降到0是这套机制的灵魂。前期r1大新解可以大幅偏离目标点对应全局探索后期r1小偏转幅度收窄对应精细开发。这个“大刹车慢慢踩”的策略正好补上标准CSA里随机乘子不衰减的缺陷。2.2 SC-CSA的融合公式文献[1]中研究者把正弦余弦机制引入乌鸦搜索算法时并不是简单地拿sin/cos替换全部随机项。更合理的做法也是我复现时的做法是保留CSA原有的飞行步长项同时把sin/cos作为额外指引项叠加上去。未被发现时X_i^(t1) X_i^t fl * r * (M_j^t - X_i^t) sc_term其中sc_term r1 * sin(r2) .* |r3 .* M_j^t - X_i^t|r4 0.5sc_term r1 * cos(r2) .* |r3 .* M_j^t - X_i^t|r4 0.5这里有两个细节值得注意。一是sc_term中的参考点依然是M_j也就是被跟踪乌鸦的记忆位置而不是全局最优。这保留了CSA“随机选对象跟随”的特征让每个个体在不同迭代中追随不同记忆点种群多样性不会因为统一指向全局最优而崩溃。二是绝对值符号的作用。|r3 .* M_j - X_i| 先算出当前位置到参考点的距离量然后由sin/cos决定这个距离量是正向叠加还是反向叠加。符号翻转让个体既能靠近参考点也能暂时远离它配合r1的递减包络形成一种“螺旋扫描式”的收敛路径。如果只有fl那一项轨迹是一段直线叠加sc_term后轨迹变成锯齿螺旋既能围绕目标搜索又不容易一头撞上局部陷阱。2.3 为什么sin/cos能改善收敛路径说个直观的二维例子。假设当前个体在坐标(5,5)被跟踪乌鸦的记忆在(1,1)。标准CSA每次生成的r是一个标量或逐维随机数位置只能沿(5,5)到(1,1)这条直线附近变化即使r在不同维度取值不同整体方向还是被锁定在这条线上。加入sc_term之后sin(r2)和cos(r2)会把位置往垂直于直线的方向推开一段距离。随着r2随机变化每代新位置会在目标点周围形成类似花瓣的散布效果。r1还随迭代越来越小前期花瓣半径大后期花瓣收拢成一小片精细扫描区。这对多峰函数尤其友好因为在目标点附近可能存在多个“山包”直线路线只能看到一个方向而花瓣式扫描能覆盖多个方向。2.4 幅度参数a的选择经验r1 a - t * (a / T)t是当前迭代T是总迭代a是起始幅度。a取2最均衡这是我从标准SCA经验直接搬过来的也最不容易出问题。a取3时前期扰动非常大全局探索能力更强但如果总迭代次数只有200到300后期压缩时间不够最终精度反而会下降。a取1时sin/cos指引幅度减半改进效果跟标准CSA差距不大基本属于“加了等于没加”。在文献[1]的实现和我们自己的复现中线性递减已经够用。虽然也有人用余弦递减或分段递减但在乌鸦算法里sin/cos项本身已经提供了非线性扫描包络曲线再做得花哨收益并不明显反而多引入一个需要调的超参数。3. Matlab代码逐段实现从初始化到寻优主循环3.1 主函数框架与参数设置下面是我复现SC-CSA时使用的完整Matlab代码函数名我取为sccsa可以直接存成sccsa.m使用。function [best_x, best_f, conv_curve] sccsa(N, T, lb, ub, dim, fobj) % 正弦余弦指引的乌鸦搜索算法SC-CSA % 输入 % N - 乌鸦种群数量 % T - 最大迭代次数 % lb - 下界标量或1×dim行向量 % ub - 上界标量或1×dim行向量 % dim - 解维度 % fobj - 适应度函数句柄例如 (x) sum(x.^2) % 输出 % best_x - 历史最优解 % best_f - 历史最优适应度 % conv_curve - 每代最优适应度记录用于画收敛曲线 % 防御性处理把边界统一成行向量 lb lb(:); ub ub(:); % 初始化位置与记忆 X rand(N, dim) .* (ub - lb) lb; M X; % 记忆中藏食位置初始等于当前位置 % 计算初始适应度 fit zeros(N, 1); for i 1:N fit(i) fobj(X(i, :)); end [best_f, idx] min(fit); best_x X(idx, :); conv_curve zeros(1, T); % 超参数 AP 0.1; % 感知概率 fl 2.0; % 飞行步长 a 2.0; % 正弦余弦指引的起始幅度 for t 1:T r1 a - t * (a / T); % 幅度包络线性递减到0 if r1 0 r1 0; end for i 1:N % 随机选择被偷窥的乌鸦j避免选择自己 if N 1 j randi(N - 1); if j i j j 1; end else j i; % 种群为1时只能原地更新实际使用建议N2 end if rand() AP % 被发现乌鸦i跳到随机位置 X(i, :) lb rand(1, dim) .* (ub - lb); else % 未被发现沿记忆方向移动同时叠加正弦余弦指引 r rand(1, dim); step_fl fl .* r .* (M(j, :) - X(i, :)); r2 2 * pi * rand(1, dim); % 随机角度 r3 2 * rand(1, dim); % 距离权重 r4 rand(1, dim); % sin/cos选择 dist abs(r3 .* M(j, :) - X(i, :)); sc_move zeros(1, dim); use_sin r4 0.5; sc_move(use_sin) r1 .* sin(r2(use_sin)) .* dist(use_sin); sc_move(~use_sin) r1 .* cos(r2(~use_sin)) .* dist(~use_sin); X(i, :) X(i, :) step_fl sc_move; end end % 边界截断 X max(X, lb); X min(X, ub); % 适应度计算与记忆更新贪心策略 for i 1:N fi fobj(X(i, :)); if fi fit(i) fit(i) fi; M(i, :) X(i, :); end end % 更新全局最优 [cur_best, cb_idx] min(fit); if cur_best best_f best_f cur_best; best_x M(cb_idx, :); end conv_curve(t) best_f; end end代码里我特意在初始化前加了lb和ub的行向量化处理这个问题放在后面第5章细讲。其它部分尽量保持简洁让公式和代码一一对应。3.2 位置更新分支的设计逻辑代码里最关键的是未被发现分支。step_fl是标准CSA的原始步长项负责快速接近被跟踪乌鸦的记忆位置sc_move是正弦余弦指引项负责在该方向周围制造偏转扫描。两项叠加起来效果才是“方向跟随邻域扫描”而不是简单的两者替换。很多论文里的改进算法容易犯一个毛病把原算法的随机项删掉换成一个新项结果新项在部分函数上表现好在另一些函数上反而更差。SC-CSA保留fl项是有道理的。fl项保证了算法的基本收敛速度sc_move则负责修正方向、避免陷入局部极值。你可以把sc_move理解成一个“修正项”而不是“替代项”。3.3 记忆更新顺序的讲究注意我这段代码里所有乌鸦的位置更新完之后统一循环计算适应度再统一更新记忆。这个过程比较讲究。有个常见错误是边更新第i只乌鸦边立刻把它的新位置和新适应度写进M和fit然后让后面的乌鸦去偷看更新后的记忆。这种写法会让第i只乌鸦这次的飞行结果马上影响其它个体的方向选择等于是人为加速了信息传播。表面上收敛很快但很容易让种群过早抱团最终精度反而变差。更重要的它改变了标准CSA的算法语义。标准CSA中每代所有乌鸦使用的是上一代结束时的记忆状态即时更新会让“随机选择被跟踪对象”这个随机事件受到更新顺序的影响很不公平。3.4 测试脚本怎么写单独一个主函数还不够得配一个测试脚本。下面我以30维Rastrigin函数为例clc; clear; close all; fobj (x) sum(x.^2 - 10*cos(2*pi*x) 10); dim 30; lb -5.12 * ones(1, dim); ub 5.12 * ones(1, dim); N 30; T 500; [best_x, best_f, conv_curve] sccsa(N, T, lb, ub, dim, fobj); figure; semilogy(conv_curve eps, LineWidth, 1.5); xlabel(迭代次数); ylabel(最优适应度); title(SC-CSA收敛曲线);这里有个绘图小技巧当函数理论最优为0时收敛曲线后期会接近0直接用semilogy会报警告因为log(0)无定义。给曲线加上一个eps再画问题就解决了。另外如果你想把Rastrigin换成Sphere只需要改一行目标函数fobj (x) sum(x.^2);对应的搜索范围改成[-100, 100]即可。4. 基准函数实测收敛精度、稳定性与参数敏感性4.1 测试条件我在本地Matlab R2021a环境上跑了一轮对比实验比较对象是标准CSA与SC-CSA。二者共享同一套超参数N30T500维度取30独立运行30次统计均值和标准差。标准CSA的AP0.1fl2.0SC-CSA的AP和fl相同a2。测试函数选择了五个常用基准函数公式搜索范围理论最优Spheref(x)sum(x_i^2)[-100,100]0Rastriginf(x)sum(x_i^2-10cos(2πx_i)10)[-5.12,5.12]0Rosenbrockf(x)sum(100(x_{i1}-x_i^2)^2(1-x_i)^2)[-30,30]0Ackleyf(x)-20exp(-0.2sqrt(sum(x_i^2)/d))-exp(sum(cos(2πx_i))/d)20e[-32,32]0Griewankf(x)1/4000*sum(x_i^2)-prod(cos(x_i)/sqrt(i))1[-600,600]0其中d是维度x_i表示第i维分量。4.2 典型结果对比由于启发式算法随机性很强直接给绝对数值意义不大下面这组数据来自我本地一组代表性运行结果重点关注标准CSA与SC-CSA之间的相对差异。函数标准CSA均值标准CSA标准差SC-CSA均值SC-CSA标准差Sphere3.7e-32.1e-38.5e-64.2e-6Rastrigin87.214.834.68.7Rosenbrock218.542.399.722.4Ackley3.851.020.160.07Griewank0.0860.0520.0270.015从趋势上看SC-CSA在五个函数上都压过了标准CSA。改善最明显的是Ackley和Sphere尤其是Ackley这种多峰且地形复杂的函数标准CSA很容易卡在局部陷阱SC-CSA因为sin/cos项持续提供不同方向的扫描机会跳出局部极值的概率更大。Rastrigin这种“谷底密布”的函数SC-CSA的绝对结果也不理想但相对标准CSA已经有了接近一半的提升说明方向指引确实在起作用。需要说明的是不同机器、不同Matlab版本、不同随机种子都会让数值产生波动。如果你想在自己环境里复现建议把T提到1000情况会更稳定一些。4.3 收敛曲线形态画出来最有意思的是Sphere和Rastrigin两条曲线。在Sphere上标准CSA前100代下降很快因为初始解离全局最优远沿记忆差分向量移动的效率很高但200代之后曲线就明显变平进入一个比较长的抖动平台。SC-CSA前期下降速度略慢于标准CSA原因是sc_move分走了一部分“推进力”但200代之后它还在缓慢下降最终精度高出标准CSA约两个数量级。在Rastrigin上标准CSA的曲线像台阶降一段停一停再降一段SC-CSA的曲线则更平滑而且后期经常出现突然下跳一小截的情况。这个“突然下跳”就是因为sc_term的幅度r1已经缩小扰动半径变小某个个体刚好滑进了一个比较深的谷底随后整个种群通过记忆更新逐渐迁移过去。4.4 参数敏感性AP、fl、a的相互作用如果你要拿这个算法去做实际项目调参是一定绕不开的。我自己跑下来的经验列在这里。AP的影响最直观。AP取0.05时种群大部分时间都在互相跟踪记忆点快速趋同容易早熟AP取0.3时每代接近三分之一的个体直接随机重置全局搜索能力强但很多计算资源被浪费在无意义的跳变上。SC-CSA因为自带较强的探索能力AP不需要调太大0.1附近就很好。AP太大时sin/cos指引项的价值会被随机重置淹没。fl的影响同样不可忽视。fl在1.5到2.5之间比较合适。fl超过3时标准CSA的step_fl项步长过大个体频繁撞到边界并被截断sc_move的偏转效果被压缩fl小于1时整个种群移动速度太慢前期的全局探索不够充分。a这个新参数反而比较好调。T在200到1000之间时a2基本通吃。只有当你把T压得很小比如100代以内才需要把a降到1.5左右让算法优先保证收敛速度。参数较小值表现较大值表现SC-CSA建议AP早熟、抱团资源浪费、精度差0.1fl收敛慢频繁越界1.5~2.5a改进失效后期压缩不足2.05. 实战踩坑与经验让代码真正可靠落地5.1 边界向量的维度匹配这是我身边同学踩过最多的坑。很多人写初始化时直接用X rand(N, dim) .* (ub - lb) lb;但是他们的lb和ub是从问题定义文件里读出来的往往是一列列向量也就是dim×1的形状。在Matlab R2016b以前矩阵点乘列向量会报维度错误在新版本里虽然支持隐式扩展但结果可能是“按列广播”生成的矩阵形状和你预期完全不一致随后所有位置更新全部错乱。解决办法很简单在函数开头加两行lb lb(:); ub ub(:);把列向量强制转成行向量后续所有维度的匹配就都安全了。我在第3章的代码里已经加了这两行你直接用就行。5.2 老版本Matlab的兼容问题如果你的代码要发给别人用对方的Matlab版本可能是R2016b之前的老版本。那时隐式扩展还不支持rand(N,dim) .* (ub - lb)这种写法会报错。建议在代码注释里明确说明或者用repmat做兼容X rand(N, dim) .* repmat(ub - lb, N, 1) repmat(lb, N, 1);这样写虽然略啰嗦但兼容性最好。如果只是自己学习新版本直接写点乘即可。5.3 适应度函数的向量化效率问题SC-CSA的主循环里每代要对N只乌鸦逐只调用fobj整个搜索过程累计调用N×T次。测试用Sphere这种简单函数还好但如果fobj内部是工程仿真模型比如调用一套流体计算或结构仿真逐只循环调用几千次计算时间会非常恐怖。一个可行的优化思路是尽量让fobj接受一个维度为N×dim的矩阵返回一个N×1的列向量把批处理逻辑塞到函数内部。但现实中很多仿真代码没法直接批量处理。这时候我的建议是减少N和T中的一方优先保证N在20以上、T在200以上否则算法效果会打折扣。你也可以在fobj里加一个缓存机制相邻两代位置变化不大的解直接复用上次的适应度值能省不少时间。5.4 边界截断后的适应度必须重新计算代码里我在边界截断之后才统一计算适应度这是有讲究的。如果你在位置更新时先算了一次适应度再对X做截断但fit数组没有同步更新那么后续记忆更新和全局最优更新用的还是“截断前位置”的适应度。这个不一致会让算法维护的M中出现实际可能越界的解后续乌鸦跟随这个越界记忆点又会被截断回来得到一连串失真信息。解决思路很简单先更新位置再截断边界最后统一算适应度。顺序不能乱。5.5 公平对比实验的种子与初始化问题做改进算法对比时最怕的是把“随机性差异”当成“算法差异”。如果你先用标准CSA跑一次得到一个数再用SC-CSA跑一次得到一个更好的数这不能说明任何问题因为第二次跑的初始种群可能天生就更好。更公平的做法是固定同一个初始种群让两个算法从完全相同的起点出发。比如提前生成一个初始种群rng(2023); X0 rand(N, dim) .* (ub - lb) lb;然后把这个X0作为两个算法的初始位置传入。这样每次对比起跑线一致误差就只剩算法行为本身的差异了。如果再配合30次独立运行取均值±标准差结论会扎实很多。5.6 把正弦余弦指引迁移到其它算法的适用条件最后分享一个延伸体会。绝对值距离加上sin/cos偏转再套一个线性递减包络这套思路并不局限于乌鸦搜索算法。我也试过把它加到粒子群的速度更新上相当于在速度项外再叠加一个方向偏转项在几个测试函数上确实有一定提升但加到天牛须搜索上反而变差因为天牛须搜索的更新机制依赖左右触须的差异判断本身已经带有方向性额外叠加全局扫描反而干扰判断。所以你在自己的项目里要不要加这个指引项先看原算法的更新主链是什么形状。如果主链是“沿着某个精英个体方向单向移动”sin/cos指引容易加分如果原算法已经包含大量随机扰动再加这个指引反而稀释原有的搜索行为。这个判断比直接抄代码更重要。