阳光生长优化算法PGA：原理、Matlab实现与参数调优指南

把“植物怎么晒太阳”变成一套优化算法，这个点子我第一次看到时觉得挺新鲜，后来动手复现了Matlab版本，发现它不仅思路有趣，代码结构也适合拿来当学习智能优化算法的入门模板。今天想聊的就是这个算法——阳光生长优化算法（Polychromatic Glow Optimization Algorithm，简称PGA），以及它对应的Matlab实现细节。

这个算法的中文名很形象，英文全称里Polychromatic Glow指的是“多色辉光”，放到算法语境里其实就是模拟植物在多色光源照射下完成光合作用、向光生长和群落光竞争的过程。和粒子群、遗传算法这类大家耳熟能详的元启发式方法相比，PGA的思路更贴近生态学：用不同波长的光代表不同的搜索策略，用植物茎秆的伸长方向代表解位置的更新方向，用光竞争机制来处理种群淘汰和多样性维护。文章后面我会把核心原理、数学公式、完整的Matlab代码框架、参数调优心得和典型问题全部拆开讲，适合刚接触智能优化算法的读者，也适合想找一个容易改造的算法框架用来做工程优化的开发者。

1. 算法背景与灵感来源

1.1 植物光照策略与优化问题的映射关系

先说算法灵感。植物生长离不开光，但不同波长的光对植物的作用千差万别：红光能促进茎的伸长和光合产物积累，蓝光能刺激叶绿素合成和细胞分化，而紫光或紫外波段的光往往参与调节植物的次生代谢和抗逆响应。植物在光照环境下并不是被动接受阳光，而是会主动调整姿态——向光性生长就是典型的例子，哪里光强，茎秆就往哪里长，同时还会通过改变叶片的朝向去争夺更好的光照位置。

把这个现象映射到优化问题里，思路就打开了。解空间可以看成一片“森林”，每个候选解就是一株“植物”，适应度值就是植物吸收光能后积累的“生物量”。全局最优解是光强最大的位置，其他个体要不断向这个光源靠拢。可是如果所有个体都只顾着冲向当前最优，种群很快会挤在一起，多样性丧失后算法就容易早熟。于是自然界的“多光谱”策略就派上了用场：让不同状态的个体使用不同的光质策略，有的负责大步探索新区域，有的负责在小范围内精修，有的负责跳出去寻找新的光照缝隙。

这样设计的好处在于，它不是靠单一公式从头到尾硬跑，而是把搜索过程拆成了拥有明确任务分工的协作模式。红光策略承担探索任务，对应算法里的大范围游走；蓝光策略承担开发任务，对应向最优解收缩；紫光策略承担精细调节任务，对应局部微调。三者比例还会随种群状态动态变化，避开“一刀切”带来的早熟或停滞。

1.2 PGA算法总体流程设计

整个PGA算法的主流程并不复杂，可以用下面几条写清楚：

1. 随机撒播初始种群，把每个个体当作一粒种子。
2. 计算每个个体的适应度，记录当前全局最优光源位置。
3. 根据某种策略，为每个个体分配一种光质（红光、蓝光或紫光）。
4. 个体按照对应光质的位置更新公式移动，形成一次“生长”。
5. 对越界个体做边界约束处理，计算新适应度，按接受准则决定是否保留。
6. 每隔固定代数执行一次光竞争淘汰，让弱势个体在最优光源附近重新播种。
7. 迭代至最大代数，输出全局最优解和收敛曲线。

从工程角度看，这个流程和大多数群体优化算法兼容性很好，不需要额外引入复杂的数据结构，Matlab实现时也比较顺手。核心工作量集中在“光质分配”和“位置更新”这两个环节上，只要把这两个模块写清楚了，整个算法的骨架就已经成型。

2. 核心数学原理

2.1 光通量模型与光照强度计算

在PGA里，个体感知到的光照强度不是恒定不变的，而是随自身位置与光源之间的距离变化。我采用了一个带衰减系数的指数模型来模拟这种衰减：

[ I_i(t) = I_0 \cdot e^{-\alpha \cdot d_i(t)} \cdot \rho(t) ]

其中：

