粒子群算法求解配电网无功优化：从IEEE33节点到MATLAB实现全流程实践

早在IEEE33节点配电网这个经典算例上跑通粒子群无功优化之前，我其实在这类项目上栽过不少跟头。最初接触“基于粒子群算法的配电网无功优化”这个题目时，我以为难点在算法推导，结果真正动手才发现，建模、编码、边界处理、潮流计算，每一步都有坑。这篇文章不是教科书式的复述，而是把我从方案选型、IEEE33节点参数整理、PSO参数调优，到MATLAB代码实现、结果分析、常见报错排查的完整过程梳理一遍。想用智能算法做配电网无功补偿的同学，或者刚接触优化算法与电力系统结合课题的工程师，应该能从这里拿走一套可以直接复用的解决思路。

先说清楚这个东西能干什么：在IEEE33节点标准配电网模型上，通过粒子群算法自动搜索无功补偿装置的安装位置和容量，让系统网络损耗最小、节点电压恢复到合格范围。很多人觉得无功优化就是个“加电容”的事，但实际求解时你会面临补偿点怎么选、装多大容量、怎么跟潮流方程耦合、怎么处理离散变量这一连串问题。这篇文章就围绕这些实际细节展开，不讲空话。

1. 为什么做这件事：无功优化的实际意义与方案选型

1.1 无功补偿到底在补什么：线损和电压的联动关系

配电网和输电网的最大区别，在于电压等级低、线路阻抗相对较大、负荷点多而分散。到了傍晚用电高峰，馈线末端的电压经常掉到0.90 p.u.以下，这时候线损率高得吓人。很多人搞不清楚无功功率对线损的影响机制，其实用一个简化公式就能说明白。

一条支路的有功损耗可以写成：

( \Delta P = \frac{P^2 + Q^2}{U^2} R )

这里P是有功功率，Q是无功功率，U是节点电压，R是线路电阻。从公式可以直观看到，无功功率Q在线路上的传输同样会占用线路容量、产生有功损耗。如果能在负荷附近就地补偿无功，那么线路上流动的Q就会大幅减少，(\Delta P)中的(Q^2 R / U^2)这一项会显著下降。

举个实际例子：一条10kV馈线，电阻R=2Ω，输送无功功率Q=1500kvar，电压U=10.5kV时，仅无功分量造成的线损大约有 (1500^2 \times 2 / 10.5^2 \approx 40.8kW)。如果把无功就地平衡掉，这部分损耗直接就消失了。再加上电压抬升后，电流减小、变压器和线路损耗进一步降低，最终网损降幅十分可观。

所以无功优化本质上是在解决两个问题：一是把系统网损降下来，二是把电压质量提上去。补偿装置（典型为并联电容器组）投入的位置和容量不同，效果差异巨大。装错了地方不仅浪费投资，还可能造成无功倒送、电压过高等新问题。这就是为什么需要优化算法，而不是拍脑袋在末端怼几组电容。

1.2 传统方法的瓶颈与粒子群算法的优势

无功优化在数学上是一个典型的非线性、多约束、混合整数优化问题。决策变量既包含连续量（无功注入），又包含离散量（电容器组数），目标函数是潮流方程的函数，无法直接写出显式表达式。

传统方法比如等网损微增率法、线性规划法、非线性规划法，在处理这类问题时条件很苛刻。等网损微增率法对网络结构有较强依赖；牛顿法和内点法需要计算梯度或海森矩阵，碰上离散变量和强非线性约束时容易发散；而且它们大多是单点搜索，一旦初始点选得不好，结果往往收敛到局部最优。

粒子群算法就不太在乎这些问题。它的优点非常明显：

• 不需要目标函数可导，适应度函数只要能通过潮流计算算出数值就行。
• 群体搜索，天然具备全局探索能力，不容易陷入局部最优点。
• 参数少、结构简单，用MATLAB几百行就能实现。
• 对连续和离散混合变量处理起来比较方便，只需在编码和更新环节稍作处理。

我用粒子群算法和传统非线性规划方法在同一套IEEE33节点系统上做过对比。粒子群算法在网损下降幅度和电压改善程度上都要优于传统方法，尤其是在补偿点数量较多时，PSO能更好地利用群体信息协同搜索。

