萤火虫算法优化实践：从理论到Matlab实现全解析

1. 萤火虫算法：不只是“追光”那么简单

大家好，我是老张，在智能优化算法这个领域摸爬滚打了十来年，用过的算法没有一百也有几十种。今天想和大家深入聊聊萤火虫算法。很多朋友第一次听说这个名字，可能会觉得它很浪漫，不就是模拟萤火虫发光吸引同伴嘛。但在我实际用它解决过电机设计、天线布局、甚至物流路径规划这些问题之后，我发现它的内核远比“追光”这个比喻要强大和精巧。

简单来说，萤火虫算法是一种模仿自然界萤火虫发光求偶行为的群智能优化算法。你可以把它想象成在一个漆黑的夜晚，有一群萤火虫在飞舞。每只萤火虫的亮度代表它找到的“食物”位置的好坏（在优化问题里，就是目标函数值的好坏，比如成本最低、效率最高）。亮度低的萤火虫会被亮度更高的萤火虫吸引，并向它靠近。通过这样一轮轮的“追逐”，整个群体最终会聚集到最亮的那几个点附近，也就是我们想找的最优解附近。

它特别适合谁呢？如果你正在用Matlab做课程设计、毕业设计，或者工作中遇到了那些参数多、关系复杂、用传统数学方法不好直接求解的优化问题，比如神经网络参数调优、工厂生产调度、图像处理中的特征选择，那么萤火虫算法很可能是一个让你眼前一亮的工具。它的优点很明显：概念直观、参数少、容易上手实现，而且对于多峰函数（有多个“山头”的最优解）的搜索能力很不错，不容易一下子就掉进局部最优的“坑”里出不来。

不过，我刚开始用它的时候也踩过坑。比如，所有萤火虫都互相吸引，计算量一下子就上去了，问题规模稍大程序就跑得慢；又比如，参数设置不合适，虫子们要么乱飞找不到北，要么过早地挤成一团，错过了更好的区域。这些实战中的细节，恰恰是教科书和论文里常常一笔带过，却又至关重要的部分。接下来，我就结合多年的使用经验，带你从最核心的原理掰开揉碎了讲，一直到一个可以直接在Matlab里跑起来的、带详细注释的完整代码。

2. 核心原理拆解：亮度、吸引力与移动的数学舞蹈

要真正玩转一个算法，死记硬背步骤是没用的，必须理解它每一步背后的“为什么”。萤火虫算法的核心，就围绕着三个关键概念展开：亮度、吸引力和移动。这三者构成了一场精妙的数学舞蹈。

2.1 亮度：你的“江湖地位”

在算法里，每只萤火虫的位置对应优化问题的一个候选解。比如我们要找函数 y = f(x) 的最小值，那么萤火虫的坐标 (x1, x2, ...) 就是自变量，而它的亮度 I就直接由目标函数值 f(x) 决定。对于最小化问题，函数值越小，解越好，我们反而让亮度越高。这是一个简单的映射关系，比如可以设定亮度 I = 1 / (f(x) + ε)，这样 f(x) 越小，I 就越大。

这里有个非常重要的细节：亮度随距离衰减。在自然界中，光越远越暗。算法中也模拟了这一点。萤火虫j在萤火虫i眼中看到的亮度，并不是j本身的绝对亮度，而是经过距离衰减后的亮度。公式通常这样表示：

I(r) = I0 * exp(-γ * r^2)

这里的I0是萤火虫自身的原始亮度（即由目标函数值决定的亮度），γ是一个叫做光吸收系数的重要参数，r是两只萤火虫之间的距离。exp(-γ * r^2)这一项就是衰减因子。γ越大，光强随距离衰减得越快，这意味着萤火虫的“视野”变短，只能看到附近的同伴；γ越小，衰减越慢，远处的萤火虫也能对其产生吸引。这个参数直接影响了算法的全局搜索与局部开发能力，我们后面调参时会重点讲。

2.2 吸引力：决定你跟谁走、走多近

看到了更亮的同伴，心动了，但到底有多想靠近它呢？这就是吸引力 β所决定的。吸引力同样是距离的函数，通常定义为：

β(r) = β0 * exp(-γ * r^2)

