前阵子处理一个电压质量分析项目，调度那边给我一张IEEE14节点系统的负荷数据，让我评估无功配置的优化空间。系统里有五台发电机、三台有载调压变压器，表面上无功容量足够，但高峰时段部分线路的有功损耗明显偏高，末端节点电压也压在下限附近。这就是电力系统里典型到不能再典型的无功优化问题：在潮流方程和各类设备约束之下，通过调节发电机机端电压、变压器变比和电容器投切容量，让系统有功网损尽可能小。我最终用粒子群算法配合Matlab/Matpower在IEEE14节点系统上把这个问题跑通了，降损效果和电压改善都比较明显。这篇博文就把整个思路、代码和数据处理的细节完整记录下来，适合做课程设计、论文复现或者刚接触电网优化的读者。

1. 无功优化到底在优化什么：从调度视角的需求拆解

1.1 线损、电压质量与设备寿命：三个被忽视的成本项

电网里的有功网损不是个小数目，尤其是轻载时段如果无功就地补偿不足，大量无功功率会在线路上来回流动，白白吃掉一部分有功容量。在IEEE14节点这种小系统上，初始网损大概在13.4MW左右，看起来不多，换算到实际大电网就是每年上千万度电的成本。调度员日常调节无功，首要目标就是把这部分损耗压下来。

第二个是电压质量。负荷节点电压偏低，电动机转矩下降、照明设备输出不足；电压偏高，变压器和电缆的绝缘老化会加速，电容器也容易过压损坏。我做项目时见过不少末端电压长期徘徊在0.95p.u.以下的变电站，用户反复投诉，但补偿容量就是投不进去——问题往往出在“不知道怎么组合最优”。

第三个容易被忽略的是设备动作次数。电容器组和变压器分接头都是有机械寿命的，频繁投切和挡位调整会显著缩短检修周期。所以在优化模型里通常要限制变量变化步长，或者尽量让解落在稳定的补偿区间。这也解释了为什么实际工程中的无功优化往往不是简单的单目标问题，不过作为算法研究，先把“网损最小+电压合格”这套单目标框架跑通，是最稳妥的起点。

1.2 为什么这个优化问题这么“棘手”

无功优化在数学上是一个带约束的非线性混合整数规划问题。难点主要有三个：

• 潮流方程是非线性的，目标函数没有闭式表达式，只能在给定控制变量下通过潮流计算得到网损值，优化和潮流是“一个萝卜一个坑”的关系。
• 控制变量里有连续量（发电机机端电压），也有离散量（变压器分接头挡位、电容器组数）。传统内点法、线性规划这类基于梯度的算法处理离散变量非常别扭，常见做法是先松弛成连续变量再四舍五入，但舍入结果经常不满足约束，还得手动校验。
• 状态变量（负荷节点电压、发电机无功出力）的越限需要通过罚函数或者约束处理机制反馈到目标函数里，反馈强弱直接影响优化方向。

我早期试过用内点法跑带变压器分接头的算例，每次求解到最后都要面临“四舍五入后潮流重新校验不通过”的问题，来回折腾非常痛苦。粒子群这类群体智能算法不需要任何梯度信息，离散变量直接通过取整/量化方式处理，天然贴合这种混合整数问题，这也是我后来把方案换成PSO的直接原因。

1.3 控制变量与状态变量的完整清单

建模之前必须先搞清楚手上有哪些“旋钮”，哪些是“观测表”。对于IEEE14节点系统，典型配置如下：

控制变量一共9维，这个维度对PSO来说非常友好。状态变量不参与优化搜索，而是在每个控制变量组合下通过潮流计算自动得到，所以算法的每次评估都要完整跑一遍潮流。

2. 粒子群算法原理：一群粒子如何协作搜索

2.1 位置-速度更新公式：PSO的直觉

粒子群算法的思路说起来很简单：把每个候选解看作搜索空间里的一个“粒子”，粒子有位置和速度两个属性。位置代表一组控制变量，速度代表当前搜索的方向和步长。每个粒子记住自己历史最好的位置pbest，整个种群共享一个全局最好的位置gbest，下一代的位置由这三个信息共同决定。

标准的更新公式是：

v(i+1) = w * v(i) + c1 * r1 * (pbest - x(i)) + c2 * r2 * (gbest - x(i)) x(i+1) = x(i) + v(i+1)

