遗传算法与粒子群算法在潮流计算中的Matlab实现与对比分析

做电力系统方向的人，大概率都绕不过潮流计算这道坎。最近有不少同行和学生在问遗传算法（GA）和粒子群算法（PSO）做潮流计算到底怎么实现、两者差在哪，正好我把Matlab代码实现完整跑了一遍，把两种算法的思路、实现细节、代码骨架和踩坑记录都整理出来。这篇东西适合两类人看：一是正在做本科或研究生课设、需要交一份“智能算法算潮流”作业的同学，二是想在实际项目中把智能算法与传统潮流计算结合的工程师。核心要回答的问题只有一个——同样的潮流方程，遗传算法和粒子群算法各自怎么解，谁的收敛更快，谁的解更稳，以及什么场景下你真的应该放弃牛顿-拉夫逊法改用它们。

1. 这东西到底在算什么：潮流计算与智能算法的交集

1.1 潮流计算的本质和传统解法的局限

潮流计算听起来很高大上，本质上就是求解一组非线性方程组。给定电网的拓扑结构、线路参数、节点注入功率，求出每一个节点的电压幅值和相角。对于稳态潮流分析来说，节点类型通常分为三类：平衡节点（Slack Bus，电压幅值和相角已知，负责吸收系统中的不平衡功率）、PV节点（有功功率和电压幅值给定，无功待求）、PQ节点（有功和无功注入给定，电压幅值和相角均待求）。

传统上最经典的做法是牛顿-拉夫逊法，它在初值合理的情况下二次收敛，迭代三到五次就能到很高的精度，在Matlab里用powerflow之类的函数几十毫秒就算完了一个IEEE 30节点的系统，快得离谱。既然传统方法这么好用，为什么还会有人研究GA和PSO算潮流？因为现实工程中总会碰到牛顿法搞不定的场景——初始值给得不好导致雅可比矩阵奇异、系统重负荷接近电压崩溃临界点、潮流无解时你根本不知道是“没有解”还是“迭代发散”，以及需要把更多复杂约束（比如线路潮流限额、安全稳定约束）塞进潮流模型里。在这些场景下，把潮流问题改写成优化问题，再扔给全局搜索算法，就成了一个实际可用的思路。

1.2 为什么不直接用牛顿法，非要用遗传算法和粒子群

我经常用一句话向别人解释：牛顿法是爬坡运动员，给个大致方向它能快速登顶，但遇到悬崖峭壁它直接掉下去；遗传算法和粒子群是闲不住的全山搜索队，效率不高但胜在不容易被地形骗过。

把潮流问题当成优化问题之后，目标函数通常是节点功率不平衡量的平方和，约束包括电压幅值上下限、发电机无功出力限制等。这种形式的优点是模型扩展非常方便——你想加任何实际约束，只需要往惩罚项里追加一条。缺点也明显：维度高、非线性强、多峰严重，恰好是普通梯度类方法最头疼的形态，却是GA和PSO这种不依赖梯度信息的元启发式算法的主场。

所以在后面所有讨论里，我默认采用的思路是：潮流方程不再当成等式方程组直接迭代求解，而是设计师一个残差最小的优化目标，然后用GA和PSO去逼近最优解。这个思路想明白了，代码逻辑就很顺了。

2. 两种算法怎么“插足”电网计算

2.1 遗传算法：向“物竞天择”借一套搜索逻辑

遗传算法（GA）这里我就不长篇大论背教科书了，用一个运输配送的例子解释。假设你要规划一条最短的快递配送路线，每一种方案就是一条完整的路线序列。你可以把这些路线方案看作一个“种群”，每一条路线是种群里的一个“个体”。有些路线短、耗油少，就是“适应度高”的个体；有些路线绕大圈，就是“适应度低”的个体。接下来模拟自然选择：优秀的个体获得更多“交叉”和“变异”的机会——交叉就是把两条路线的片段拼接出新的路线，变异就是把某一段顺序随机重排一下。每迭代一代，低质量的个体被淘汰，种群的平均质量逐步上升，直到逼近最优路线。

把同样的逻辑搬到潮流计算里，一个个体就是一组电压幅值和相角的组合，适应度函数就是潮流方程残差的倒数（或负值）。GA的搜索逻辑不依赖任何梯度信息，也不要求目标函数光滑，所以即便潮流方程在大范围的初值空间里呈现出高度非线性和多峰特征，GA也不会因为初值差而直接崩掉。但代价是慢，而且越到后期，通过交叉变异提升解的质量就越有限。