• ( I_0 ) 是基准光强，通常取1。
• ( \alpha ) 是光吸收系数，控制光强随距离衰减的速度。
• ( d_i(t) = |X_i(t) - X_{best}(t)|_2 ) 是个体与当前最优解之间的欧几里得距离。
• ( \rho(t) ) 是随迭代进程变化的全局光照系数，用来模拟季节或太阳高度角的变化。

这个模型的直观含义是：距离最优解越远的个体，感知到的光强越弱，它们的状态越差，越有必要进行大范围探索；距离越近的个体，光强越强，适合在局部做精细化搜索。( \rho(t) ) 我一般设计成从1缓慢衰减到0的余弦函数，前期保证较强的探索强度，后期逐步收敛到精细开发。

2.2 三种光质下的位置更新公式

三种光质对应三类不同的位置更新机制。

红光策略模拟植物在遮蔽环境下“徒长”的行为，目标是跳出当前位置，向远处探索。更新公式为：

[ X_i(t+1) = X_{best}(t) + L \cdot (X_{best}(t) - X_i(t)) ]

这里 ( L ) 是Levy飞行随机向量。Levy飞行的特点是步长分布存在重尾，偶尔会出现大步长跳跃，非常适合在搜索空间中制造“变异”，避免种群陷入局部极值。

蓝光策略模拟正常光合作用下的趋光生长，个体朝最优光源方向收缩，同时引入种群中其他个体的差异信息来保持多样性：

[ X_i(t+1) = X_i(t) + r_1 \cdot (X_{best}(t) - X_i(t)) + r_2 \cdot \beta \cdot I_i(t) \cdot (X_{p}(t) - X_{q}(t)) ]

其中 ( r_1, r_2 ) 是0到1之间的随机向量，( X_p(t) ) 和 ( X_q(t) ) 是从种群中随机选择的两个不同个体，( \beta ) 是生长系数。这一项借鉴了差分进化中“差异向量”的思想，让个体既能朝最优靠拢，又能从种群内部的差异中获得扰动。

紫光策略模拟植物在富光环境下的精细分化，只在小邻域内做微小的高斯扰动：

[ X_i(t+1) = X_i(t) + \sigma \cdot N(0,1) ]

其中 ( \sigma ) 是扰动标准差，一般设置为搜索区间宽度的0.05到0.1倍。这样可以让适应度较好的个体在已有位置附近做精修，提升收敛精度。

2.3 光竞争淘汰机制的数学表达

光竞争的生态意义是：高大植物会遮挡矮小植物的光照，被严重遮阴的植物生长受阻甚至死亡，但死亡留下的空隙又会让新种子有机会萌发。在PGA中，我每隔固定周期淘汰适应度最差的一部分个体，然后在当前最优光源附近通过高斯扰动重新生成新个体：

[ X_{new} = X_{best}(t) + \sigma_w \cdot N(0,1) ]

其中 ( \sigma_w ) 是重新播种的扰动宽度，我通常设置为搜索区间宽度乘0.2。这个机制的作用是双重的：一方面清退无用个体，减少无效计算；另一方面在最优解附近补充新样本，增强局部开发能力。淘汰比例一般取20%左右，周期取5到10代比较合适。

3. Matlab代码实现与关键参数解析

3.1 主函数框架结构

先给出PGA主函数的完整骨架。这个函数接口风格和经典粒子群算法的封装方式相似，方便替换成自己的目标函数。

