★ FEATURED ARTICLE 差分进化算法Differential Evolution, DE是我在整理单目标连续优化测试问题时几乎必测的算法。它不要求目标函数可导不依赖梯度信息几十行Matlab就能写完但性能却经常能跟一些复杂算法掰手腕。很多人第一次接触DE都是拿Sphere函数来跑因为Sphere是教科书级的基准函数表达式简单、单峰、全局最优已知。今天这篇东西就围绕“差分进化算法matlab源码(单目标sphere函数)”这个主题把一份能直接跑的源码拿出来逐段拆解顺便把我在调参时踩过的坑也一起写清楚。如果你是刚学进化计算的新手或者正需要用DE做baseline对比、拿它当入门级求解器用这篇文章应该能省你不少事儿。我尽量把每个算子的“为什么这么写”都讲透而不是只贴一份没有注释的代码。1. 问题定义与算法设计思路1.1 为什么拿Sphere函数当试金石Sphere函数的数学形式是 f(x)sum_{i1}^D x_i^2全局最小值在原点 x0 处最优值 f(0)0。它没有任何局部最优陷阱或者说它压根就是一个“标准答案”明摆着的函数。你可能会觉得拿这么简单的函数测算法特别没意思但简单函数恰恰是最重要的试金石一个DE实现如果连Sphere都优化不到足够小的精度那代码里多半有逻辑错误反过来如果能顺利收敛到接近0至少说明变异、交叉、选择这条流水线是对的。在单目标优化领域Sphere函数通常会被放在benchmark测试集的第一位很多论文用“Sphere DE”作为入门组合。它的搜索范围一般定义为[-100,100]维度可以从10、30一直测到1000。搜索范围大、最优在中心这对算法是一种很好的压力测试。如果初始化的时候个体没有覆盖到原点附近DE能不能靠差分变异自己“游”过去这是检验算法探索能力的第一步。对于初学者来说Sphere也是一个绝佳的调试工具无论你的代码怎么改最优适应度都应该稳定下降如果曲线不动了那你就能很快定位是变异还是交叉出了问题。1.2 DE的算法流程初始化、变异、交叉、选择DE的标准流程可以拆成四个步骤。第一步是初始化种群随机生成NP个个体每个个体就是一个D维向量对应一个候选解。第二步是变异对当前个体i从种群中随机挑出三个互不相同的个体r1、r2、r3用公式 V X_r1 F * (X_r2 - X_r3) 生成一个变异向量。这个扰动方向完全是随机的但又被种群分布约束着所以它既有点“盲目”又有点“智能”。第三步是交叉把变异向量V和当前个体X_i按维度混合得到一个试验向量U。第四步是选择比较U和X_i的适应度谁更好谁就存活到下一代。循环执行这四步直到达到最大迭代次数或者适应度降到预设阈值。整个过程完全不需要求导也不需要知道目标函数的长相只需要能计算出每个候选解的适应度。也正因为如此DE在工程优化里非常受欢迎目标函数可能是黑箱甚至带有噪声但只要有个评分函数DE就能跑起来。1.3 为什么“差分”变异能自动调整步长差分变形的核心是 V X_r1 F * (X_r2 - X_r3)这个式子里 X_r2 - X_r3 是种群中两个个体的间距乘以F之后作为扰动步长。有意思的是这个步长会随着种群逐渐收敛而自动变小早期种群分散差分向量大变异扰动也大搜索范围自然就大后期种群集中到最优区域附近个体之间差异变小差分向量也变小变异扰动就自动收窄相当于步长自适应衰减。这个特性和梯度下降里的学习率衰减有点像但DE是纯粹从种群状态里“长”出来的不需要额外设计衰减策略。它也不会像梯度下降那样刚起步就陷入局部点因为它是群体里的随机差分天然带有随机性。这一点是DE最容易被人忽略的优势。2. matlab源码实现与核心算子详解2.1 完整源码与运行方式下面这份Matlab源码是主脚本形式变量设置、主循环、画图全在同一个文件里。我把每一段都写了注释方便你在自己的机器上跑通。%% 差分进化算法求解单目标Sphere函数 clear; clc; close all; % ---------------- 参数设置 ---------------- NP 50; % 种群规模 D 30; % 决策变量维度 F 0.5; % 缩放因子 CR 0.9; % 交叉概率 MaxIt 1000; % 最大迭代次数 lb -100; % 变量下界 ub 100; % 变量上界 % ---------------- 初始化 ---------------- X lb rand(NP, D) .* (ub - lb); % 随机生成初始种群 fit sum(X.^2, 2); % Sphere函数适应度 [best_fitness, best_index] min(fit); best_solution X(best_index, :); convergence zeros(1, MaxIt); % 记录每代的最优适应度 % ---------------- 主循环 ---------------- for it 1:MaxIt for i 1:NP % 1. 从种群中随机选择三个互不相同且不等于i的个体 index_pool setdiff(1:NP, i); rp randperm(length(index_pool)); r index_pool(rp(1:3)); r1 r(1); r2 r(2); r3 r(3); % 2. 变异DE/rand/1 V X(r1, :) F * (X(r2, :) - X(r3, :)); % 边界处理越界量直接拉到边界 V max(V, lb); V min(V, ub); % 3. 交叉二项式交叉 U X(i, :); jrand randi(D); % 保证至少一维来自变异向量 for j 1:D if rand CR || j jrand U(j) V(j); end end % 4. 选择贪心策略 f_new sum(U.^2, 2); if f_new fit(i) X(i, :) U; fit(i) f_new; end end % 记录全局最优 [current_best, best_index] min(fit); if current_best best_fitness best_fitness current_best; best_solution X(best_index, :); end convergence(it) best_fitness; end % ---------------- 输出结果 ---------------- fprintf(最优解\n); disp(best_solution); fprintf(最优适应度%.2e\n, best_fitness); % ---------------- 画出收敛曲线 ---------------- figure; semilogy(convergence); grid on; xlabel(迭代次数); ylabel(最优适应度); title(DE优化Sphere函数收敛曲线);把这段代码保存成de_sphere.m直接运行就能看到结果。convergence数组记录了每一代的全局最优适应度最后用semilogy画出来因为适应度变化经常跨越多个数量级用对数坐标才看得清楚。如果你用的是旧版Matlabrandperm(n,3)这种写法可能不支持所以我特意用了先randperm(length(index_pool))再取前三个的写法兼容性更好。2.2 变异算子随机三选一的细节变异是DE区别与遗传算法的最大特点。标准DE/rand/1要求 r1、r2、r3三个个体互不相同并且都不能等于当前个体i。为什么一定要不等于i因为如果 r2 等于 r3X_r2 - X_r3 就会变成零向量变异就失效了如果 r1 等于 i那变异向量会以当前个体为基础加一个扰动可能会削弱种群的探索能力虽然某些DE变体会故意这么设计但标准DE/rand/1不推荐。代码里用setdiff(1:NP, i)把i从候选池中剔除再用randperm随机打乱并取前三个索引这样能保证选取的个体一定满足约束。F的取值通常在0.4到0.9之间。F太大会让变异向量经常跳出搜索边界虽然边界处理会拉回来但大量个体被压在边界上会影响收敛F太小会让差分扰动不足种群慢慢变成“一潭死水”。我自己的经验是在Sphere这类光滑单峰问题上F0.5是一个很稳的起点但换到复杂多峰函数时可能需要适当调大。2.3 交叉算子为什么必须保证至少一维来自变异交叉的目的是把变异向量和当前父代个体做信息混合产生试验向量。二项式交叉的做法是遍历每一维按概率CR决定这一维是取变异向量的值还是父代的值。这里有一个极其容易漏掉的细节必须强制至少一个维度来自变异向量代码里的jrand randi(D)就是干这个的。如果少了这个条件当随机数恰好都大于CR时试验向量会完全等于父代个体这次变异等于白做了选择步骤也无法更新个体白白浪费一次适应度计算。交叉概率CR决定了试验向量中有多少比例的维被“刷新”。CR偏低时每次只有少数维度改变收敛会比较慢但对变量独立的问题来说反而稳定CR偏高时大量维度一起改变搜索步幅更大适合变量之间有关联的问题。Sphere函数的各个维度彼此独立所以用较高的CR通常效果不错我在代码里默认设了0.9。如果你想感受CR的影响可以把代码里的CR从0.1到0.9都跑一遍观察收敛曲线会非常直观。2.4 选择算子的“贪心”逻辑选择这一步最简单但也很重要。代码写的是if f_new fit(i)也就是说只要试验向量的适应度不差于当前个体就替换。为什么不写从纯贪心的角度和都能保证最优适应度不增但允许在适应度相等时用新个体替换旧个体这能给种群带来一点位置上的随机性避免某些个体长期驻留。在稀疏问题里这种细微的差别可能影响多样性所以标准DE一般写。另一个容易犯的错是只更新了X(i,:)却忘了更新fit(i)。下一次迭代再用到这个个体的适应度时拿到的还是旧值整个算法的逻辑就乱了。我见过不少新手在这里踩坑最后收敛曲线来回抖动。所以代码里刻意把fit(i)和X(i,:)绑定更新这个习惯非常值得养成。3. 参数设置与实验观察3.1 三个关键参数的作用范围DE核心参数只有三个种群规模NP、缩放因子F、交叉概率CR。虽然号称参数少但每个参数对性能的影响都不小。我整理了一个常用范围表方便直接对照。参数常用范围推荐起点作用NP5D ~ 10D50D30时控制种群多样性F0.4 ~ 0.90.5控制差分扰动幅度CR0.3 ~ 0.90.9控制交叉维度比例NP太小会让种群多样性不足容易早熟NP太大会增加每一代的适应度计算量收敛速度并不一定变快。对于D30的Sphere问题NP30~100都算是合理区间。F和CR的搭配也很讲究我见过有些优化问题F取0.8、CR取0.2效果出奇地好所以不要迷信默认值要拿自己的问题去测。一个比较稳的调试流程是先固定NP然后分别在F的网格和CR的网格上各跑几次用收敛曲线和最终精度来选参数。3.2 维度变化对收敛的影响我实际跑的时候D从10升到30再到50会发现同样的NP和迭代次数下高维问题的收敛速度明显变慢。D10时几百代已经能收敛到1e-12以下D30时往往需要上千代才能到类似精度D50时想达到同样的精度就需要更多迭代或者更大NP。这个现象的本质是搜索空间体积随着维度指数增长也就是常说的“维度灾难”。因此如果你需要把这个问题改造成高维场景建议把最大迭代次数和维度挂钩比如“迭代次数至少取500到1000”起步再根据收敛曲线判断是否要继续。Sphere函数的理论最优值是0但浮点数计算下通常只能逼近到1e-14左右就到底了这是机器精度导致的不是算法的问题。看收敛曲线时如果最后一大段是平直线说明已经达到当前种群条件下的极限。3.3 收敛曲线怎么读、怎么判断收敛用semilogy画出来的曲线一眼就能看出算法的收敛行为。正常情况下适应度会先快速下降然后慢慢变缓最终趋于平直。如果曲线呈现平滑下降说明变异和交叉配合得比较顺如果曲线呈阶梯状说明每次只降低一点可能是个体更新频率较低可以适当提高CR或F。如果曲线在后半段突然再次明显下降那往往是某个个体偶然跳到了更好的区域这说明算法仍然在探索还没完全收敛。反过来说如果曲线一开始下降很快随后立刻变成一条平线而最终适应度离0还很远那多半是参数不合适或边界处理有问题。这时候不要急着加迭代次数而是先检查F是否太小、NP是否太低。还有一种情况是曲线不断上下震荡这通常说明F过大或者CR过高导致种群一直在最优区域附近“反复横跳”适当调小F会让曲线稳定下来。4. 常见问题与调试技巧4.1 越界处理的三种方案在变异步骤之后V可能超出搜索边界。上面代码用的是边界吸收也就是把越界的分量直接钳制到边界上。这种办法简单直接对Sphere这种最优在搜索范围内部的问题效果很好因为边界吸收不会把个体推离原点太远。但如果最优值本身就是搜索边界边界吸收会导致大量个体挤在边界上种群多样性下降。处理方式实现思路适用场景边界吸收越界值直接设为边界最优在搜索范围内部时稳定随机重置越界后重新在边界内均匀随机取值需要保持种群多样性时反射越界后按镜面方式反弹回边界内想要保留差分方向信息要做随机重置只需要把V max(V, lb); V min(V, ub);换成V lb rand(size(V)) .* (ub - lb);。要做反射需要逐维判断并计算。反射有一个小坑如果跨步太大反射后可能还是越界需要再处理一次。我在调试复杂测试函数时最常用的是边界吸收加随机重置结合这样实现简单多样性也不会太差。4.2 早熟与多样性的调整早熟是进化算法的老问题。DE的差分变异天然有步长自适应的特性但也会在种群高度集中后失去探索能力。比如所有个体都挤在某个区域内个体之间的差分向量趋近于零变异就成了微扰这时候算法几乎只能做局部精修。如果当前区域不是全局最优那就成了“早熟收敛”。遇到这种情况最直接的调整是增大NP和F。NP变大能维持更多样性F变大能让差分向量注入更强的扰动。还有一个很有效的策略是换用DE/rand/2也就是把变异公式改成 V X_r1 F*(X_r2-X_r3) F*(X_r4-X_r5)。两次差分叠加可以得到更大范围的扰动搜索能力更强。代价是每一代需要多选两个个体计算量略增但通常能明显改善早熟问题。代码只需要把变异部分的选取从3个个体改成5个个体其他逻辑不变有兴趣的人可以自己改一下试试。4.3 代码运行效率的优化技巧上面这份源码为了教学清晰用了两重循环每代里对每个个体还要做D维交叉循环。在Matlab里这种循环方式跑D30、MaxIt1000是完全没问题的但如果你想把D拉到500甚至1000循环就会明显变慢。一个简单的优化是把交叉的for j 1:D循环向量化crossover_mask rand(1, D) CR; crossover_mask(jrand) true; U V .* crossover_mask X(i, :) .* (~crossover_mask);这样可以一次完成所有维度的交叉不需要内层for循环。另一个容易被忽视的性能瓶颈是数组动态增长。如果每次迭代都往一个变量后面追加结果Matlab会反复分配内存特别慢。所以代码里初始化了convergence zeros(1, MaxIt)每次迭代直接赋值这是很标准的做法。4.4 把代码抽成通用函数的复用技巧Sphere函数只是测试用例真正用DE时你可能要解的是自己的目标函数。如果每次换函数都复制整段脚本很容易改错。我习惯把DE主流程封装成一个函数只把适应度函数作为参数传进去。比如function [best_solution, best_fitness, convergence] de_solver(fitness_func, NP, D, F, CR, MaxIt, lb, ub)主循环里的fit fitness_func(X);和f_new fitness_func(U);只要改成调用函数句柄整套算法就能复用到几乎所有单目标连续优化问题。Sphere的适应度函数写成(x) sum(x.^2, 2)就可以了。这样一来你以后想测Rastrigin、Ackley只需要换函数句柄其他代码都不用动。我后来就是用这个套路快速搭了一个测试工具箱省了非常多的重复劳动。最后再分享一点个人经验DE虽然参数少但并不意味着可以无脑套默认参数。Sphere函数上F0.5、CR0.9很强换到多峰函数未必还是最优配置。我调试的时候习惯先固定NP用网格搜索快速扫一遍F和CR不需要太精细一个大方向的设置就能让性能差出好几个数量级。从这个角度讲Sphere更像是一面镜子它能照出你的DE实现是否标准但真正检验算法实力的永远是更难的问题。
阅读完成 · 觉得有帮助？