看，它和亮度衰减用了同一个衰减因子。β0是最大吸引力，也就是两只萤火虫距离为0时的吸引力（通常设为1）。这个公式意味着，吸引力也随距离增加而呈指数衰减。一只很亮的萤火虫，如果离你太远，对你的吸引力可能还不如一只亮度稍差但就在你身边的萤火虫。这种机制非常符合直觉，也保证了种群不会盲目地全部涌向当前全局最优解，有助于维持种群的多样性，避免过早收敛。

2.3 移动：位置更新的具体一步

这是算法最“动”起来的一步。对于任意一只萤火虫i，它会遍历种群中所有比它亮的萤火虫j（注意，这里通常是全互连比较，计算量就来源于此），计算每一只对自己的吸引力，然后综合地更新自己的位置。位置更新公式是核心：

xi_new = xi_old + β(r) * (xj - xi_old) + α * (rand - 0.5)

我们来拆解这个公式的三部分：

1. xi_old：这是萤火虫i当前的位置。
2. β(r) * (xj - xi_old)：这是确定性移动部分。它驱使i向着更亮的j移动。移动的步长由吸引力β(r)控制，距离越近，吸引力越大，这一步的“拉力”就越强。这体现了算法的“开发”能力，即在有希望的区域进行精细搜索。
3. α * (rand - 0.5)：这是随机扰动部分。α是步长因子，rand是一个在[0,1]内均匀分布的随机数，所以(rand - 0.5)就是一个在[-0.5, 0.5]之间的随机扰动。这一项至关重要！它为移动增加了随机性，使得萤火虫有机会跳出当前吸引它的那个点的“势力范围”，去探索更广阔的空间。这是算法“勘探”能力的主要来源。α的大小控制着探索的力度。

和GSO算法的一个小对比：有些朋友可能也听过另一种萤火虫优化算法（GSO）。它们最大的区别在于“社交网络”不同。FA（本文介绍的）是“全球社交”，每只虫子都和所有更亮的虫子比较。而GSO是“局部社交”，每只虫子只关注自己“决策半径”内的邻居。FA的全局信息交流更充分，但计算更慢；GSO计算快，但容易陷入局部最优。在实际应用中，FA因其结构简单、性能稳定，使用更为广泛。

3. 手把手Matlab实现：从零搭建你的第一个FA模型

理论说得再多，不如一行代码。下面我就带你用Matlab从头实现一个标准的萤火虫算法，用来求解一个经典的多峰测试函数——Ackley函数的最小值。这个函数全局最优点在(0,0,...,0)，值为0，但有很多局部极值点，很适合检验算法的性能。

3.1 准备工作与参数设置

首先，我们明确要解决的问题。Ackley函数在二维形式下长这样： f(x, y) = -20 * exp(-0.2 * sqrt(0.5*(x^2+y^2))) - exp(0.5*(cos(2πx)+cos(2πy))) + e + 20 我们的目标是在定义域内，比如 x, y 都在 [-5, 5] 之间，找到使 f 最小的 (x, y)。

打开Matlab，新建一个脚本。我们先定义算法运行所需的全部参数。这些参数就像厨师的调料，搭配好了才能做出好菜。

%% 萤火虫算法参数设置 clear all; close all; clc; % 问题定义 n = 2; % 问题维度 (Ackley函数是2维) lb = -5 * ones(1, n); % 变量下界 ub = 5 * ones(1, n); % 变量上界 % 算法参数 pop_size = 30; % 萤火虫种群大小 max_gen = 100; % 最大迭代次数 alpha = 0.2; % 步长因子 (通常0.1~0.5) beta0 = 1.0; % 最大吸引力 (通常设为1) gamma = 1.0; % 光吸收系数 (关键参数！通常0.1~10)

参数选择经验谈：

• pop_size：种群数量。一般20-50。问题越复杂、维度越高，种群可以适当增大，但计算量也会增加。30是个不错的起点。
• max_gen：迭代次数。根据问题难度调整，可以先设100-500次观察收敛情况。
• alpha：随机步长。这是平衡“探索”与“开发”的关键。我习惯在迭代初期设得大一点（如0.5），鼓励探索；后期逐渐减小（如降到0.01），进行精细开发。这里我们先用一个固定值。
• gamma：最关键的参数！它控制着吸引力和亮度随距离衰减的速度。gamma很小（如0.01）时，衰减慢，远处虫子也能相互吸引，全局搜索能力强，但收敛慢；gamma很大（如10）时，衰减快，只有很近的虫子相互影响，算法很快局部收敛，但可能陷入局部最优。通常从1附近开始调试。

