很多刚接触智能优化算法的朋友都喜欢直接跑开源代码，却不深究算法内部到底做了什么。我自己最初做差分进化算法（Differential Evolution，DE）实验时，也是先把各种测试函数的标准结果表格跑出来，再回头补的数学推导。后来带新人时发现，与其直接丢一堆参考资料，不如用一段足够简单、又能完整展示DE工作流程的代码逐行拆开讲。这篇文章就基于单目标Sphere函数，手写一套结构清晰的Matlab源码，并把我调试过程中踩过的坑、调参数的经验一起写出来。内容不要求你有很强的编程功底，但至少得知道Matlab的矩阵、循环和函数文件是怎么组织的。

1. 选Sphere函数当试验场：不只是因为简单

1.1 Sphere函数长什么样

Sphere函数的标准数学形式是：

f(x) = Σ x_i^2，i = 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）包含以下步骤：

1. 初始化种群，每个个体是一个D维实数向量。
2. 计算所有个体的适应度值。
3. 对每个个体执行：
   • 从种群中随机选择三个与当前个体不同的个体；
   • 用其中两个的差分向量加权后加到第三个个体上，得到变异向量；
   • 将变异向量与当前个体按维度交叉，生成试验向量；
   • 试验向量越界时进行边界修复；
   • 比较试验向量与当前个体的适应度，保留更优者。
4. 重复步骤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); end

3.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 常用参数范围与实测经验

根据我自己的实验和平时参考的资料，给出一组比较常用的起点参数：

这里推荐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('F=0.3', 'F=0.5', 'F=0.9'); xlabel('迭代次数'); ylabel('最优适应度值'); grid on; hold off;

这个实验在10维Sphere上能很直观地看到：F=0.3时前期收敛慢但后期精度不差，F=0.9时前期长得快但可能出现震荡。多做几次这样的实验，比单纯记一套“最优参数”有用得多。

4.4 不要忽略随机种子

调参时一定要固定随机种子，否则两次实验之间种群的初始分布不同，收敛曲线会带有非常大的随机波动。我在自己实验脚本里习惯先加一句：

% 固定随机种子，使实验可复现 rng(0);

这样调参时才能真正确认是参数变了导致结果变化，而不是随机性在捣乱。这个习惯在我后面跑Rastrigin函数时帮助极大。

5. 跑起来之后的常见问题：不收敛、早熟、NaN的排查链路

5.1 收敛速度远低于预期

如果你按上面的代码跑10维Sphere，800代收敛到1e-20是很正常的。如果实际结果连1e-5都到不了，首先检查以下几项：

1. 边界范围是否正确。例如把下界和上界设反了，初始化种群生成在[100, -100]之外的区域，大部分个体直接越界，搜索几乎没有效果。
2. 交叉概率是否过小。CR = 0.3以下时，每代能更新的维度很少，收敛会明显变慢。
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，不要急着瞎猜，按这个顺序排查：

1. 边界中是否引入了NaN。例如lb或ub包含NaN，初始化时rand .* (ub - lb)就会产生NaN。
2. F是否为0。如果F = 0，变异向量直接退化为population(r1, :)，试验向量一旦被选择，大量个体变成同一个解，后续差分为0，种群僵死。这虽然不是NaN的来源，但会让适应度曲线一动不动。
3. 是否在自定义适应度函数里出现了除以0或者对数负值。如果将来把sphere_func改成其他测试函数，尤其要小心那些包含log(x)或1/x项的函数，参数越界时很容易产生Inf。

Sphere本身绝对不可能产生NaN，所以一旦出现，问题一定出在代码逻辑或参数设置上，不要怀疑函数本身。

5.4 画图观察种群演化

只记录最优适应度曲线是观察算法全貌的最基础方式。进阶一点的做法是画种群分布的散点图，看个体如何在搜索空间里移动。在二维Sphere函数上画这个效果特别明显：

% 修改主循环内部，在每次迭代后画出种群分布（仅限dim=2时） 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_{i=1}^{D-1} [100*(x_{i+1}-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”论文，一眼就能看穿作者到底改了哪一步、改动了哪里。调试算法就像拆解一台机器，只有亲手拼过一遍，才知道每个齿轮的咬合关系。希望这套代码能成为你读代码、写代码、调代码的第一步。