打个比方，想象一群鸟在陌生的山谷里找食物最密集的地方。每只鸟都知道自己飞过的最高食物点，也听得见群体里目前发现的最高点哨声，于是每次都朝这两个方向的合成方向飞一点，同时保留一部分自己原来的惯性。群体就这样慢慢收敛到食物最密集的区域。

在无功优化里，“食物”就是有功网损最低的控制变量组合。粒子每飞到一个新位置，就调用一次潮流计算去“尝”这个位置的网损高低，然后更新自己的记忆。这个过程完全不需要目标函数的梯度，这是它面对非线性潮流方程时最大的优势。

2.2 惯性权重、学习因子与边界处理

参数是PSO的命根子，我复现这个算例时固定用下面这套参数：

惯性权重线性递减是我最推荐的策略：迭代前期保持较大的w，粒子能在大范围内探索不同无功补偿组合；后期w变小，粒子速度降低，能在gbest附近细挖。固定w的话，要么前期收敛慢，要么后期震荡不收敛，我实测固定w=0.5后期经常出现适应度反复跳动的情况。

边界处理有两个关键点：位置越界时钳制到边界值，同时把对应速度置零，避免粒子反复“撞墙”；速度向量本身也要做限幅，用max(min(v, vmax), -vmax)包一层。这两个细节不做，粒子很容易飞到电压或变比的非法区间，潮流计算还会出现不收敛。

2.3 离散变量：变压器分接头与电容器组怎么编码

PSO更新出来的位置是连续实数，但变压器变比和电容组数是离散的。我的做法是在每次位置更新后做一次量化：

% 变压器变比按0.0125步长取整 x(6:8) = round(x(6:8) / 0.0125) * 0.0125; % 电容容量按5Mvar步长取整 x(9) = round(x(9) / 5) * 5;

注意量化之后位置可能越过边界，所以量化完还要再夹一次上下限。另外初始化粒子的时候也必须做同样的离散化，否则第0代粒子就会落在非法离散点上。遇到过有人把离散化只放在迭代循环里，初始化直接随机生成，结果前几代适应度忽高忽低，后面才慢慢追上来——这就是初始化没做量化导致的。

3. IEEE14节点系统的Matlab建模：从case14数据到可求解优化模型

3.1 为什么要用Matpower而不是手写牛顿-拉夫逊

手写潮流程序并不难，打开任何一本电力系统分析教材，按牛顿-拉夫逊法一步步列方程就行，但工程效率上完全没必要。Matpower是开源工具箱，自带case14标准数据文件，一个loadcase('case14')就能把母线、线路、发电机、变压器参数全部读进来。我要做的核心工作不是重写潮流，而是把粒子群输出的控制变量准确翻译成Matpower能识别的数据格式。

用Matpower还有一个好处：它内置了发电机无功越限处理和潮流收敛判断，省掉很多矩阵构建的细节。代价是它的数据组织方式是“矩阵+列索引”，修改控制变量时列号写错是最常见的低级错误，我下面会专门讲。

3.2 把控制变量“塞进”潮流数据的三个修改点

IEEE14节点里的5台发电机分布在节点1、2、3、6、8，其中节点1是平衡节点，其余是PV节点。粒子给出的前5维是机端电压，对应到Matpower的gen矩阵第6列（VG列），不是第2列（PG列）。我第一次写这段代码就误改过PG列，结果潮流怎么算都爆量，检查了很久才发现修改错位置。正确写法是用节点号反查gen行：

gen_bus = [1; 2; 3; 6; 8]; [~, idx_gen] = ismember(gen_bus, mpc.gen(:, 1)); mpc.gen(idx_gen, 6) = VG; % 第6列是机端电压

变压器变比在branch矩阵第9列（TAP列），但这里有个坑：不是所有支路都是变压器支路，不能用固定行号硬编。正确做法是先用find找出非零变比的支路：

tap_branch = find(mpc.branch(:, 9) ~= 0); % 确认一下是不是3条，且顺序和你的变量定义一致 mpc.branch(tap_branch, 9) = tap;

case14里变压器支路通常是4-7、4-9、5-6三条，但不同版本的Matpower在支路排列顺序上可能有差异，一定要先打印出来核对，再决定tap变量的对应顺序。