1.3 各类智能算法的横向对比与选题考量

还有同学问我，既然智能算法这么多，遗传算法、差分进化、蜂群算法也能做啊，为什么非选粒子群？

这个问题问得很有价值。我把几种常见算法的特点列一下，方便大家做选型参考。

最终选PSO，是因为它对无功优化这个场景的匹配度最高。本项目里决策变量维度并不高，不超过10个补偿节点，PSO在这个维度区间收敛速度明显快于GA，实现复杂度又低于DE和ABC。而且PSO在迭代前期因为有群体信息交互，能快速锁定最优解的区域，这对后续调参和结果分析都很友好。

2. IEEE33节点配电网建模与无功优化问题描述

2.1 IEEE33节点系统参数与拓扑结构：一份必须核对的数据清单

IEEE33节点系统是配电网研究里最经典的辐射状算例，由一个根节点（通常编号0或1）通过33个节点、32条支路呈辐射状网络向外供电。基准电压取12.66kV，基准功率取10MVA，系统总负荷约3715kW+2300kvar。

这里要特别强调：网上能找到的IEEE33节点参数版本很多，有的节点编号从0开始，有的从1开始，支路参数矩阵格式也有差异。我第一版程序就因为在导入数据时把节点编号错位，导致潮流结果完全对不上。后来我花了一个晚上把支路数据和负荷数据手工对齐，才算把坑填上。

支路参数表只截取前几行给大家感受一下格式：

注意最后一个节点的负荷是60kW/20kvar，而总负荷3715kW/2300kvar是由所有节点负荷之和验证出来的。你拿到数据后，第一步就要算一遍总负荷，如果对不上，后面的所有结果都不可信。

基准值换算也是个常见坑。潮流计算里如果使用标幺值，需要把实际线路阻抗除以基准阻抗。12.66kV系统、10MVA基准下，基准阻抗是 (Z_{base} = U_{base}^2 / S_{base} = 12.66^2 / 10 = 16.0249\Omega)。如果不做换算直接用有名值计算，后面的前推回代迭代极容易不收敛。

2.2 目标函数、约束条件和决策变量的数学表述

无功优化的目标是综合网损和电压质量。最常见的目标函数是：

( \min F = P_{loss} + \lambda \sum_{i=1}^{n} (\Delta U_i)^2 )

其中 (P_{loss}) 是系统总有功网损，(\Delta U_i) 是节点i的电压偏差标幺值，(\lambda) 是电压惩罚系数。(P_{loss}) 的计算由潮流结果得出，是各支路损耗之和。

等式约束其实就是潮流方程本身，即每个节点必须满足有功、无功功率平衡，这个在前推回代潮流中会被自动满足，不需要显式作为约束处理。

不等式约束包括：

• 节点电压约束：(0.95 \le U_i \le 1.05)
• 无功补偿容量约束：(0 \le Q_{c,i} \le Q_{c,i}^{max})
• 补偿节点数量约束：总补偿节点数不超过预设上限

决策变量就是各个候选补偿节点的无功补偿容量。这里有个典型的“矛盾”点：如果允许所有33个节点都装补偿，那决策变量就是33维，PSO搜索空间太大、收敛慢；如果只锁定三五个节点，又可能漏掉最优解。我在项目中采用了两步策略：先用电压灵敏度分析初步筛选电压偏低的末端节点作为候选补偿点，再让PSO在候选集内搜索最优容量组合。该方法既能压缩搜索空间，又能保证解的质量。

常用候选补偿节点可以选择18、17、32、25、33这类末端或电压偏低区域。实测在这些节点上做补偿，电压改善效果最明显。

2.3 补偿容量离散化与PSO编码策略

实际电容器组是按标准容量分档投切的，比如每档50kvar或者100kvar，不可能连续调节。这个特性必须反映到优化模型中，否则求出来的“0.173Mvar”在实际工程里没法落地。

我的做法是把每个补偿候选点的容量设置为0到400kvar之间、步长20kvar的离散值。粒子位置向量在更新时保持连续值，但在计算适应度时对粒子位置进行就近取整，转换为实际可投切的电容器组数，再乘以单组容量得到实际补偿无功值。

向量编码大致长这样：