2.2 粒子群算法：鸟群觅食的数学化

粒子群算法（PSO）的逻辑比GA更简单直观——想象一群鸟在一片区域里找食物，每只鸟不知道食物在哪，但知道当前位置离食物有多远。每只鸟一边自己探索，一边观察群体里哪个伙伴目前离食物最近，然后不断调整自己的飞行方向和速度：既保留一点自己的惯性，又向“个人历史最佳位置”靠拢，再向“全局历史最佳位置”靠拢。三个参数——惯性权重w、个体学习因子c1、社会学习因子c2——把这三股力量拧在一起。

在潮流计算里，粒子的位置向量就是一组状态变量（电压幅值和相角），速度向量就是每次迭代状态变量的修正量。每一次飞行结束后，用潮流方程残差评估新位置的适应度，不断更新个体最优pbest和全局最优gbest。相比GA，PSO最大的特点是收敛速度肉眼可见地快，往往前期十几代就能逼近一个相当好的解。但隐患也很明显：如果群体过早地聚集到某个局部最优附近，后面的搜索基本就僵住了，表现为“早熟收敛”，这在潮流这种多峰问题上特别常见。

2.3 把潮流问题改写成优化问题的工程细节

如果你之前没写过智能算法求潮流的代码，这里有一个最容易踩坑的认知误区：不要直接把所有的电压幅值和相角全部丢给算法去搜索。因为平衡节点的电压幅值和相角是定死的，PV节点的电压幅值也是定死的，如果把它们也当作变量自由搜索，最终结果一定会偏离物理边界。

正确的做法是做“自由度归约”。设系统总节点数为n，PV节点数为npv，平衡节点数为1，则参与优化的变量数量是(n-1)个相角（平衡节点相角固定为0）加上(n-1-npv)个PQ节点的电压幅值，注意PV节点的电压幅值不参与优化。这样处理之后，算法搜索空间的维度显著降低，而且每一轮评估适应度之前，直接把定值节点赋回去，保证候选解始终落在物理可接受的范围内。

3. Matlab代码实现的完整思路（含核心代码）

3.1 节点导纳矩阵与潮流目标函数

第一次上手建议不要直接用IEEE 14或IEEE 30节点系统，先用一个小型系统把流程跑通。我用的是一个5节点系统，包含1个平衡节点、1个PV节点和3个PQ节点。第一步是把节点导纳矩阵Ybus建好，这步通常用支路参数矩阵直接组装。

function Ybus = build_ybus(branchData, busData) n = size(busData, 1); Ybus = zeros(n, n); % branchData: [from, to, r, x, b] for k = 1:size(branchData, 1) f = branchData(k, 1); t = branchData(k, 2); r = branchData(k, 3); x = branchData(k, 4); b = branchData(k, 5); z = r + 1j * x; y = 1 / z; Ybus(f, f) = Ybus(f, f) + y + 1j * b / 2; Ybus(t, t) = Ybus(t, t) + y + 1j * b / 2; Ybus(f, t) = Ybus(f, t) - y; Ybus(t, f) = Ybus(t, f) - y; end end

接下来是目标函数。核心任务是：输入一组状态变量，计算各节点注入功率的实部和虚部，与给定的P、Q期望值求差，把所有偏差的平方加总。

function f = power_flow_fitness(x, Ybus, busType, Psp, Qsp, fixedV) % x = [theta_pq; theta_pv; V_pq] % 根据优化的索引结构，恢复完整V和theta向量 n = length(Psp); V = ones(n, 1); theta = zeros(n, 1); % 把优化变量填入对应位置 V(pqIdx) = x(1:numPQ); theta(pqIdx) = x(numPQ+1:numPQ+numPQ); theta(pvIdx) = x(numPQ+numPQ+1:end); % 固定电压幅值：平衡节点与PV节点 V(slackIdx) = fixedV(slackIdx); V(pvIdx) = fixedV(pvIdx); theta(slackIdx) = 0; % 计算复功率 Vc = V .* exp(1j * theta); I = Ybus * Vc; S = Vc .* conj(I); Pcal = real(S); Qcal = imag(S); % 潮流残差平方和 f = sum((Pcal - Psp).^2) + sum((Qcal - Qsp).^2); % 可以在这里追加罚函数，处理PV节点无功越限或电压越限 end

注意看，目标函数返回的不是潮流计算通常意义上的向量残差，而是一个标量。这个标量就是GA和PSO优化的“适应度”核心——值越小，说明当前的电压状态越接近真实潮流解。当残差小于1e-6时，基本可以认为已经找到了一个高精度的潮流解。