并联电容器我选择加在节点9，这是无功紧缺较明显的负荷节点。修改的是bus矩阵第5列（Bs列），注意单位换算和叠加问题：

bus9 = find(mpc.bus(:, 1) == 9); % 节点9原有并联导纳保留，叠加补偿容量（100MVA基准下转为标幺值） mpc.bus(bus9, 5) = mpc.bus(bus9, 5) + QC / 100;

很多初学者直接mpc.bus(bus9, 5) = QC / 100，等于把系统原有的固定补偿给覆盖掉了，优化出来的结果自然不对劲。叠加而不是覆盖，这个细节非常关键。

3.3 适应度函数：网损、罚函数与收敛判断

每次粒子位置确定后，目标函数分成三部分：

cost = Ploss + alpha * sum(电压越限平方和) + beta * sum(无功越限平方和)

有功网损从潮流结果里取：

res = runpf(mpc, mpoption('verbose', 0)); Ploss = sum(res.branch(:, 14) + res.branch(:, 16));

这个写法利用了branch矩阵第14列（PF，支路始端有功）和第16列（PT，支路末端有功），两者相加就是该支路的有功损耗。注意不是第18列，第18列在标准支路表里不存在。

电压罚项关注所有PQ负荷节点的电压越限，无功罚项关注每台发电机的出力越限。如果潮流不收敛，直接返回一个很大的惩罚值，比如1e6，把这个粒子的解直接判死。罚函数系数的量级选择我在第5节详细说，这里先记住alpha和beta不能拍脑袋乱取。

4. PSO无功优化的核心代码实现：主循环与适应度函数拆解

4.1 主循环骨架：初始化、迭代、记录

整体代码结构非常固定，我直接给出可以运行的骨架：

%% PSO主程序 clc; clear; close all; mpc = loadcase('case14'); % PSO参数 nPop = 40; % 种群规模 MaxIter = 100; % 最大迭代次数 w_max = 0.9; w_min = 0.4; c1 = 1.5; c2 = 1.5; % 控制变量维度与范围 % x = [VG1, VG2, VG3, VG6, VG8, tap1, tap2, tap3, QC] nVar = 9; lb = [0.95*ones(1,5), 0.90*ones(1,3), 0]; ub = [1.05*ones(1,5), 1.10*ones(1,3), 60]; vmax = 0.1 * (ub - lb); % 初始化 particle = repmat(lb, nPop, 1) + rand(nPop, nVar).*repmat(ub-lb, nPop, 1); velocity = vmax .* (2*rand(nPop, nVar) - 1); % 初始化时做一次离散化 for i = 1:nPop particle(i, 6:8) = round(particle(i, 6:8) / 0.0125) * 0.0125; particle(i, 9) = round(particle(i, 9) / 5) * 5; end pbest = particle; pbest_cost = inf(nPop, 1); cost = calFitness(particle(1, :), mpc); gbest = particle(1, :); gbest_cost = cost; % 迭代 record = zeros(MaxIter, 1); for iter = 1:MaxIter w = w_max - (w_max - w_min) * iter / MaxIter; for i = 1:nPop r1 = rand(1, nVar); r2 = rand(1, nVar); velocity(i, :) = w * velocity(i, :) ... + c1 * r1 .* (pbest(i, :) - particle(i, :)) ... + c2 * r2 .* (gbest - particle(i, :)); % 速度限幅 velocity(i, :) = max(min(velocity(i, :), vmax), -vmax); % 位置更新 particle(i, :) = particle(i, :) + velocity(i, :); % 边界钳制 particle(i, :) = max(particle(i, :), lb); particle(i, :) = min(particle(i, :), ub); % 离散化后再钳制一次 particle(i, 6:8) = round(particle(i, 6:8) / 0.0125) * 0.0125; particle(i, 9) = round(particle(i, 9) / 5) * 5; particle(i, :) = max(particle(i, :), lb); particle(i, :) = min(particle(i, :), ub); cost = calFitness(particle(i, :), mpc); if cost < pbest_cost(i) pbest_cost(i) = cost; pbest(i, :) = particle(i, :); end if cost < gbest_cost gbest_cost = cost; gbest = particle(i, :); end end record(iter) = gbest_cost; end