function [Best_pos, Best_fitness, Convergence_curve] = PGA(SearchAgents, Max_iter, lb, ub, dim, fobj) % PGA: 阳光生长优化算法 % 输入: % SearchAgents - 种群规模 % Max_iter - 最大迭代次数 % lb, ub - 搜索下界和上界 % dim - 决策变量维度 % fobj - 目标函数句柄 % 输出: % Best_pos - 全局最优位置 % Best_fitness - 全局最优适应度 % Convergence_curve - 收敛曲线 % 基础参数 alpha = 0.5; % 光吸收系数 I0 = 1.0; % 基准光强 beta_growth = 1.5; % 生长系数 sigma_ratio = 0.1; % 紫光扰动宽度占区间宽度比例 sigma_restart = 0.2; % 重播种扰动宽度占区间宽度比例 elimination_period = 10;% 光竞争淘汰周期 elimination_rate = 0.2; % 淘汰比例 % 初始化种群: 种子随机撒播 Positions = lb + rand(SearchAgents, dim) .* (ub - lb); Fitness = zeros(SearchAgents, 1); for i = 1:SearchAgents Fitness(i) = fobj(Positions(i, :)); end [Best_fitness, best_idx] = min(Fitness); Best_pos = Positions(best_idx, :); Convergence_curve = zeros(1, Max_iter); for t = 1:Max_iter % 全局光照系数: 余弦衰减 rho = cos(pi / 2 * t / Max_iter); if rho < 0.05 rho = 0.05; end % 按适应度排序, 确定每个个体的光质分配 [~, sorted_idx] = sort(Fitness); rank_map = zeros(1, SearchAgents); for r = 1:SearchAgents rank_map(sorted_idx(r)) = r; end for i = 1:SearchAgents old_position = Positions(i, :); old_fitness = Fitness(i); % 计算个体感知到的光通量 distance = norm(Positions(i, :) - Best_pos); light_intensity = I0 * exp(-alpha * distance) * rho + 0.01; % 根据排名分配光质: 前30%紫光, 中间40%蓝光, 后30%红光 rank_i = rank_map(i); if rank_i <= round(SearchAgents * 0.3) % 紫光策略: 局部精修 sigma = sigma_ratio * (ub - lb); Positions(i, :) = Positions(i, :) + randn(1, dim) .* sigma; elseif rank_i <= round(SearchAgents * 0.7) % 蓝光策略: 趋光生长 + 种群差异扰动 r1 = rand(1, dim); r2 = rand(1, dim); idx1 = randi(SearchAgents); idx2 = randi(SearchAgents); while idx2 == idx1 idx2 = randi(SearchAgents); end Positions(i, :) = Positions(i, :) + ... r1 .* (Best_pos - Positions(i, :)) + ... r2 .* beta_growth .* light_intensity .* ... (Positions(idx1, :) - Positions(idx2, :)); else % 红光策略: Levy飞行探索 levy_step = levy_flight(dim); Positions(i, :) = Best_pos + levy_step .* (Best_pos - Positions(i, :)); end % 边界约束 Positions(i, :) = max(Positions(i, :), lb); Positions(i, :) = min(Positions(i, :), ub); % 评估新位置 new_fitness = fobj(Positions(i, :)); % 接受准则: 更优则接受, 更差则按概率接受 if new_fitness < old_fitness Fitness(i) = new_fitness; else scale = abs(old_fitness) + 1; T = 1 - t / Max_iter; accept_prob = exp(-(new_fitness - old_fitness) / (scale * T + eps)); if rand < accept_prob Fitness(i) = new_fitness; else Positions(i, :) = old_position; Fitness(i) = old_fitness; end end end % 更新全局最优 [current_best, best_idx] = min(Fitness); if current_best < Best_fitness Best_fitness = current_best; Best_pos = Positions(best_idx, :); end % 光竞争淘汰机制 if mod(t, elimination_period) == 0 && t < Max_iter [~, sort_idx] = sort(Fitness); num_eliminate = max(1, floor(SearchAgents * elimination_rate)); sigma_w = sigma_restart * (ub - lb); for j = 1:num_eliminate idx = sort_idx(end - j + 1); Positions(idx, :) = Best_pos + randn(1, dim) .* sigma_w; Positions(idx, :) = max(Positions(idx, :), lb); Positions(idx, :) = min(Positions(idx, :), ub); Fitness(idx) = fobj(Positions(idx, :)); end end Convergence_curve(t) = Best_fitness; end end

这段代码里我做了三处容易被忽略但又很重要的设计。第一个是rank_map的构建方式，我用一次排序就能完整得到每个个体对应的档次，避免了每次循环里反复计算相对适应度导致的除零隐患。第二个是接受准则里的scale归一化，直接对适应度差值做概率计算时，遇到数值量级很大的目标函数很容易让概率失真，除以一个与当前适应度同量级的scale后稳定得多。第三个是淘汰周期避开最后一轮，防止在算法即将结束时突然重播种，破坏已经收敛的解。

3.2 光质分配策略的实现逻辑