3.2 GA求解潮流的关键参数与代码骨架

Matlab从R2010a开始提供了全局优化工具箱，内置了遗传算法函数ga，不需要自己造轮子。但必须强调的是，内置ga是一个通用框架，你要做的事情是把目标函数封装成它要求的格式。一个完整的最小实现骨架如下：

nvars = numPQ + numPQ + numPV; % 优化变量总数 lb = [0.85 * ones(numPQ, 1); -pi * ones(numPQ, 1); -pi * ones(numPV, 1)]; ub = [1.15 * ones(numPQ, 1); pi * ones(numPQ, 1); pi * ones(numPV, 1)]; options = optimoptions('ga', ... 'PopulationSize', 150, ... 'MaxGenerations', 300, ... 'MaxStallGenerations', 50, ... 'Display', 'iter', ... 'UseParallel', true, ... 'PlotFcn', @gaplotbestf); [x_ga, fval_ga] = ga(@(x) power_flow_fitness(x, Ybus, busType, Psp, Qsp, fixedV), ... nvars, [], [], [], [], lb, ub, [], options);

这里有几个我实测出来的结论想多说一句。种群规模不要盲目追求大，100到200之间就比较合适。太大了计算量翻倍但精度提升有限；太小了容易在一开始就丢掉某些关键区域的搜索机会。最大代数设成300，配合MaxStallGenerations设置50代没有明显改进就提前停止，可以省掉大量无效迭代。相角的上下界设成-pi到pi，电压幅值设成0.85到1.15，这是工程上可接受的电压偏差范围，同时也是在帮算法缩小搜索空间。

GA在这种配置下，通常到100代左右残差能进入10的负4次方量级，想进一步降到1e-8就很吃力了。这背后的原因是GA靠交叉变异做局部精调的能力比较弱，当种群整体都快收敛时，变异步长不好控制，导致最后阶段是“大炮打蚊子”。

3.3 PSO求解潮流的关键参数与代码骨架

PSO在Matlab里同样有内置函数particleswarm，从R2014b开始提供。不过我建议自己写一个实现，原因有两个：一是你很容易在简历和报告里说清楚内部逻辑，二是实现对参数的掌控更细。内置函数的骨架和GA类似：

options = optimoptions('particleswarm', ... 'SwarmSize', 100, ... 'MaxIterations', 300, ... 'InertiaRange', [0.4, 0.9], ... 'SelfAdjustmentWeight', 1.5, ... 'SocialAdjustmentWeight', 1.7, ... 'Display', 'iter'); [x_pso, fval_pso] = particleswarm(@(x) power_flow_fitness(x, Ybus, busType, Psp, Qsp, fixedV), ... nvars, lb, ub, options);

如果是自己写标准PSO，核心更新公式就是这么几行：

% 参数设置 w = 0.8; % 惯性权重 c1 = 1.5; % 个体学习因子 c2 = 1.7; % 社会学习因子 vmax = 0.1 * (ub - lb); % 速度上限 % 初始化 x = repmat(lb, nSwarm, 1) + rand(nSwarm, nvars) .* repmat(ub - lb, nSwarm, 1); v = zeros(nSwarm, nvars); pbest = x; fitness_pbest = evaluate_population(x); [gbest_val, gbest_idx] = min(fitness_pbest); gbest = x(gbest_idx, :); for iter = 1:maxIter for i = 1:nSwarm v(i, :) = w * v(i, :) ... + c1 * rand * (pbest(i, :) - x(i, :)) ... + c2 * rand * (gbest - x(i, :)); % 限速 v(i, :) = max(min(v(i, :), vmax), -vmax); x(i, :) = x(i, :) + v(i, :); % 越界处理 x(i, :) = min(max(x(i, :), lb), ub); new_fit = evaluate_single(x(i, :)); if new_fit < fitness_pbest(i) pbest(i, :) = x(i, :); fitness_pbest(i) = new_fit; end if new_fit < gbest_val gbest = x(i, :); gbest_val = new_fit; end end end

PSO前期收敛是真快，我跑到两三百个粒子、不到50代的时候残差就能压到1e-5左右，这是GA前五六十代做不到的。它的危险在于后期容易原地打转，如果惯性权重不衰减，粒子会在全局最优附近来回振荡，很难压到更高的精度；如果惯性权重衰减得太快，又会在还没搜完整个空间时就收缩到局部最优。内置particleswarm默认会在后半段逐渐收缩惯性范围，这个设计实际上很有讲究，自写实现的时候建议也做线性递减：w从0.9线性降到0.4。