3.2 初始化种群与计算初始亮度

接下来，我们在搜索空间内随机生成初始的萤火虫种群，并计算它们的初始亮度（即目标函数值）。

%% 初始化种群 % 在上下界范围内随机生成种群位置 pop = lb + (ub - lb) .* rand(pop_size, n); % 计算每个萤火虫的初始目标函数值（对于最小化问题，值越小越好） fitness = zeros(pop_size, 1); for i = 1:pop_size fitness(i) = ackley_func(pop(i, :)); % 调用Ackley函数 end % 对于最小化问题，我们将适应度值转换为亮度：值越小，亮度越高 % 这里采用简单直接的倒数关系（加个小常数防止除零） I = 1 ./ (fitness + 1e-10); % 记录历史最优解 [best_I, best_index] = max(I); % 注意，亮度I越大越好 best_solution = pop(best_index, :); best_fitness = fitness(best_index); history_best = zeros(max_gen, 1); % 记录每一代的最优值 history_best(1) = best_fitness;

这里我封装了一个ackley_func.m函数文件，内容如下：

function [y] = ackley_func(x) % 计算Ackley函数值 % 输入x是一个行向量 d = length(x); sum1 = sum(x.^2); sum2 = sum(cos(2 * pi * x)); y = -20 * exp(-0.2 * sqrt(sum1/d)) - exp(sum2/d) + 20 + exp(1); end

3.3 核心迭代循环：吸引、移动与更新

这是算法的引擎部分，对应原理中的移动公式。我们需要一个双循环：对于每只萤火虫i，检查所有其他萤火虫j，如果j比i亮，则i向j移动。

%% 开始迭代优化 for gen = 2:max_gen % 从第2代开始，第1代已经初始化了 % 复制当前种群，用于计算距离和更新 new_pop = pop; for i = 1:pop_size for j = 1:pop_size % 如果萤火虫j比i亮（对于最小化问题，适应度值更小，亮度更高） if fitness(j) < fitness(i) % 注意这里比较的是适应度值，越小越好 % 计算萤火虫i和j之间的欧氏距离 r_ij = norm(pop(i, :) - pop(j, :)); % 计算基于距离的吸引力 beta = beta0 * exp(-gamma * r_ij^2); % 更新萤火虫i的位置（向j移动） % 公式: new = old + 吸引力 * (xj - xi) + 随机扰动 random_vec = rand(1, n) - 0.5; % 生成[-0.5, 0.5]的随机向量 new_pop(i, :) = pop(i, :) + beta * (pop(j, :) - pop(i, :)) + alpha * random_vec; % 确保新位置在边界内（一种简单的边界处理方式：越界则拉回边界） new_pop(i, :) = max(new_pop(i, :), lb); new_pop(i, :) = min(new_pop(i, :), ub); % 计算新位置的适应度 new_fitness = ackley_func(new_pop(i, :)); % 如果新位置更好，则接受这次移动 if new_fitness < fitness(i) fitness(i) = new_fitness; pop(i, :) = new_pop(i, :); I(i) = 1 / (new_fitness + 1e-10); % 更新亮度 else % 如果移动后效果更差，则保持原位置（也可以有其他策略） new_pop(i, :) = pop(i, :); % 这里我们简单拒绝此次移动 end end end end % 每一代结束后，更新全局最优解 [current_best_fitness, idx] = min(fitness); if current_best_fitness < best_fitness best_fitness = current_best_fitness; best_solution = pop(idx, :); end history_best(gen) = best_fitness; % 每10代显示一次进度 if mod(gen, 10) == 0 fprintf('迭代次数: %d, 当前最优值: %.6f\n', gen, best_fitness); end end

代码细节与坑点提示：