光质分配是整个算法的灵魂。我采用的是“精英精修、劣势探索”的生态位分化逻辑，而不是简单的随机分配。这样做的好处是可以让适应度靠前的个体稳定地对当前最优区域做局部挖掘，同时让适应度靠后的个体持续向外探索，维持种群的全局覆盖能力。

代码里用前30%个体走紫光、中间40%走蓝光、后30%走红光的分段方式。实际调试时，这三个比例是可以调整的。如果目标函数多峰性强，需要更强的探索能力，可以把红光的比例提高到40%，紫光压缩到20%。相反，如果目标函数相对平坦、局部极值少，可以加大紫光比例，让算法更快收敛到高精度区域。

排名分配还有一个隐性的好处：它天然给种群引入了“自适应”特性。随着迭代推进，种群整体适应度不断提升，同样一个个体可能在前期排名靠后，只能走红光探索；后期排名进入前30%，自动切换到紫光精修。这种角色转换不需要额外判断，完全由种群自身状态驱动，实现起来非常简洁。

3.3 Levy飞行生成函数

Levy飞行在红光策略中承担的是重尾随机游走任务。Matlab里生成Levy随机数有现成的方法，下面是常用的实现：

function L = levy_flight(dim) % 生成Levy飞行随机步长向量 beta = 1.5; sigma = (gamma(1 + beta) * sin(pi * beta / 2) / ... (gamma((1 + beta) / 2) * beta * 2^((beta - 1) / 2)))^(1 / beta); u = randn(1, dim) * sigma; v = randn(1, dim); L = u ./ abs(v).^(1 / beta); % 防止步长过大造成越界震荡 L = L ./ max(abs(L)); end

这里有一个常见的坑：直接用原始Levy步长时，可能产生幅值极大的值，导致个体一下子飞出搜索边界很远。虽然边界约束会把它拉回来，但频繁越界会浪费大量评估次数。我在代码末尾加了归一化，把最大步长压缩到1以内，这样就算出现重尾跳跃，也控制在相对合理的范围内。

3.4 调用示例与目标函数封装

PGA主函数写好后，调用方式和粒子群几乎一样。拿经典的Sphere函数做测试：

%% 测试函数 function f = Sphere(x) f = sum(x.^2); end %% 主程序调用 clear; clc; dim = 30; lb = -100 * ones(1, dim); ub = 100 * ones(1, dim); SearchAgents = 60; Max_iter = 500; [Best_pos, Best_fitness, Convergence_curve] = PGA(SearchAgents, Max_iter, lb, ub, dim, @Sphere); fprintf('最优适应度: %.6e\n', Best_fitness); plot(Convergence_curve, 'LineWidth', 1.5); xlabel('迭代次数'); ylabel('最优适应度'); grid on;

4. 基准函数实验与参数调优建议

4.1 不同测试函数的实验现象

用Sphere、Rosenbrock、Rastrigin这三个典型测试函数分别跑实验时，观察到的收敛行为差异很大。Sphere函数是单峰光滑函数，PGA在前期就能迅速下降，50代以内基本收敛到1e-10量级；Rosenbrock函数的谷底是一条弯曲的狭长通道，PGA的蓝光策略里加入了种群差异扰动，让个体在通道中能沿着谷底方向移动，避免了粒子群算法里常见的“之”字形震荡；Rastrigin函数布满大量局部极值，这时红光策略和光竞争淘汰机制的作用就非常明显了，每10代淘汰一批弱势个体并在最优附近重播种，能够让种群反复跳出局部陷阱。

不过要提醒的是，Rastrigin这类多峰函数上，PGA也不是每次都能跑出全局最优。它的成功率大概在七成左右，和种群规模、淘汰周期、光质分配比例都有关系。实际操作时可以把种群规模放大到100，最大迭代次数提高到1000，成功率会明显上升。

4.2 关键参数灵敏度分析

我把PGA里几个关键参数的调试经验整理成了表格，方便对照：

从我自己的调试经验看，最值得优先调的是淘汰周期和光质分配比例。这两个参数直接决定了算法在后期能否稳定收敛。如果发现收敛曲线出现明显台阶式下降，通常是淘汰周期过短，种群刚聚集就被重新打散；如果曲线长期不动，则说明淘汰周期太长，种群多样性已经崩溃，重播也救不回来。

