很多刚接触智能优化算法的朋友都喜欢直接跑开源代码却不深究算法内部到底做了什么。我自己最初做差分进化算法Differential EvolutionDE实验时也是先把各种测试函数的标准结果表格跑出来再回头补的数学推导。后来带新人时发现与其直接丢一堆参考资料不如用一段足够简单、又能完整展示DE工作流程的代码逐行拆开讲。这篇文章就基于单目标Sphere函数手写一套结构清晰的Matlab源码并把我调试过程中踩过的坑、调参数的经验一起写出来。内容不要求你有很强的编程功底但至少得知道Matlab的矩阵、循环和函数文件是怎么组织的。1. 选Sphere函数当试验场不只是因为简单1.1 Sphere函数长什么样Sphere函数的标准数学形式是f(x) Σ x_i^2i 1, 2, ..., D其中D是维度x_i是第i个决策变量。全局最优解是x* (0, 0, ..., 0)此时f(x*) 0。如果你在二维平面画它的等高线会得到一组同心圆中心在原点。三维空间里就是一组同心球面“Sphere”这个名字就是这么来的。很多刚接触智能算法的朋友会觉得这个函数太简单用来验证算法有点浪费。其实恰恰相反Sphere是测试优化算法收敛速度、稳定性和实现正确性的“基准中的基准”。因为我见过太多人在Rastrigin函数上跑出了看似漂亮的结果结果一检查发现代码里变量越界处理写错了导致部分维度根本没参与搜索。1.2 Sphere函数不可替代的作用Sphere函数具备几个很“友好”的性质特别适合作为验证代码正确性的第一步它是凸函数只有一个全局最优点没有局部最优干扰。算法收敛不到全局最优点时可以断定是算法本身的问题而不是函数地形太复杂。梯度方向永远指向原点最优点位于搜索空间的几何中心。这让我们可以直观地观察种群是否在向中心汇聚。最优值已知0不受数值误差影响。收敛精度可以直接用适应度值衡量非常方便。即使维度升高函数性质也不会突变。10维和100维的Sphere只是计算量线性增长没有本质变化。我从第一版代码开始就用Sphere做调试基准等确认DE实现没有逻辑错误之后才去跑其他更复杂的函数。这样做的好处是一旦后面出现非预期结果我能迅速判断是函数特性导致的还是算法实现的问题。1.3 单目标Sphere在各种实验中的定位在发表实验对比或者做课程作业时Sphere函数往往作为benchmark套件中的第一个测试函数。它给出的基线收敛图非常漂亮几乎是指数级下降用来直观展示你的算法是否能正确找到最优点。这篇博客里我也用它来展示DE的完整实现后面自己写算法或者给别人讲DE时这套代码都能直接在实验报告中复用。2. 差分进化算法的核心机制先理解再写代码2.1 从遗传算法到差分进化为什么需要差分策略遗传算法GA使用选择、交叉和变异三个算子核心在于编码和交叉操作的设计。差分进化算法属于进化算法的一个分支但它最大的不同在于变异操作DE用种群中个体的差分向量来产生扰动不需要单独的变异算子和复杂的编码方案直接在实数空间中操作。打个比方来说GA的变异更像是在当前解附近随意做一次小抖动而DE的变异是让一群个体告诉你“该往哪个方向走”再由方向向量放缩后形成新解。差分向量自带种群分布信息这在搜索前期能保留多样性在后期又能自动缩小步长收敛行为非常自洽。2.2 DE的标准流程分解一个经典的DE算法通常记为DE/rand/1/bin包含以下步骤初始化种群每个个体是一个D维实数向量。计算所有个体的适应度值。对每个个体执行从种群中随机选择三个与当前个体不同的个体用其中两个的差分向量加权后加到第三个个体上得到变异向量将变异向量与当前个体按维度交叉生成试验向量试验向量越界时进行边界修复比较试验向量与当前个体的适应度保留更优者。重复步骤3直到达到最大迭代次数或满足停止条件。其中每一步都不复杂但组合在一起之后整段代码的可读性和调试难度取决于你怎么组织数据结构。我一直习惯用二维矩阵存储整个种群行是种群个体列是决策变量维度。这样做的好处是与适应度向量分开管理迭代时逻辑清晰。2.3 三个算子背后的“为什么”变异算子Mutation。标准形式是v_i x_r1 F * (x_r2 - x_r3)其中F是缩放因子通常取0.4到0.9。差分向量(x_r2 - x_r3)提供了种群当前分布的方向信息。F决定了扰动的强度F过大会让搜索变成完全随机游走F过小则种群很快僵化。交叉算子Crossover。二项式交叉binomial crossover是最常见的选择对试验向量的每一个维度以概率CR决定是取变异向量的分量还是取当前个体的分量。但至少要保证一个维度来自变异向量否则试验向量就等于当前个体白白浪费一次评估。选择算子Selection。DE采用一对一贪婪选择试验向量不比当前个体差就替换。这里有个容易被忽略的细节比较适应度时是试验向量与当前个体比而不是与种群最优比。这种局部竞争机制让DE有种群内部的竞争压力又不至于让所有个体过早汇聚到同一个点。2.4 DE为何能兼顾探索与开发DE的强大之处在于差分向量会随着种群收敛而自动缩小。搜索初期种群散布范围大差分向量的模长也大相当于在做全局探索搜索后期所有个体集中到最优点附近差分向量的模长变得很小变异步长自动细化相当于在做局部精化。这个特性不需要额外设计衰减参数完全是自适应的。理解到这里代码怎么写就已经很明朗了。下面给出完整实现。3. Matlab源码逐段拆解一套可直接运行的DE实现3.1 函数文件与参数设置我习惯把DE的主逻辑封装成一个函数文件而不是写一个大型脚本。这样既方便通过命令行传入不同的Sphere维度、种群规模、迭代次数等参数做实验也方便后续扩展到其他测试函数。function [best_solution, best_fitness, fitness_history] de_sphere(pop_size, dim, F, CR, max_iter, lb, ub) % 差分进化算法DE/rand/1/bin求解Sphere函数 % 输入 % pop_size - 种群规模 NP % dim - 决策变量维度 D % F - 缩放因子 % CR - 交叉概率 % max_iter - 最大迭代次数 % lb, ub - 变量下界和上界可以是标量或维度相同的向量 % 输出 % best_solution - 最优解向量 % best_fitness - 最优适应度值 % fitness_history - 每代最优适应度记录 % 对输入边界做统一处理 if numel(lb) 1 lb lb * ones(1, dim); end if numel(ub) 1 ub ub * ones(1, dim); end3.2 种群初始化与适应度评估初始化阶段我使用均匀随机分布% 初始化种群 population lb rand(pop_size, dim) .* (ub - lb); fitness zeros(pop_size, 1); for i 1:pop_size fitness(i) sphere_func(population(i, :)); end % 记录初始最优 [best_fitness, best_idx] min(fitness); best_solution population(best_idx, :); fitness_history zeros(max_iter, 1); fitness_history(1) best_fitness;这里有一个很多人会忽视的点rand(pop_size, dim)生成的是0到1之间的均匀分布随机数用到边界缩放时注意点乘的顺序不要写反成lb rand .* (ub - lb)之外的奇怪表达式。边界向量化之后二维矩阵的每一列都会被正确映射到对应维度区间。3.3 主循环变异、交叉、选择循环体内的代码是DE的核心% 主循环 for g 1:max_iter for i 1:pop_size % 随机选择三个与当前个体不同的个体 r1,r2,r3 candidates 1:pop_size; candidates(i) []; idx candidates(randperm(pop_size - 1, 3)); r1 idx(1); r2 idx(2); r3 idx(3); % 变异v x_r1 F * (x_r2 - x_r3) mutant population(r1, :) F * (population(r2, :) - population(r3, :)); % 交叉二项式交叉至少一维来自变异向量 trial population(i, :); jrand randi(dim); for j 1:dim if rand CR || j jrand trial(j) mutant(j); end end % 越界处理采用随机重生策略 for j 1:dim if trial(j) lb(j) || trial(j) ub(j) trial(j) lb(j) rand * (ub(j) - lb(j)); end end % 选择贪婪保留更优个体 trial_fitness sphere_func(trial); if trial_fitness fitness(i) population(i, :) trial; fitness(i) trial_fitness; end end % 记录当前代全局最优 [best_fitness, best_idx] min(fitness); best_solution population(best_idx, :); fitness_history(g) best_fitness; end end这段代码我做了几个刻意的设计选择值得展开说明随机选取r1、r2、r3的方式。先用candidates 1:pop_size删除当前索引i再从剩余个体中随机抽取3个。这样能绝对保证变异个体不与当前个体重合。如果不做这一步DE的变异算子里就可能出现“自己和自己差分”的情况种群更新会受到错误引导收敛速度会明显受到干扰。交叉中jrand的作用。二项式交叉为了保证试验向量与当前个体不完全相同至少要有一个维度来自变异向量。这个维度由jrand randi(dim)随机指定。有初学朋友试过删掉jrand判断只保留rand CR结果算法在很多函数上依旧能跑但收敛稳定性会下降原因就是偶尔会出现一代完全不更新的情况。边界处理用了随机重生而不是简单的截断。截断处理最简单但一个明显副作用是多量个体贴到边界上尤其在高维问题里边界处的解通常不是最优区域。随机重生让越界个体回到搜索空间内部能保留更好的种群多样性。3.4 Sphere函数单独封装单独建立一个sphere_func.m文件function val sphere_func(x) % Sphere测试函数f(x) sum(x.^2) val sum(x.^2); end虽然当前最简单实现一行就够但单独封装有实实在在的好处后续切换到Rosenbrock或者Ackley函数时只需要改这一个文件DE主循环一行都不用动。这在做benchmark对比实验时非常实用。3.5 一个标准的调用示例假设我们要在10维空间、变量范围[-100, 100]内搜索种群60个个体迭代800次clear; clc; pop_size 60; dim 10; F 0.5; CR 0.9; max_iter 800; lb -100; ub 100; [best_solution, best_fitness, fitness_history] ... de_sphere(pop_size, dim, F, CR, max_iter, lb, ub); % 绘制收敛曲线 figure; semilogy(1:max_iter, fitness_history, LineWidth, 1.5); xlabel(迭代次数); ylabel(最优适应度值对数坐标); title(DE求解Sphere函数收敛曲线); grid on; fprintf(最优解\n); disp(best_solution); fprintf(最优适应度值%.12e\n, best_fitness);我用的是semoilogy而不是plot因为Sphere函数上DE收敛非常快200代之后最优适应度可能已经降到1e-10以下线性坐标下曲线会直接贴到横轴看不出任何细节。对数坐标能看到收敛曲线中“快速下降”和“缓慢精修”两个阶段这对判断算法状态至关重要。4. 参数到底怎么调种群规模、F、CR的搭配逻辑4.1 三个关键参数对搜索行为的影响DE的参数少是它流行的原因但少不代表不重要。三个主要参数分别是种群规模NP、缩放因子F、交叉概率CR每个参数对搜索行为的影响都不一样。种群规模NP。决定了差分向量的多样性来源。NP太小差分向量所能描述的搜索方向太少容易过早陷入局部最优NP太大每代的计算开销成倍增长。对于Sphere这类单峰函数10维情况下NP取40到80已经足够我到目前没见过有人用1000的大种群去解10维Sphere——纯属浪费算力。缩放因子F。控制差分扰动的步长。F过大变异向量会频繁跳出边界大量个体靠随机重生维持搜索行为退化成随机搜索F过小差分向量带来的更新步长太短种群可能在某个区域附近磨蹭很久收敛速度令人难以忍受。交叉概率CR。控制试验向量从变异向量继承多少比例的分量。CR越接近1试验向量越接近纯变异向量搜索的自由度更高CR越接近0试验向量越接近当前个体全局搜索能力下降但局部精细化能力上升。4.2 常用参数范围与实测经验根据我自己的实验和平时参考的资料给出一组比较常用的起点参数参数常用范围针对Sphere函数的推荐值说明NP30 ~ 20060维度低时取小值维度高时取大值F0.4 ~ 0.90.5无特殊需求时取0.5最稳妥CR0.3 ~ 0.950.9单峰连续函数建议取偏大的CRmax_iter300 ~ 200080010维Sphere收敛非常快不必取太大这里推荐CR 0.9是因为Sphere函数的各维度是独立的单维上不需要太多保留原信息的操作。相反在Rosenbrock这种维度间强耦合的函数上CR取值太大会导致搜索效率下降我一般会从0.5开始测试。4.3 一个快速参数实验观察法如果你想自己观察不同参数组合的差异可以写一个参数扫描循环把每个参数组合下的收敛曲线画在同一张图里figure; hold on; colors lines(3); test_F [0.3, 0.5, 0.9]; for k 1:3 [~, ~, history] de_sphere(60, 10, test_F(k), 0.9, 800, -100, 100); semilogy(1:800, history, Color, colors(k, :), LineWidth, 1.5); end legend(F0.3, F0.5, F0.9); xlabel(迭代次数); ylabel(最优适应度值); grid on; hold off;这个实验在10维Sphere上能很直观地看到F0.3时前期收敛慢但后期精度不差F0.9时前期长得快但可能出现震荡。多做几次这样的实验比单纯记一套“最优参数”有用得多。4.4 不要忽略随机种子调参时一定要固定随机种子否则两次实验之间种群的初始分布不同收敛曲线会带有非常大的随机波动。我在自己实验脚本里习惯先加一句% 固定随机种子使实验可复现 rng(0);这样调参时才能真正确认是参数变了导致结果变化而不是随机性在捣乱。这个习惯在我后面跑Rastrigin函数时帮助极大。5. 跑起来之后的常见问题不收敛、早熟、NaN的排查链路5.1 收敛速度远低于预期如果你按上面的代码跑10维Sphere800代收敛到1e-20是很正常的。如果实际结果连1e-5都到不了首先检查以下几项边界范围是否正确。例如把下界和上界设反了初始化种群生成在[100, -100]之外的区域大部分个体直接越界搜索几乎没有效果。交叉概率是否过小。CR 0.3以下时每代能更新的维度很少收敛会明显变慢。变异向量的差分方向是否被正确计算。建议在代码里加一行调试输出手动验证三代内变异向量是否符合r1 F * (r2 - r3)的公式。5.2 种群出现“早熟”如何判断与处理Sphere函数没有局部最优一般来说不会出现早熟。但如果你用同样的代码去跑Rastrigin函数就会遇到种群所有个体聚集到一起、差分向量长度趋近于0、种群完全停止移动的情况。判断早熟的标准不是“适应度不再下降”而是“种群个体之间的差异几乎为0”。我一般会计算每代种群的平均成对距离用来监测种群多样性。代码里可以这样加一段% 计算种群平均成对距离用于监测多样性 D pdist(population, euclidean); avg_distance mean(D);如果avg_distance在几十代之后掉到接近0基本上可以确认种群已经失去多样性了。应对早熟的手段很多比如增大F值、增加NP、或者在变异算子中引入当前最优个体的引导项变成DE/rand-to-best策略但那是后面进阶话题这里先不展开。5.3 出现NaN或Inf时的排查思路代码里如果出现NaN和Inf不要急着瞎猜按这个顺序排查边界中是否引入了NaN。例如lb或ub包含NaN初始化时rand .* (ub - lb)就会产生NaN。F是否为0。如果F 0变异向量直接退化为population(r1, :)试验向量一旦被选择大量个体变成同一个解后续差分为0种群僵死。这虽然不是NaN的来源但会让适应度曲线一动不动。是否在自定义适应度函数里出现了除以0或者对数负值。如果将来把sphere_func改成其他测试函数尤其要小心那些包含log(x)或1/x项的函数参数越界时很容易产生Inf。Sphere本身绝对不可能产生NaN所以一旦出现问题一定出在代码逻辑或参数设置上不要怀疑函数本身。5.4 画图观察种群演化只记录最优适应度曲线是观察算法全貌的最基础方式。进阶一点的做法是画种群分布的散点图看个体如何在搜索空间里移动。在二维Sphere函数上画这个效果特别明显% 修改主循环内部在每次迭代后画出种群分布仅限dim2时 figure; for g 1:max_iter plot(population(:, 1), population(:, 2), o); xlim([lb, ub]); ylim([lb, ub]); title([迭代次数: , num2str(g)]); drawnow; end你能清楚看到种群从大范围随机散布逐渐收拢到原点附近的过程。这个“收拢”的过程是否平滑、是否有个体反复被边界弹回都是调试时的宝贵信息。6. 从Sphere走向更复杂的测试函数代码扩展的两种思路6.1 直接替换适应度函数因为前面把sphere_func单独封装了扩展到其他测试函数只需要新建一个函数文件并修改主函数调用即可。比如换到Rosenbrock函数在rosenbrock_func.m中写function val rosenbrock_func(x) % Rosenbrock函数f(x) sum_{i1}^{D-1} [100*(x_{i1}-x_i^2)^2 (x_i-1)^2] val sum(100 * (x(2:end) - x(1:end-1).^2).^2 (x(1:end-1) - 1).^2); end注意Rosenbrock的全局最优点在x* (1, 1, ..., 1)这就给代码提了个醒不同测试函数有不同可行域和最优点初始化边界要选得合适才行。我的做法是写一个小脚本统一测试多个函数每个函数注册自己的上下界和最优值。6.2 添加约束惩罚函数法真实工程优化问题几乎都是有约束的。最简单的方式是把约束写成惩罚项加到适应度上function val constrained_func(x) % 原目标函数 f sum(x.^2); % 约束例如 x(1) x(2) 1 g 1 - (x(1) x(2)); penalty 1e6 * max(0, g)^2; val f penalty; end惩罚系数设多大来平衡约束违反量和目标函数值是一门学问。这里不展开但至少提供了从Sphere这样无约束测试函数迈向工程问题的自然路径——先把DE主逻辑跑通后面加惩罚项只是换一个适应度函数而已。6.3 高维问题与记录策略维度升高之后运行时间会明显增加。10维Sphere几百代瞬间跑完100维可能就要等上一会儿了。这之后有一个实用技巧早停策略。当最优适应度连续几十代都没有显著改善时提前终止迭代省时间十分明显。% 早停示例连续50代最优适应度下降不足1e-12则终止 patience 50; improve_threshold 1e-12; wait 0; for g 1:max_iter % ... 主循环 ... if abs(fitness_history(g) - fitness_history(max(1, g-1))) improve_threshold wait wait 1; if wait patience disp([迭代, num2str(g), 代后早停]); break; end else wait 0; end end注意早停阈值不要设太苛刻否则可能错过后期大幅精修的机会。在Sphere上误差1e-12以下时早停几乎无风险但在更复杂的函数上可以放宽到1e-6。6.4 多轮独立运行统计正式实验不能只看一次运行结果。我一般会把算法跑N次比如20次记录每次的最优值、收敛代数然后输出均值、标准差和中位数。这样才能说明算法性能稳定的说法是有数学依据的。Matlab脚本中用for r 1:repeats包一层即可不必修改DE主函数本身。最后再分享一点个人体会这套DE代码前前后后帮我带过好几批学生也作为模块嵌进过我自己的对比实验中。虽然现在很多工具箱里都有现成的ga或particleswarm函数但手写一遍DE的价值仍然不可替代——你会透彻理解每个算子的执行顺序、每个参数的作用以及边界处理选择对结果的影响。到后面你再去看那些所谓“改进的DE”论文一眼就能看穿作者到底改了哪一步、改动了哪里。调试算法就像拆解一台机器只有亲手拼过一遍才知道每个齿轮的咬合关系。希望这套代码能成为你读代码、写代码、调代码的第一步。