4. 实测对比：GA vs PSO在哪一步拉开了差距

4.1 测试环境与对比指标

我的实验环境是Matlab R2023b，一台i7-12650H、16GB内存的笔记本。算例用的是IEEE 14节点系统，节点导纳矩阵直接读Matpower的标准数据，不做任何简化。为了公平比较，GA和PSO的种群/粒子数都设为100，迭代代数都设为300，电压上下限和相角上下限完全一致，目标函数用同一个函数文件。每组实验独立重复20次，记录每次的最优残差、达到残差1e-4所需的代数、单次运行总耗时、20次结果的标准差。这个“重复多次看标准差”的习惯很重要，因为智能算法本质上是随机搜索，只看一次结果完全说明不了问题。

4.2 收敛曲线与结果表格

直接说数据吧。20次独立运行中，GA每次都能收敛到可行解，最终那个最小残差的中位数大概是6.3e-6，最好的一次到了9.8e-7，最差的一次则在4.1e-5附近。PSO的表现有点两极分化：如果粒子群前期收到了正确的信号，50代左右残差就能到1e-5，最终甚至能摸到5e-7这个量级，比GA最好的成绩还漂亮；但20次里有4次明显早熟，群体聚集在一个残差约1e-3的局部区域出不来，表现为迭代曲线长时间水平。单次运行耗时方面，PSO大概比GA快15%左右，因为GA每一代要做的选择、交叉、变异操作本身有一定额外成本。

两种算法求解质量的典型情况：

这份数据指向一个非常典型的结论：PSO适合在有限迭代代数和精度要求不极端的情况下快速拿到一个可用的解；GA适合通过多次运行积累信心，追求稳定输出。如果只跑一次就交差，PSO的早熟风险是实实在在存在的；如果工程上允许重复运行多次取最优值，那么PSO的上限比GA要高。

4.3 参数敏感性：一个常见但容易被忽略的问题

很多同学抄完代码直接跑，觉得GA和PSO“不好用”，十有八九是参数没调对。参数敏感性的对比其实很能说明两种算法的性格差异。

GA里最影响表现的是种群规模和交叉率。我试过把种群从100降到30，最终残差直接退化了一个量级；把交叉率从0.8调到0.3，收敛速度明显变慢。但GA的好处是宽容度大——即便参数不是最优的，它通常还是能收敛到某个可用的解，只是快慢问题。

PSO的脾气就差很多了。惯性权重w从0.8改到0.6，可能就从“顺利收敛”变成“严重早熟”。self adjustment weight（个体学习因子）设得过大，粒子会过度自嗨，群体协作感丧失；社会学习因子设得过大，粒子又会疯狂涌向当前最好位置，多样性迅速崩塌。网上有很多论文给了“标准参数组合”：w=0.8，c1=c2=2，但这真的只是一个起点，具体到每一个算例都应该跑一个小参数扫描。我个人的经验是：先用w=0.9、c1=2、c2=2跑一遍，看收敛曲线有没有明显平台期，再按每次增减15%的幅度调整，记录结果，不要凭感觉拍脑袋。

5. 常见问题与排查技巧实录

5.1 明明代码逻辑没问题，结果却很不稳定

这是被问得最多的一个问题。很多人第一版代码跑出来的最优残差忽高忽低，一次到1e-8，下一次直接卡在1e-2。排查思路按以下顺序来。

先看随机种子。GA和PSO都是随机算法，不加rng(seed)控制的话，每次运行的初始种群完全随机，结果波动很正常。这不算bug，但如果你要对比两者的性能，必须在完全相同的数据和初始化条件下测试，并且保证两种算法使用同样的随机流。

再看速度限制和边界处理。PSO粒子飞出边界后会怎样处理，直接影响搜索稳定性。最简单的钳位法（超出边界就按边界值算）和反射法（超出边界就按镜像位置弹回）表现差异不小，我实测反射法在电压幅值这种有物理边界的变量上效果更好，因为它让粒子有机会在边界附近“滑行”，而不是全部堆死在边界上。GA里则要注意变异算子的选择，如果变异步长固定为0.01，前期搜索和后期精调用的都是同一种扰动强度，后期收敛自然很不稳定。使用高斯变异、变异幅值随代数线性收缩，能显著改善末期稳定性。