[ x_i = [cap_1, cap_2, cap_3, ..., cap_m] ]

其中m是候选补偿节点数量，cap值表示该节点投入的补偿容量。粒子群更新得到的是连续位置值，再通过取整函数转换成实际容量。这种编解码方式既保留了PSO的搜索连续性，又让最终答案符合工程可行性。

适应度函数内部调用潮流计算时，需要把补偿容量以“负的无功负荷”叠加到对应节点上，即该节点的无功负荷变为 ((Q_{load,i} - Q_{c,i}))。

注意：如果补偿容量大于节点无功负荷，就会造成无功倒送。在IEEE33节点这种小系统里，倒送可能会让某些节点电压升高超过1.05 p.u.，目标函数里的惩罚项这时候就起作用了。

3. 粒子群算法求解无功优化的关键细节

3.1 从鸟群觅食到数学公式：PSO工作原理

粒子群算法的思想源于对鸟群觅食行为的模拟。鸟群在寻找食物时，每只鸟有自己发现的最佳位置，同时群体共享信息知道整个鸟群目前找到的最佳位置，个体在飞行时同时受到“自己想去的方向”和“同伴推荐的方向”共同影响。

在无功优化里，每只“鸟”就是一个候选补偿方案，每只鸟的位置就是一个补偿容量向量，适应度就是该方案下的网损和电压惩罚值。

速度更新公式：

( v_i(t+1) = w v_i(t) + c_1 r_1 (pbest_i - x_i(t)) + c_2 r_2 (gbest - x_i(t)) )

位置更新公式：

( x_i(t+1) = x_i(t) + v_i(t+1) )

式子中w是惯性权重，控制上一代速度对当前速度的影响；c1是自我认知学习因子，c2是社会学习因子；r1和r2是[0,1]之间的随机数；pbest是粒子自身历史最优位置，gbest是整个群体的历史最优位置。

这三个项的分工很有意思。(w v_i(t)) 是粒子的“惯性”部分，它保持当前运动方向，负责开发；(c_1 r_1 (pbest_i - x_i)) 把粒子拉向自己的历史最优，负责个体记忆；(c_2 r_2 (gbest - x_i)) 把粒子拉向群体最优，负责协作。三股力量相互制衡，让粒子既不会漫无目的地乱飞，也不会过早抱团失去探索能力。

3.2 惯性权重和学习因子的调参心得

这部分是PSO能不能收敛好、收敛快的关键所在。我跑了大量实验后得出的经验是：

• w大则全局搜索能力强：w接近0.9时，粒子受自身惯性影响大，会跑得很远，适合前期快速扫描整个解空间。
• w小则局部开发能力强：w接近0.4时，粒子运动速度迅速衰减，在局部精细搜索。
• 线性递减策略最实用：从0.9线性降到0.4，让算法前期“广撒网”，后期“精打细磨”。

[ w = w_{max} - (w_{max} - w_{min}) \times \frac{t}{T_{max}} ]

其中t是当前迭代次数，Tmax是最大迭代次数。

学习因子方面，经典取法c1=c2=2，实测效果稳定。也有人做改进，比如让c1从2.5降到0.5、c2从0.5升到2.5，前期的“自我探索”多些，后期的“社会跟随”多些。我实际对比下来，对于IEEE33节点的无功优化问题，经典参数加线性递减w，已经足够得到很好的结果，没必要一开始就上花哨的改进策略。

vmax也要重点提一下。如果不对最大速度做限制，粒子容易一次迭代跑出可行域十万八千里，后面再飞回来就难了。我一般把vmax设为变量取值范围的10%~20%。比如补偿容量0到0.4Mvar，速度上限就取0.04到0.08之间。

3.3 工程落地时的三个关键细节

第一，随机种子影响结果。PSO每次运行的结果都不完全一样，因为r1和r2是随机数。不能跑一次看到好结果就下结论。我在实际项目中会固定随机种子跑10次，记录平均值和标准差，这样评估算法性能才有可靠性。

第二，越界处理。粒子位置超出约束边界时，直接用边界值截断。有的教程喜欢用“边界反弹”策略，让粒子被弹回可行域内。但反弹策略在靠近边界的最优解附近容易造成反复震荡，实际效果不如直接截断稳。