这段循环没有用并行，因为每个粒子都要修改mpc并调用runpf，串行反而省心。如果以后升级到IEEE39节点或者更大系统，可以考虑parfor替换最外层粒子评估，但要注意parfor里mpc和随机数流的独立性问题。

4.2 适应度函数代码：潮流求解＋惩罚项

function cost = calFitness(x, mpc) % 解析控制变量 VG = x(1:5); tap = x(6:8); QC = x(9); % 1) 设置发电机端电压 gen_bus = [1; 2; 3; 6; 8]; [~, idx_gen] = ismember(gen_bus, mpc.gen(:, 1)); mpc.gen(idx_gen, 6) = VG; % 2) 设置变压器变比 tap_branch = find(mpc.branch(:, 9) ~= 0); if length(tap_branch) ~= 3 error('变压器支路数量不是3条，请检查case14版本'); end mpc.branch(tap_branch, 9) = tap; % 3) 设置节点9并联电容（叠加，不覆盖） bus9 = find(mpc.bus(:, 1) == 9); mpc.bus(bus9, 5) = mpc.bus(bus9, 5) + QC / 100; % 4) 潮流计算 res = runpf(mpc, mpoption('verbose', 0, 'out.all', 0)); if ~res.converged cost = 1e6; return; end % 5) 有功网损 Ploss = sum(res.branch(:, 14) + res.branch(:, 16)); % 6) 负荷节点电压越限惩罚 pq_idx = find(res.bus(:, 2) == 1); V = res.bus(pq_idx, 8); Vmin = 0.95; Vmax = 1.05; Vpen = sum(max(0, V - Vmax).^2) + sum(max(0, Vmin - V).^2); % 7) 发电机无功越限惩罚 Qg = res.gen(:, 3); Qmax = res.gen(:, 4); Qmin = res.gen(:, 5); Qpen = sum(max(0, Qg - Qmax).^2) + sum(max(0, Qmin - Qg).^2); % 8) 加权目标 alpha = 1000; % 电压罚系数 beta = 5000; % 无功罚系数 cost = Ploss + alpha * Vpen + beta * Qpen; end

这里alpha和beta的单位需要理解一下：Ploss单位是MW，Vpen单位是p.u.的平方，Qpen单位是Mvar的平方，所以罚系数天然带着“把越限代价放大到和网损可比”的含义。数值太小，粒子会为了压低网损在非法电压区间游荡；数值太大，网损项完全被淹掉，算法退化成“只找可行解”的搜索器。

4.3 可视化与结果导出

迭代结束后，我习惯做三件事：画收敛曲线、打印最优控制变量、对比优化前后电压分布。

figure; plot(record, 'LineWidth', 2); xlabel('迭代次数'); ylabel('最优适应度'); title('PSO收敛曲线'); grid on; fprintf('最优控制变量:\n'); fprintf('发电机端电压(p.u.): %.3f %.3f %.3f %.3f %.3f\n', gbest(1:5)); fprintf('变压器变比: %.4f %.4f %.4f\n', gbest(6:8)); fprintf('补偿容量(Mvar): %.1f\n', gbest(9));

收敛曲线是判断算法是否稳定的第一张图：如果曲线一直台阶式跳变，说明vmax或罚系数有问题；如果曲线早早平坦但数值较高，可能陷入局部最优，需要回归参数或者重新初始化。

5. 复现中反复踩到的坑：从变压器变比到罚函数系数

5.1 Matpower的PV→PQ切换会把你的罚函数“架空”

我在第3节写的适应度函数里加入了发电机无功越限罚项，但实际跑起来你会发现，绝大多数粒子的Qpen都是0。原因不在罚函数写错，而在于Matpower的runpf内部会自动处理发电机无功越限：某个PV节点的无功出力一旦超过Qmax或低于Qmin，该节点会被自动改判为PQ节点重新求解潮流，结果里返回的Qg就被“修正”到限制值附近。

这带来一个隐蔽问题：看似约束满足，实际系统状态已经变了，网损也会随之变化。我的处理建议是保留Qpen罚项，同时额外检查潮流结果中是否有节点被转换，如果res.gen里某台发电机的Qg恰好顶在Qmax上，就要留个心眼，说明当前解已经把该发电机无功推到极限，继续增加补偿可能才是正确方向。对比手动实现牛顿-拉夫逊的程序，Matpower这个自动转换机制常常让新手误以为“罚函数失效”，其实它是帮你在潮流层面兜底了。