最后别忘了检查自适应度函数本身有没有bug。我有一个让人头秃的经历——优化变量和节点恢复矩阵的索引没对齐，导致电压幅值错位，结果算法一直在追逐一个根本不合理的解。这种问题用视觉化排查最容易发现：画一次电压幅值分布图，如果出现相邻节点电压一个在1.1、一个在0.9这种怪异的锯齿形态，基本可以断定是索引问题而不是算法问题。

5.2 罚函数权重该怎么选

把潮流等式约束放进目标函数之后，你通常还要考虑不等式约束——例如PV节点的无功出力越限、PQ节点的电压幅值越限。处理不等式约束最朴素的办法是加罚函数项，也就是在目标函数后面加上一个“越限程度乘以权重系数”的惩罚项。这个权重系数的选取其实相当微妙。

权重设小了，算法会无视约束，最终解电压0.5也能被接受；权重设大了，又会喧宾夺主，算法拼命压低约束违反量而忽略了潮流残差，最终可能得到一个不满足潮流方程的方案。我试过几种定权思路后，比较推荐“分级递增”的策略：前几十代用一个相对小的权重，让搜索尽量广地探索；在最后几十代把权重逐步调大，逼迫算法把解拉回到可行域内。这个方案在实际算例中效果不错。

另外强烈建议不要在第一次迭代就把罚函数加上。先用纯等式约束跑通GA/PSO，确认基架没问题之后，再把不等式约束的惩罚项加上去。分模块调试，是我这么多年调智能算法养成的最重要的习惯。

5.3 边界越界、粒子飞出可行域怎么办

很多人设定了lb和ub，但实际搜索过程中仍然会出现变量取值超出边界。这往往不是因为边界设定本身有问题，而是因为边界钳位之后没有同步考虑变量的物理含义。比如你把电压幅值限制在[0.85, 1.15]，但在这个区间里依然存在大量电压组合不满足PV节点的无功约束或线路潮流约束，此时解虽然在“边界内”，却在物理上不可行。

更隐蔽的一种越界发生在PSO的速度更新里：如果vmax设得太大，粒子一次飞行就从1.15冲到了0.7，即便你后续做钳位，粒子的飞行轨迹已经穿过了大片不可行区域，搜索效率大打折扣。vmax的经验取值是变量范围的10%到20%。例如电压幅值宽度为0.3，vmax就设0.03到0.06。这个细节我至少见过三四个项目圈进去了，代码看着没问题，就是不收敛，把最大飞行速度降下来之后马上就顺了。

5.4 有没有必要GA+PSO混合用

很多论文和博客都提到混合策略，这里说点不一样的。GA和PSO的优势劣势是互补的：GA全局探索能力好、后期稳定，但前期进展慢；PSO前期收敛快、局部开发能力强，但容易早熟。混合的思路无外乎两种。

一种是“接力式”——先用GA跑几十代找到一个有潜力的区域，把得到的gbest作为PSO的初始粒子，再让PSO快速精调。这种组合在IEEE 14节点上确实比任一单算法表现更好，代价是实现复杂度上升、调试变量翻倍。另一种是“并行式”——种群内一部分个体按照GA逻辑更新，一部分按照PSO逻辑更新，信息共享相同的pbest/gbest。这种思想在很多论文里被包装成“混合粒子群遗传算法”，实际工程价值见仁见智。

如果只是做课设或者验证算法性能，我不建议一上来就搞混合。先把两种单一算法的参数调明白，再考虑组合。连单一算法的行为都摸不清楚，混合之后出了任何异常你都无从下手去定位问题。

6. 一点使用心得

说句实在话，用遗传算法或粒子群算法替代牛顿-拉夫逊法去做标准系统的常规潮流计算，效率上没有任何优势，甚至可以说是杀鸡用牛刀。但如果你面临的是一个高维度、强约束、对初值极度敏感的复杂潮流问题，或者在潮流计算基础上还要做安全约束优化、无功优化、可用输电能力评估，那么GA和PSO的灵活性就体现出来了——毕竟你只需要改动目标函数和约束条件，不用重新推导任何偏导数和雅可比矩阵。我个人的经验是：在一开始写代码的时候就把目标函数封装成接口，案例测试用GA和PSO分别跑通，之后再嫁接各种约束和改进策略就会很顺。对比实验请一定记录不同随机种子下的多次运行结果，不要只看单次最优值，这个习惯能让你少写一整篇“让人困惑的结论”。最后再分享一个小技巧：如果你只是想快速验证自己的想法，先在Matlab里用ga和particleswarm的内置实现跑通流程，再考虑自写改进算法，站在封装好的轮子上出发，能省下不少调试时间。