1. 亮度比较：我们存储的是适应度值fitness（目标函数值，越小越好）。在比较时，直接判断if fitness(j) < fitness(i)，这意味着j的位置更好，i应该被j吸引。
2. 距离计算：norm()函数计算欧氏距离，这是最常用的方式。对于超高维问题，计算距离可能成为瓶颈。
3. 边界处理：更新位置后，一定要检查是否超出了搜索空间的边界。我用了最简单的“反射”或“夹紧”方式，直接将其设置为边界值。更复杂的方法可以是用随机值替代，或者进行反射。
4. 移动接受准则：上述代码采用了“贪婪接受”，只有移动后解变好才更新位置。这是标准做法。你也可以尝试总是移动，然后在整个一代结束后统一评估，但前者更简单高效。
5. 随机扰动：alpha * (rand - 0.5)这里的rand是生成一个标量，然后与向量(xj - xi)相乘？不对！注意，随机扰动应该是一个与位置维度相同的向量。所以代码中我用了rand(1, n) - 0.5来生成n维随机扰动向量。这是一个新手常犯的错误，会导致搜索方向偏颇。

3.4 结果可视化与性能分析

迭代完成后，我们当然要看看效果如何。

%% 结果可视化 fprintf('\n优化结束！\n'); fprintf('找到的最优解: [%.6f, %.6f]\n', best_solution(1), best_solution(2)); fprintf('对应的最优函数值: %.10f\n', best_fitness); % 绘制收敛曲线 figure(1); plot(1:max_gen, history_best, 'b-', 'LineWidth', 1.5); xlabel('迭代次数'); ylabel('最优适应度值'); title('萤火虫算法收敛曲线'); grid on; % 绘制搜索过程（二维问题可以画） if n == 2 figure(2); % 绘制Ackley函数曲面 [X, Y] = meshgrid(linspace(lb(1), ub(1), 100), linspace(lb(2), ub(2), 100)); Z = zeros(size(X)); for i = 1:size(X,1) for j = 1:size(X,2) Z(i,j) = ackley_func([X(i,j), Y(i,j)]); end end surf(X, Y, Z, 'EdgeColor', 'none', 'FaceAlpha', 0.6); hold on; % 绘制最终种群分布 scatter3(pop(:,1), pop(:,2), fitness, 40, 'r', 'filled'); plot3(best_solution(1), best_solution(2), best_fitness, 'yp', 'MarkerSize', 20, 'MarkerFaceColor', 'y'); xlabel('x1'); ylabel('x2'); zlabel('f(x)'); title('最终萤火虫种群在解空间中的分布'); hold off; colorbar; end

运行这段代码，你会在命令行看到迭代过程中最优值的下降情况，并得到两张图：一张是收敛曲线，可以看到算法是如何一步步逼近最优值的；另一张是种群分布图（针对二维问题），可以看到最终的萤火虫们是否聚集在了全局最优点 (0,0) 附近。第一次跑，你可能发现结果不一定完美，这很正常，因为算法中有随机性，而且参数是初始设置。这恰恰引出了我们下一个最重要的话题：调参优化。

4. 性能调优与实战技巧：让FA在你的问题上发光

如果你直接运行上面的代码，可能会发现有时候效果很好，能非常接近0，有时候却停滞在某个局部最优。这说明算法的性能对参数和问题本身很敏感。下面分享我积累的一些调优经验和实战技巧。

4.1 关键参数的影响与动态调整策略

1. 光吸收系数 γ：这是最重要的参数。你可以把它理解为萤火虫的“视力”范围。

   • γ 太小（如0.01）：衰减很慢，吸引力几乎不随距离减弱。这意味着远处的亮萤火虫也能对当前萤火虫产生很强的吸引力。好处是全局搜索能力极强，种群信息交流充分，不容易陷入局部最优。坏处是收敛速度会非常慢，因为所有个体都受到全局最优个体的强烈牵引，可能一直在震荡。
   • γ 太大（如10）：衰减极快，吸引力随距离急剧下降。萤火虫基本上只能看到和吸引紧挨着自己的邻居。这导致算法局部搜索能力强，收敛速度快，但极易陷入局部最优，因为种群被分割成多个小团体，缺乏全局信息交流。
   • 调参建议：没有一个放之四海而皆准的值。对于搜索范围大的问题（如我们的[-5,5]），可以从1开始尝试。一个有效的策略是动态调整γ：在迭代初期，设置较小的γ值，增强全局探索；随着迭代进行，逐渐增大γ值，加强局部开发，加快收敛。例如：gamma = 0.1 + (1.0 - 0.1) * (gen/max_gen);