5.2 变压器支路顺序与变比赋值错位

前面代码里用find自动找变压器支路，这个方法能解决行号问题，但有一个新风险：很多版本的case14里，三条变压器支路的顺序并不是按4-7、4-9、5-6排列的。如果你的tap向量依次对应的是“4-7、4-9、5-6”，而find返回的是“5-6、4-7、4-9”，那么赋值后变比全部错位，优化结果毫无意义。

我每次写新数据文件时都会先跑一句：

[match, idx] = ismember([4 7; 4 9; 5 6], mpc.branch(tap_branch, 1:2), 'rows');

确认顺序后再赋值。这种“校验一遍数据顺序”的习惯，能省掉后面至少一个小时的排查时间。

5.3 罚函数系数取多少：从量纲角度理解

罚系数不是玄学，可以从量纲反推。IEEE14节点初始网损约13MW，电压越限0.02p.u.时平方项是0.0004，如果alpha取1000，罚项贡献是0.4，相当于网损的3%，这个量级正好能推动粒子避开越限同时不吞掉优化信号。无功越限罚项的量级也类似，beta取5000时，无功越限2Mvar贡献20kW当量，对总成本的影响比较明显。

我实测下来，电压罚系数alpha在1000~10000之间都能得到可行解，beta在5000~20000之间比较稳。如果发现优化结果里电压违规点很多，就先加大alpha；如果发现曲线非常平但是数值偏高，就适当下调beta，让网损项重新成为主导。这个调参循环通常两三轮就能收敛到稳定区间。

6. 仿真结果分析：网损、电压与收敛曲线到底说明什么

6.1 有功网损：优化前后的账

我按上面的参数跑了一组实验，初始潮流下系统有功网损约13.4MW，优化后降到11.3MW左右，降幅约15.6%。这个数字在同类文献里属于正常区间，具体数值会因为case14版本、补偿节点选择、罚系数和随机种子不同而有差异，但趋势是一致的：PSO能找到明显优于初始状态的无功组合。

有一个现象值得留意：最优解往往不是把所有补偿都投满，而是把节点9的电容器组投到15~20Mvar左右，配合变压器变比调整。因为过补偿会导致无功倒送，抬升电压的同时反而增大损耗，这正是“无功优化”区别于“无功补偿越多越好”的直观证明。

6.2 负荷节点电压：越界点如何被“拉回”

初始状态下，节点14电压约0.98p.u.，虽然没到0.95的下限，但裕度已经很小。优化后节点14电压提升到1.01左右，系统内所有PQ节点的电压基本上都落在0.98~1.05区间，电压分布整体上移且更均匀。

这说明罚函数起了作用——粒子在搜索过程中并不是只顾压低网损，一旦某个粒子试图通过降低末端电压来获得更低的网损，电压越限罚项就会把它的适应度拉高，让它被淘汰。最终保留的gbest是网损和电压质量两个目标之间的平衡点。

6.3 收敛曲线与多次运行的稳定性经验

我保存的收敛曲线显示，gbest在前30代快速下降，从初始的15左右一路降到约12，之后进入约50代的小幅优化，最终稳定在11.3附近。中间偶尔会出现小幅跳变，那是粒子在探索新的离散组合时短暂找到更优解的标记。

PSO是随机算法，单次运行不能代表全部。我连续跑10次的记录里，gbest波动范围大约在11.0~11.6MW之间，相对稳定。如果你在复现时发现每次都差很多，多半是vmax或者惯性权重设置不合理，导致某些运行过早收敛到局部最优。稳妥做法是固定一个随机种子，或者取10次运行中的最优结果作为最终结论——这也是论文里常见的处理方式。

如果后续想把这套流程迁移到IEEE30、IEEE39甚至实际电网模型，只需要换loadcase的数据文件，重新梳理控制变量范围和初始shunt，主循环和适应度函数的骨架基本不用动。我在迁移过程中最大的体会是：粒子群算法在中等规模系统上完全够用，真正决定项目成败的往往不是算法本身，而是数据接口、离散化时机和罚函数量级这些“配角”细节。