第三，离散量的取整时机。如果每个粒子的适应度计算都先对位置取整再算潮流，会导致多个连续位置映射到同一个离散档位上，适应度函数变成阶梯函数。粒子在平面上移动时可能陷入“平台区”，不知道该往哪走。解决办法是让位置更新始终在连续空间进行，只有真正计算实际补偿方案时才取整，算法收敛后再把最终gbest位置取整输出。

4. 实操过程与结果展示

4.1 前推回代潮流计算子程序

IEEE33节点是辐射状配电网，用前推回代法计算潮流比牛顿-拉夫逊法简单一个数量级。这种方法的思路很直观：从末端节点开始向根节点回推，计算各支路流过的功率；再从根节点向末端前推，算出各节点电压。反复迭代直到前后两次电压差值满足精度要求。

前推回代算法步骤整理如下：

1. 初始化所有节点电压为1.0 p.u.。
2. 回代：从末端支路开始，根据节点负荷和下游支路功率累加出当前支路末端流出的功率。
3. 前推：从根节点开始，用支路首端电压减去压降得到支路末端电压。
4. 收敛判断：前后两次迭代的电压最大差值小于1e-6时停止。

核心MATLAB代码可以按这个思路写：

function [V, Ploss] = my33BackwardForward(branch, load, Qc) % branch: 支路矩阵，每行[首端节点, 末端节点, R, X] % load: 节点负荷矩阵，每行[P, Q]，下标从1开始 % Qc: 节点无功补偿容量向量 n = size(load, 1); V = ones(n, 1); Q = load(:, 2) - Qc(:); iter = 0; while iter < 100 V_old = V; % 回代：从支路末端向前累加功率 S = zeros(n, 1); S(:, 1) = load(:, 1) + 1j * Q; for k = size(branch, 1):-1:1 e = branch(k, 2); s = branch(k, 1); S(s) = S(s) + S(e); % 加上支路损耗 dl = (abs(S(e)) / V(e))^2 * (branch(k, 3) + 1j * branch(k, 4)); S(s) = S(s) + dl; end % 前推：从根节点往末端更新电压 for k = 1:size(branch, 1) s = branch(k, 1); e = branch(k, 2); I = conj(S(e)) / conj(V(s)); V(e) = V(s) - I * (branch(k, 3) + 1j * branch(k, 4)); end if max(abs(abs(V) - abs(V_old))) < 1e-6 break; end iter = iter + 1; end % 计算网损：所有支路损耗之和 Ploss = 0; for k = 1:size(branch, 1) s = branch(k, 1); e = branch(k, 2); I = (V(s) - V(e)) / (branch(k, 3) + 1j * branch(k, 4)); Ploss = Ploss + abs(I)^2 * branch(k, 3); end end

这里用的是标幺值计算，支路阻抗需要提前除以基准阻抗。如果直接用法，电压收敛精度可以到1e-6，完全满足后续优化迭代的要求。

4.2 PSO主程序与适应度函数

主程序负责调度整个求解流程，包括数据导入、参数设置、种群初始化、迭代搜索、结果输出。提前说一句：这里给的是结构清晰的参考代码思路，实际使用时要按照自己定义的变量类型做适配。

% 参数设置 popsize = 30; % 种群规模 maxIter = 100; % 最大迭代次数 w_max = 0.9; w_min = 0.4; c1 = 2.0; c2 = 2.0; dim = numel(candidateNodes); % 补偿节点数 xmin = zeros(1, dim); % 补偿容量下限 xmax = 0.4 * ones(1, dim); % 补偿容量上限 vmax = 0.08 * ones(1, dim); % 初始化粒子位置和速度 x = repmat(xmin, popsize, 1) + rand(popsize, dim) .* repmat(xmax - xmin, popsize, 1); v = -vmax + 2 * vmax .* rand(popsize, dim); % 计算初始适应度 for i = 1:popsize fitness(i) = calcFitness(x(i, :), branch, load); end pbest = fitness; pbestX = x; [gbest, idx] = min(pbest); gbestX = pbestX(idx, :); % 主迭代 for t = 1:maxIter w = w_max - (w_max - w_min) * t / maxIter; for i = 1:popsize v(i, :) = w * v(i, :) + c1 * rand(1, dim) .* (pbestX(i, :) - x(i, :)) ... + c2 * rand(1, dim) .* (gbestX - x(i, :)); % 速度越界处理 v(i, :) = min(max(v(i, :), -vmax), vmax); % 位置更新与越界处理 x(i, :) = x(i, :) + v(i, :); x(i, :) = min(max(x(i, :), xmin), xmax); % 计算适应度 fitness(i) = calcFitness(x(i, :), branch, load); % 更新个体最优 if fitness(i) < pbest(i) pbest(i) = fitness(i); pbestX(i, :) = x(i, :); end end % 更新全局最优 [gbest, idx] = min(pbest); gbestX = pbestX(idx, :); % 记录收敛曲线 convergence(t) = gbest; end