2. 步长因子 α：这是探索的“发动机”。

   • 固定 α：像我们代码中那样，简单但效果有限。前期可能探索不足，后期可能扰动过大破坏好解。
   • 动态衰减 α：这是最常用且有效的策略。在迭代初期，给一个较大的α（如0.5），让萤火虫有足够大的随机步长去探索未知区域；随着迭代进行，逐渐减小α（如线性衰减到0.01），让算法后期专注于在最有希望的区域进行精细搜索。实现很简单：alpha = alpha_init * (1 - gen/max_gen)^delta，其中delta是衰减系数，通常取1或2。

3. 种群大小 pop_size 和最大吸引力 β0：

   • pop_size：越大，搜索能力越强，但每次迭代的计算量也越大（O(N^2)复杂度）。对于简单问题，20-30足够；对于复杂、多峰问题，可以增加到50-100。需要在效果和效率间权衡。
   • beta0：通常设为1即可。它代表了吸引力的基准强度。有些改进算法会动态调整β0，但标准FA中固定为1效果就不错。

4.2 改进思路：提升效率与精度

标准FA有两个明显的缺点：计算复杂度高（O(N^2)）和后期收敛慢。这里提供几个我实践过的改进思路，你可以尝试融入你的代码：

1. 邻域搜索结构（降低复杂度）：不让每只萤火虫都和所有其他虫子比较，而是只和一定“拓扑结构”内的邻居比较。比如：

   • 环形拓扑：每只虫子只和左右各k只虫子比较。
   • 随机邻域：每代随机为每只虫子选择m个邻居。
   • 精英引导：只让普通萤火虫被当前最优的少数几个“精英”萤火虫吸引。 这样可以大幅减少距离计算和比较的次数，尤其在大规模种群时效果显著。

2. 自适应参数与混合策略：

   • 结合其他算法的优点：比如，在FA的更新公式中引入粒子群算法（PSO）的“惯性”概念，让萤火虫的移动不仅受吸引力和随机扰动影响，还保留一部分上一代的速度。
   • 莱维飞行（Lévy Flight）：用莱维飞行分布代替简单的均匀随机扰动alpha * (rand-0.5)。莱维飞行具有长距离跳跃的特点，能显著增强算法的全局探索能力，避免陷入局部最优。Matlab中可以通过Mantegna算法生成莱维飞行步长。

3. 针对具体问题的编码与约束处理：

   • 离散优化问题：FA最初是为连续优化设计的。对于旅行商问题（TSP）、调度问题等离散问题，需要重新设计位置表示和移动算子。例如，位置可以表示为城市的排列，移动操作可以定义为两种排列之间的转换（如交换、插入、逆序）。
   • 约束处理：对于有约束的优化问题（比如x+y<10），需要在更新位置后处理约束。常用方法有罚函数法（将违反约束的程度加到目标函数值上，使其变“差”）、可行解保持法（只向可行域内的萤火虫移动）等。

4.3 一个简单的动态参数示例

让我们修改一下之前的代码，加入动态的α和γ，看看效果：

% 在迭代循环内部（gen循环内）替换掉固定的alpha和gamma alpha_current = alpha_init * (1 - gen/max_gen); % 线性衰减 gamma_current = gamma_init + (gamma_final - gamma_init) * (gen/max_gen); % 线性增加 % 然后在移动更新公式中使用这两个动态参数 beta = beta0 * exp(-gamma_current * r_ij^2); random_vec = rand(1, n) - 0.5; new_pop(i, :) = pop(i, :) + beta * (pop(j, :) - pop(i, :)) + alpha_current * random_vec;

你可以设置alpha_init=0.5,alpha_final=0.01,gamma_init=0.1,gamma_final=1.0。多跑几次对比一下，你会发现收敛曲线通常会更平滑、更稳定，找到全局最优的概率也更高。

最后，我想说，算法学习没有捷径。理解原理是第一步，动手实现是第二步，而根据具体问题调整、改进甚至创新，才是真正掌握它的标志。萤火虫算法就像一个灵活的工具箱，基础版本已经能解决很多问题，而它的各种变体和改进思路则为你应对更复杂的挑战提供了可能。我建议你先把标准版本的代码吃透，跑通，然后尝试用我上面提到的改进方法去修改代码，观察性能变化。这个过程里遇到的每一个报错和每一次性能提升，都是最宝贵的经验。