4.3 收敛曲线的判读习惯

拿到一组实验结果，我会先看前10%迭代阶段的下降斜率，这个阶段主要体现红光策略的探索效率；再看中后段是否出现持续的平台期，平台期越短说明光质分配和淘汰机制的配合越好；最后看最终收敛值能否达到目标函数在该维度下的已知最优精度量级。

一个实用的小技巧是同时绘制最优适应度曲线的y轴对数坐标。线性坐标下早期下降太快，后期微小改进几乎不可见，对数坐标能把收敛精度的变化展示得更清楚。如果对数曲线在后半段不再下降，基本可以判断算法已经收敛到当前参数配置下的极限。

5. 常见问题与避坑指南

5.1 代码运行报错排查

使用这套Matlab代码时，最常遇到的问题集中在行列向量维度不匹配上。初始化时rand(SearchAgents, dim) .* (ub - lb)，如果ub和lb是列向量，而randn(1, dim)是行向量，点乘时会触发维度错误。我的建议是统一用一行向量传入边界，也就是ones(1, dim)乘标量边界值的方式。另一个频发错误是fobj写成脚本而不是函数句柄，导致fobj(Positions(i, :))无法调用。排查方法很简单，在命令行手工执行一次feval(fobj, rand(1, dim))就能确认。

5.2 早熟收敛与多样性丢失

PGA最需要警惕的是种群多样性快速衰退。因为我给蓝光策略加了强趋光项，如果所有个体都冲着一个方向收缩，整个种群会在较短时间内挤成一团，这时候光竞争淘汰机制也起不了多大作用，因为重播种的位置也在最优附近，跳不出局部区域的引力范围。

我的应对办法是两个方向同时下手：一是把红光策略的最低比例锁死在20%，保证任何时候都有至少五分之一的个体在做全局探索；二是把重播种的扰动宽度从0.2适当提高到0.3，让新种子稍微分散一些。遇到高多峰问题，甚至可以给重播种环节加一个随机切换逻辑，有50%概率从整个搜索空间随机撒点，而不是都堆在最优附近。

5.3 实际工程应用注意事项

在实际项目里使用PGA，我最想强调的是一条经验：不要把算法参数固定死，至少要针对不同的目标函数保留几套预设参数组。比如在某个机构的配电网优化模拟项目里，我用的是种群80、淘汰周期8、淘汰比例0.25这组参数，跑一次大概需要3分钟左右；而处理一个图像处理Demo里的参数寻优问题时，同样的参数组合出现了收敛过慢的现象，换成淘汰周期15、紫光比例40%后，效果明显改善。这说明参数和问题结构是强相关的，跑实验时先做一组小规模参数扫描，比直接堆算力划算得多。

另外，工程上把PGA嵌入实际的优化流程时，建议在目标函数里加一个评估计数器。PGA单次迭代会重复评估大量个体，特别是重播种阶段会密集调用目标函数，如果目标函数本身计算量很大，整个优化耗时可能变得不可接受。加上缓存机制，对重复出现的位置直接返回历史评估结果，能省不少时间。

5.4 一套我常用的参数初始化模板

如果你想把PGA跑起来又不想从零开始调，可以照着我这套模板起步。针对30维以下的连续优化问题，种群规模取60，最大迭代次数取500，光吸收系数取0.5，生长系数取1.5，淘汰周期取10代，淘汰比例取20%，紫光策略的扰动宽度取搜索区间宽度的0.08倍。这套参数跑经典单峰函数和多峰函数都有不错的基线表现，后续再根据实际结果微调即可。

6. 关于后续扩展的想法

PGA这套框架的可扩展性很好，我最近尝试的方向是把三种光质策略替换成多种群分工：一个子种群专门负责全局探索，另一个子种群专门负责局部精修，两个子种群之间定期交换个体，模拟自然界中不同光照环境下植物的分化适应。从初步实验看，这种方式在多目标优化问题上有潜力，后续可以考虑把PGA扩展到多目标排序和拥挤度距离计算的方向上。对已经掌握基本PGA实现的读者来说，这个方向值得尝试。