适应度函数内部做的事情是：把粒子的连续位置向量取整为离散容量档位，然后叠加到候选补偿节点上，调用潮流计算函数，最后返回网损与电压惩罚之和。一个容易忽略的点是，如果某节点电压越限，惩罚项必须给得足够大，否则粒子会肆无忌惮地让个别节点电压超限来换取网损微降。

function f = calcFitness(x, branch, load) % 离散化处理 Qc = round(x / 0.02) * 0.02; % 初始化补偿向量 Qc_vec = zeros(size(load, 1), 1); Qc_vec(candidateNodes) = Qc; % 潮流计算 [V, Ploss] = my33BackwardForward(branch, load, Qc_vec); % 电压越限惩罚 V_abs = abs(V); penalty = sum((min(V_abs, 0.95) - 0.95).^2) + sum((max(V_abs, 1.05) - 1.05).^2); f = Ploss + 100 * penalty; end

这里的candidateNodes是预先选出的候选补偿节点编号数组。需要提前声明为全局变量或通过闭包方式传给适应度函数，否则MATLAB会报未定义变量错误。

4.3 典型优化结果与收敛性、稳定性分析

我在IEEE33节点系统上按上述方法完成求解后，得到了一组比较有代表性的结果，这里展示给大家参考（由于PSO有随机性，数值会有少量波动，但在合理范围内）。

先看不加补偿时的潮流结果：系统总有功网损为202.68kW左右，这个数值与公开文献中的IEEE33节点常规潮流结果一致，可以用这个值来验证潮流程序是否正确。

采用4个候选补偿节点、最大单点补偿容量400kvar、步长20kvar、种群30、迭代100次的条件下，我得到的优化结果如下：

最佳补偿方案中，补偿点主要集中在电压偏低的末端区域，其中节点18的补偿容量最大，节点17、32、33也都有不同容量的补偿配置。这个结果符合工程直觉：末端节点离电源远，电压支撑最弱，无功补偿的边际效益最高。

电压分布对比也很典型。优化前，越靠近末端电压下降越明显，节点18附近电压只有0.9131，已经明显越限；优化后，全网络节点电压都回到0.95 p.u.以上，电压曲线整体抬升且更加平缓。

再看收敛曲线。典型运行中，前20代群体最优适应度快速下降，从202kW左右降到160kW附近；30~50代之间继续缓降；到60代左右基本收敛到151kW附近，之后不再有明显变化。这说明100次迭代对这个规模的问题来说是足够的，再增加迭代次数收益有限。

稳定性方面，我连续运行10次，记录最终网损最优值在150.5kW到153.2kW之间波动，波动幅度不到2%。这个稳定性说明PSO对这个问题的搜索能力是可靠的，没有出现有的批次严重陷入局部最优的情况。

5. 常见问题与避坑日记

5.1 潮流不收敛？先从三处排查

前推回代潮流不收敛或无思路，通常不是算法问题，而是数据问题。我自己的排查顺序如下：

第一，检查支路数据的节点编号顺序。IEEE33节点是从0开始编号还是从1开始编号，直接决定程序和数据的对应关系。如果根节点不是0或1，或者支路首末端节点定义反了，潮流一定会发散或得到离谱结果。

第二，检查基准阻抗换算。用标幺值计算时必须把有名值除以基准阻抗。如果忘了这一步，电阻电抗会大几十倍，迭代中电压会出现明显波动，甚至越迭代越发散。

