把“植物怎么晒太阳”变成一套优化算法这个点子我第一次看到时觉得挺新鲜后来动手复现了Matlab版本发现它不仅思路有趣代码结构也适合拿来当学习智能优化算法的入门模板。今天想聊的就是这个算法——阳光生长优化算法Polychromatic Glow Optimization Algorithm简称PGA以及它对应的Matlab实现细节。这个算法的中文名很形象英文全称里Polychromatic Glow指的是“多色辉光”放到算法语境里其实就是模拟植物在多色光源照射下完成光合作用、向光生长和群落光竞争的过程。和粒子群、遗传算法这类大家耳熟能详的元启发式方法相比PGA的思路更贴近生态学用不同波长的光代表不同的搜索策略用植物茎秆的伸长方向代表解位置的更新方向用光竞争机制来处理种群淘汰和多样性维护。文章后面我会把核心原理、数学公式、完整的Matlab代码框架、参数调优心得和典型问题全部拆开讲适合刚接触智能优化算法的读者也适合想找一个容易改造的算法框架用来做工程优化的开发者。1. 算法背景与灵感来源1.1 植物光照策略与优化问题的映射关系先说算法灵感。植物生长离不开光但不同波长的光对植物的作用千差万别红光能促进茎的伸长和光合产物积累蓝光能刺激叶绿素合成和细胞分化而紫光或紫外波段的光往往参与调节植物的次生代谢和抗逆响应。植物在光照环境下并不是被动接受阳光而是会主动调整姿态——向光性生长就是典型的例子哪里光强茎秆就往哪里长同时还会通过改变叶片的朝向去争夺更好的光照位置。把这个现象映射到优化问题里思路就打开了。解空间可以看成一片“森林”每个候选解就是一株“植物”适应度值就是植物吸收光能后积累的“生物量”。全局最优解是光强最大的位置其他个体要不断向这个光源靠拢。可是如果所有个体都只顾着冲向当前最优种群很快会挤在一起多样性丧失后算法就容易早熟。于是自然界的“多光谱”策略就派上了用场让不同状态的个体使用不同的光质策略有的负责大步探索新区域有的负责在小范围内精修有的负责跳出去寻找新的光照缝隙。这样设计的好处在于它不是靠单一公式从头到尾硬跑而是把搜索过程拆成了拥有明确任务分工的协作模式。红光策略承担探索任务对应算法里的大范围游走蓝光策略承担开发任务对应向最优解收缩紫光策略承担精细调节任务对应局部微调。三者比例还会随种群状态动态变化避开“一刀切”带来的早熟或停滞。1.2 PGA算法总体流程设计整个PGA算法的主流程并不复杂可以用下面几条写清楚随机撒播初始种群把每个个体当作一粒种子。计算每个个体的适应度记录当前全局最优光源位置。根据某种策略为每个个体分配一种光质红光、蓝光或紫光。个体按照对应光质的位置更新公式移动形成一次“生长”。对越界个体做边界约束处理计算新适应度按接受准则决定是否保留。每隔固定代数执行一次光竞争淘汰让弱势个体在最优光源附近重新播种。迭代至最大代数输出全局最优解和收敛曲线。从工程角度看这个流程和大多数群体优化算法兼容性很好不需要额外引入复杂的数据结构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(t1) X_{best}(t) L \cdot (X_{best}(t) - X_i(t)) ]这里 ( L ) 是Levy飞行随机向量。Levy飞行的特点是步长分布存在重尾偶尔会出现大步长跳跃非常适合在搜索空间中制造“变异”避免种群陷入局部极值。蓝光策略模拟正常光合作用下的趋光生长个体朝最优光源方向收缩同时引入种群中其他个体的差异信息来保持多样性[ X_i(t1) 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(t1) 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里几个关键参数的调试经验整理成了表格方便对照参数影响对象建议范围备注种群规模全局探索能力和计算成本20~100维度越高越取大值最大迭代次数收敛深度300~1000可视目标函数复杂度增加光吸收系数alpha光通量随距离衰减速度0.1~1.0值越大远距离个体越趋向探索生长系数beta蓝光策略扰动幅度1.0~2.5过大容易发散过小收敛慢淘汰周期多样性恢复频率5~15周期过短会导致振荡淘汰比例多样性恢复强度0.15~0.3比例过大会丢失有效信息紫光扰动宽度精细搜索步长区间宽度的0.05~0.1倍后期决定收敛精度从我自己的调试经验看最值得优先调的是淘汰周期和光质分配比例。这两个参数直接决定了算法在后期能否稳定收敛。如果发现收敛曲线出现明显台阶式下降通常是淘汰周期过短种群刚聚集就被重新打散如果曲线长期不动则说明淘汰周期太长种群多样性已经崩溃重播也救不回来。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实现的读者来说这个方向值得尝试。