第三，检查负荷数据的维度对齐。节点负荷必须和节点编号一一对应，某个节点漏了或重复会导致功率不平衡。以前我就因为复制粘贴时把两个节点的负荷弄混，导致潮流结果比真实值差出一个数量级。

判断潮流程序是否正确的标准很简单：不投入补偿时，IEEE33节点的总网损应当稳定在202.68kW左右。如果偏差超过5%，优先怀疑前面三类问题。

5.2 PSO搜索散漫、收敛曲线不降怎么办

这是最让人抓狂的情况。明明程序没报错，但收敛曲线一路平坦，网损降不下去。我做过的有效调整如下：

首先看候选补偿节点选得对不对。如果选的节点全在电源附近，电流本来就小，无功补偿意义不大，网损自然降不下去。我一开始把所有候选节点放在前10个节点里，结果网损只降了几个千瓦，后来改成末端节点，效果立刻不一样了。

然后检查vmax是否过大。vmax如果设置成和变量范围一样大，粒子的位置更新会像“随机乱跳”，完全丧失局部搜索能力。把vmax降到变量范围的15%左右，收敛通常会明显改善。

最后检查适应度函数里是否真的反映了补偿效果。有一种情况是代码里Qc没有正确叠加到潮流计算中，PSO再怎么搜，目标函数都不变，收敛曲线自然是一条水平线。这种“隐式bug”很难发现，需要人为打印补偿容量代入后的节点电压和网损，确认是否随补偿变化。

5.3 结果反复横跳、多次运行差异大

如果连续运行多次得到的gbest差异超过5%，说明算法稳定性不够。我从实际经验总结了几种缓解手段：

• 增加种群规模。30个粒子在4维优化问题上没问题，但如果候选节点增加到8个以上，种群最好加到50以上。
• 增加迭代次数。后期收敛曲线还在下降就说明还没收敛，增加迭代数看变化。
• 引入早停机制。如果gbest连续20代没有更新，可以提前终止，避免浪费时间。
• 多次运行取最优。工程上最终方案建议取10次运行中网损最小的一次，并验证该方案下电压约束全部满足。

另外，补偿容量的离散化步长也会影响结果稳定性。步长取得太大，比如100kvar一档，算法在几个离散档位之间切换，收敛曲线容易出现锯齿状波动；步长太小，比如5kvar一档，搜索空间又会变得很大。20kvar是我在IEEE33节点系统上试下来比较合适的折中点。

5.4 那些不起眼但让人卡壳的小问题

MATLAB环境下的几个高频问题也值得多说一句：

• 脚本和函数命名冲突，比如脚本叫PSO.m，函数也叫PSO，会导致调用混乱。尽量给主脚本一个独立名字，函数名保持唯一。
• 候选补偿节点编号与矩阵下标混用。IEEE33节点有节点0，但MATLAB矩阵下标从1开始，手动转换时容易差一位。建议用一个数组保存“节点编号→矩阵索引”的映射关系，所有涉及节点编号的地方统一查表转换。
• 保存在循环里被覆盖。如果每次迭代只保存gbest数值而不保存gbestX向量，最后拿到最优适应度却丢了对应的方案。我习惯在每次gbest更新时同时记录gbestX和对应的补偿方案明细表。

结尾：一点个人复盘

这个项目做下来，我最深的感受是：优化算法本身不是最大的门槛，建模的细致程度才决定结果优劣。粒子群算法一遍就能跑通，但把IEEE33节点数据、前推回代潮流、离散化补偿策略、惩罚函数这几个模块可靠地衔接到一起，才是我真正花时间的地方。如果后面你想在这个基础上继续扩展，我建议优先尝试三个方面：一是把单目标换成“网损+电压偏差+补偿投资费用”的多目标优化；二是在系统中接入分布式电源，观察分布式电源无功调节与传统电容器的协同效果；三是把静态优化改成基于时序负荷曲线的动态无功优化。这套PSO框架在这些方向上改起来并不困难，核心的粒子更新机制和适应度计算流程都能复用。最后再分享一个实操小技巧：每次跑完优化，把最佳补偿方案手动代入潮流，重新算一遍电压分布和网损，核对结果和程序输出一致。别嫌麻烦，这个习惯帮我揪出过不少隐蔽的逻辑错误。