前阵子搭了一个优化演示工程，帮一个朋友把 IEEE 30 节点系统的网损降下来。算法选的粒子群，底层潮流计算交给 Matpower 引擎，整套跑在 Matlab 里。借这个机会，我把整个项目的建模思路、代码骨架、以及调试过程中踩过的坑完整记录下来。严格来说，IEEE 30 节点系统是一个输电网测试算例，但在实际工作中，很多做配电网算法验证的同学也拿它当作标准网架来用，所以这个标题组合非常经典。文章不空谈原理，重点讲清楚我怎么建模、怎么把粒子群和 Matpower 串起来、以及最终如何分析和判断结果，希望能给正在做电力系统智能优化方向的同学一个直接可参考的模板。
1. 这个优化问题到底在算什么：有功-无功联合优化的数学模型与工程动因
1.1 为什么要把"有功"和"无功"放在同一个问题里
很多刚开始接触电网优化的同学会有一个困惑：有功功率和无功功率不是两回事吗？有功管系统频率，无功管电压水平，为什么非要放在一起优化？
其实这是从运行机理上就耦合在一起的。电网里任何一条线路、一台变压器，它的损耗都可以写成 (P_{loss} = I^2 R) 的形式，而电流 I 的幅值同时由有功潮流和无功潮流共同决定。也就是说，线路上的无功流动同样会带来发热损耗。如果你只调整有功出力，不管无功分布，那么某些线路上的无功环流会一直存在，网损降不下去；反过来，只调无功，有功的不合理分配又会让重载线路继续超负荷。所以真正要降网损，必须同时调整有功出力和无功分布，这就是"有功-无功联合优化"的由来。
配电网的场景里这个问题更明显。配电网的线路 R/X 比值大，无功流动造成的电压跌落和线损占比比输电网高得多。分布式电源接入之后，有功和无功的耦合关系又进一步复杂化。所以在配电网规划、降损、电压治理这些课题里，有功-无功联合优化基本上是一个绕不开的基础问题。
1.2 目标函数与决策变量怎么选
这次工程里，目标函数取的是系统总网损最小。原因很朴素：网损是可观测、可计算的直接经济指标，也是最好和基态潮流做对比的结果量。你要是想换成发电成本最小、电压偏差最小，方法完全一样，只是改动适应度函数里的一个表达式而已。
决策变量的选取是建模里最关键的一步。IEEE 30 节点系统一共有 30 个节点、41 条支路，6 台发电机。我采用的是最常见的变量组合方式：
- 除平衡节点外的 5 台发电机有功出力；
- 6 台发电机的机端电压幅值。
也就是说，一个粒子是一个 11 维向量，里面的每个分量被提交给 Matpower 潮流计算，最终得到一组潮流结果和一个网损值。为什么用机端电压代替无功出力？因为这两个量在潮流计算里是强耦合的，你设定了机端电压，潮流解出来后自然会得到对应的无功注入。直接拿无功出力当决策变量反而麻烦，因为 Matpower 在计算时会有 PV 节点和 PQ 节点之间的自动切换逻辑，直接用无功出力很容易撞上这个切换边界。这一点后面踩坑部分我会单独说。
1.3 约束条件从哪来：潮流等式与设备不等式
优化问题不能只有目标函数，约束条件才是让它"像电网"而不是"像数学游戏"的关键。约束分两类：
等式约束就是潮流方程本身，配电网或者输电网的任何运行状态都必须满足节点功率平衡。这个我不会自己手写，而是委托给 Matpower 的runpf函数来完成，跑通了就代表等式约束被满足了。
不等式约束包括：发电机有功上下限、发电机无功上下限、节点电压幅值上下限、线路传输容量上限。这些限制在 Matpower 的case30数据文件里都有初始值。bus表里有每个节点的正常电压上限Vmax和下限Vmin，gen表里有每台发电机的Pmax、Pmin、Qmax、Qmin。在粒子群算法框架下，不等式约束大部分靠罚函数处理，小部分靠边界截断处理，这个设计展开讲就是第 4 章的内容。
让我点一下数学上的关键逻辑：目标函数要优化的是网损，决策变量是发电机有功和机端电压，约束是潮流方程加设备上下限。把这三句话写清楚，这个优化问题才算是被正确建模了。很多同学代码写不明白，往往不是粒子群没掌握，而是决策变量和约束没有从 Matpower 数据结构上想透。
2. 为什么用粒子群而不是直接调用内点法
2.1 Matpower自带最优潮流的局限
Matpower 自带的runopf函数是用内点法求解最优潮流（OPF）的，精度高、速度快，在干净的标准算例上效果相当好。那为什么还要自己做粒子群？
一个很现实的原因是，实际工程项目里的优化模型往往和教科书上的标准 OPF 不完全一样。你可能要加一个不光滑的约束（比如某个节点的电压合格率分段惩罚），可能要把目标函数改成某种启发式指标，可能要对某些设备动作次数做惩罚。这些非标准、非光滑、甚至不可导的改动，对传统内点法是致命的，因为内点法需要构造拉格朗日函数并求梯度或者海森矩阵，一旦模型变得不光滑，数值上非常容易发散。
粒子群这类群智能算法则是另一套逻辑：它不要求目标函数可导，不要求约束光滑，它只把目标函数当成一个黑盒来调用。潮流算得出来，就给一个具体的网损值；算不出来，就给一个巨大的惩罚值。这种"黑盒优化"的特性，让它在科研探索阶段显得特别灵活。
2.2 粒子群适合这种问题的真正原因
工程上选择粒子群，主要有三个原因。
第一，实现成本低。一个有速度-位置更新公式加适应度计算的粒子群核心代码，写得好也就一两百行。相比于内点法复杂的矩阵求导和修正方程，粒子群的代码量要小得多，调起来也更直观。
第二，天然支持并行。粒子群每一代要计算多个粒子的适应度，而这些计算彼此独立，完全可以用parfor并行跑。在我这个例子里，每个粒子都要调用一次runpf，用并行池可以把整体耗时压缩到原来的三分之一左右。
第三，全局搜索能力比传统梯度类算法好。电网优化问题本质上是有多个局部最优点的非线性非凸问题，内点法非常依赖初值，初值不好就容易收敛到局部解或者压根不收敛。粒子群通过种群协作的方式在解空间里撒点搜索，虽然不保证找到全局最优，但工程实践里找到的解往往是可用的，而且稳定性好。
2.3 它与遗传算法的简单对比
也有人问我，为什么不直接上遗传算法？两个算法我都实现过，结论是：在"电网潮流作为黑盒"这个场景下，粒子群通常更省事。
遗传算法的选择、交叉、变异三个算子，每个都需要调试合理的概率参数，否则要么收敛太慢，要么种群过早失去多样性。粒子群的参数相对集中，核心就是惯性权重和学习因子，尤其惯性权重可以按迭代次数线性递减，整体调参路径清晰得多。另外，遗传算法的二进制编码在边界连续优化问题里会有精度和编码长度之间的矛盾，而粒子群直接用实数向量操作，不需要解码，天然贴合"发电机出力用连续实数表示"的需求。
当然，如果问题本身是离散变量主导的（比如电容器分组投切、变压器分接头挡位），遗传算法的离散处理反而更方便。连续变量主导的有功-无功优化，粒子群是性价比最高的选择。
3. 环境准备与数据摸底：Matpower、case30和配电网改造
3.1 十分钟搞定Matpower环境
Matpower 是一个纯 Matlab 实现的开源电力系统分析工具包，不需要编译，下载解压就能用。从官网或者 GitHub 下载最新稳定版本之后，把解压目录整个加入 Matlab 路径即可。
%% 将Matpower路径加入Matlab环境 addpath(genpath('D:/matpower7.1')); savepath; % 保存路径，避免每次启动重新设置 mpver('ver'); % 查看并确认Matpower版本运行mpver('ver')能看到默认的case30等标准算例已经内置。我建议每次用 Matpower 之前都显式确认版本号，因为不同版本对runpf返回结果的字段略有差异，确认版本可以少很多莫名其妙的"数组越界"报错。
3.2 读懂case30数据结构：bus、gen、branch
Matpower 的核心数据结构是mpc结构体，它包含bus、gen、branch三个大矩阵，以及baseMVA这个基准容量。case30.m文件本质上就是这些矩阵的文本初始化。
bus矩阵的关键列：
- 第 1 列节点编号；
- 第 2 列节点类型：1 是 PQ 节点，2 是 PV 节点，3 是平衡节点；
- 第 3、4 列是有功负荷和无功负荷；
- 第 8 列是节点电压幅值初始值；
- 第 12、13 列是正常运行时电压上限和下限。
gen矩阵的关键列：
- 第 1 列是发电机所在节点号；
- 第 2 列是有功出力；
- 第 3 列是无功出力；
- 第 4、5 列是无功上限和下限；
- 第 6 列是机端电压幅值设定值；
- 第 9、10 列是有功上限和下限。
branch矩阵记录线路或者变压器的参数，包括首末端节点、电阻、电抗、电纳、长期允许载流量、变压器变比等。
第一件事永远是用loadcase('case30')载入数据，然后用runpf跑一次基态潮流，看看初始网损是多少。基态潮流的结果是整个优化项目的比较基准，后面所有优化效果的评估都靠它说话。
mpc = loadcase('case30'); % 载入IEEE30节点数据 res0 = runpf(mpc); % 基态潮流 loss0 = sum(res0.branch(:, 14) + res0.branch(:, 16)); % 基态网损，单位MW fprintf('基态网损: %.4f MW\n', loss0);res.branch矩阵的第 14 列是首端有功注入，第 16 列是末端有功注入，两者相加就是该支路的损耗。Matpower 返回的success字段等于 1 才代表潮流收敛，这个字段在后续迭代里非常重要。
3.3 再说一遍：这是输电网算例，怎么改成配电网风格
我必须把话说清楚：IEEE 30 节点系统本身是输电网络，线路阻抗、电压等级、负荷水平都不是典型配电网。但是在算法开发和验证阶段，它的拓扑结构、发电机数量、约束类型足够丰富，完全可以通过改造来模拟配电网特征。
常见的改造思路有三种：
- 把部分发电机改成分布式电源模型，降低其出力上限，并修改无功调节范围；
- 将线路参数调大电阻分量，提高 R/X 比值，让电压问题更突出；
- 增加部分节点的无功负荷，让无功分布更加紧张。
这样改完之后，还是在case30的框架下做优化，但结果更接近配电网降损的真实需求。如果你已经有某个具体的配电网拓扑数据，也可以用自己的潮流数据替换掉case30，粒子群的代码框架完全不用动，只需要改loadcase那一个地方。
4. 粒子群与潮流计算耦合：编码、适应度函数、主循环
4.1 粒子编码设计：一串向量对应一台发电机的运行方式
粒子群优化的第一个核心步骤是编码。我的设计是：把每台非平衡发电机的有功出力，和全部发电机的机端电压幅值，拼接成一维向量。
以case30为例，6 台发电机节点编号分别是 1、2、5、8、11、13，其中节点 1 是平衡节点。那么决策向量定义为：
- 前 5 个分量：节点 2、5、8、11、13 这 5 台发电机的有功出力；
- 后 6 个分量：节点 1、2、5、8、11、13 这 6 台发电机的机端电压幅值。
总维度是 11。这里有一个细节，有些代码会把平衡节点的有功也拿来做变量，这是不对的。平衡节点的有功在潮流里是"补差"角色，系统总有功不平衡时由它自动平衡，把它当决策变量会和潮流方程矛盾。
为了让粒子初始化范围合理，我直接从mpc.gen里读取Pmax、Pmin、Vmax、Vmin来生成上下界。这样粒子群搜索空间和物理设备允许范围自动对齐，省去手动维护边界数据的麻烦。
% 读取设备上下限 genBus = mpc.gen(:, 1); isSlack = zeros(size(genBus)); % 平衡节点由bus矩阵type列为3的节点确定 slackIdx = find(mpc.bus(:, 2) == 3); isSlack = ismember(genBus, slackIdx); % 有功决策量范围: 非平衡发电机 PgMin = mpc.gen(~isSlack, 10); PgMax = mpc.gen(~isSlack, 9); % 电压决策量范围: 全部发电机 VgMin = 0.94; VgMax = 1.06;4.2 适应度函数：网损加罚函数，一个值得反复调试的核心
适应度函数是粒子群和电网模型之间的唯一接口。我的设计分三层。
第一层：把粒子的 11 个分量写回mpc.gen，调用runpf计算潮流。如果潮流不收敛，直接返回一个很大的惩罚值，比如1e6。这一步一定要包在try-catch里，因为粒子群探索阶段经常会产生不合理的控制变量组合，让潮流计算直接报错退出，如果不捕获异常，整个优化循环就崩了。
第二层：计算网损。在潮流收敛的情况下，用首末端有功之和得到系统网损。
第三层：罚函数处理越限。我主要检查两个地方：节点电压幅值是否越界，以及发电机无功出力是否越限。节点电压在res.bus第 8 列，发电机无功在res.gen第 3 列。越限越严重，惩罚值越大，这样粒子群在迭代过程中会自然避开不合理区域。
function loss = psoFitness(x, mpcBase) mpc = mpcBase; genBus = mpc.gen(:, 1); slackIdx = find(mpc.bus(:, 2) == 3); isSlack = ismember(genBus, slackIdx); % 前5维写到非平衡机有功 mpc.gen(~isSlack, 2) = x(1:5); % 后6维写到全部机端电压 mpc.gen(:, 6) = x(6:11); try res = runpf(mpc); if res.success ~= 1 loss = 1e6; return; end % 网损 loss = sum(res.branch(:, 14) + res.branch(:, 16)); % 电压越限惩罚 vmin = min(res.bus(:, 8)); vmax = max(res.bus(:, 8)); if vmin < 0.95 loss = loss + 2000 * abs(vmin - 0.95); end if vmax > 1.05 loss = loss + 2000 * abs(vmax - 1.05); end % 发电机无功越限惩罚 if min(res.gen(:, 3) - res.gen(:, 4)) < 0 || min(res.gen(:, 5) - res.gen(:, 3)) < 0 loss = loss + 3000; end catch loss = 1e6; end end这个适应度函数里，罚系数 2000、3000、1e6 都是我反复试验后比较稳的组合。罚函数设计是这个项目里最需要耐心调试的部分，权重太小约束形同虚设，权重太大粒子群会过早向可行域收缩，丧失搜索能力。后面我会单独讲怎么判断罚函数是否合适。
4.3 速度位移更新与边界处理
粒子群的标准更新公式不用担心，直接用经典版本：
v = w * v + c1 * r1 * (pbest - x) + c2 * r2 * (gbest - x) x = x + v其中w是惯性权重，c1是自我认知学习因子，c2是社会学习因子，r1、r2是 0 到 1 的随机数。我的参数设定如下表：
| 参数 | 取值 | 说明 |
|---|---|---|
| 惯性权重 w | 0.9 到 0.4 线性递减 | 前期全局搜索，后期局部精细搜索 |
| 学习因子 c1 | 2.0 | 保留自我经验 |
| 学习因子 c2 | 2.0 | 群体信息共享 |
| 种群规模 | 40 | 网络规模不大，40 个粒子足够 |
| 最大迭代次数 | 200 | 配合收敛曲线判断是否提前结束 |
| 速度上限 vmax | 0.3 | 每个决策量按范围归一化后限制步长 |
线性递减的惯性权重是工程实践里非常稳定的配置。迭代初期w大，粒子跑得快，可以探索大范围；迭代后期w小，粒子在局部精细搜索。很多失败的粒子群应用案例，问题都出在全程用同一个w，导致搜索行为要么太粗放要么太保守。
边界处理上我的做法是：粒子越界时，直接把对应分量拉回到边界值，并且把该维速度设为 0。这个"重置边界"策略比"随机反弹"策略更容易让粒子在可行域边界附近稳定下来，因为在电力系统优化里，很多最优解往往就在发电机出力上限附近。
4.4 主循环骨架与计算量预估
主循环的骨架大致是这样：
% 初始化粒子群 nPop = 40; nVar = 11; maxIter = 200; x = repmat(xMin, nPop, 1) + rand(nPop, nVar) .* repmat(xMax - xMin, nPop, 1); v = zeros(nPop, nVar); pbest = x; pbestVal = arrayfun(@(i) psoFitness(x(i, :), mpc), 1:nPop); gbestIdx = find(pbestVal == min(pbestVal), 1); gbest = pbest(gbestIdx, :); gbestVal = pbestVal(gbestIdx); trace = zeros(maxIter, 1); % 主迭代 for iter = 1:maxIter w = 0.9 - 0.5 * (iter / maxIter); % 线性递减 for i = 1:nPop r1 = rand(nVar, 1); r2 = rand(nVar, 1); v(i, :) = w * v(i, :) + 2.0 * r1' .* (pbest(i, :) - x(i, :)) + 2.0 * r2' .* (gbest - x(i, :)); v(i, v(i, :) > vmax) = vmax; v(i, v(i, :) < -vmax) = -vmax; x(i, :) = x(i, :) + v(i, :); x(i, x(i, :) > xMax) = xMax; x(i, x(i, :) < xMin) = xMin; fval = psoFitness(x(i, :), mpc); if fval < pbestVal(i) pbest(i, :) = x(i, :); pbestVal(i) = fval; end end [minVal, idx] = min(pbestVal); if minVal < gbestVal gbestVal = minVal; gbest = pbest(idx, :); end trace(iter) = gbestVal; end计算量方面，40 个粒子乘以 200 代，就是 8000 次潮流计算。在普通笔记本上，一次case30潮流大约几十毫秒，所以整个优化过程大概一到三分钟。如果需要跑大规模配电网，建议改用parfor对粒子循环做并行化，实测能缩短到原来的四分之一左右。
5. 跑通之后怎么分析结果：迭代曲线、电压分布、网损对比
5.1 从基态潮流开始：先记住优化前的网损
不管优化结果多漂亮，没有基准值就毫无意义。我的习惯是把基态潮流结果、基态网损、基态最低最高电压全部打印出来，作为后续对比的 Baseline。
基态网损由runpf直接得到，然后优化结束之后，把gbest重新写入mpc.gen，再跑一次潮流，同样统计网损和电压。两者做差除基态，就是降损率。
这里有一个细节，优化停止时gbest对应的潮流并不一定已经用最终粒子再算过一遍，因为pbestVal可能是粒子在某一代的位置计算出来的。所以最终结果必须先重算潮流，再统计指标，否则容易出现"结果和实际不一致"的尴尬。
5.2 迭代曲线与收敛判断
粒子群迭代过程里，我把每一代的最优适应度记录在trace数组里，用semilogy画收敛曲线。注意这里我通常用对数坐标，因为早期适应度可能从几千快速下降到十几个 MW，线性坐标会让后期曲线贴在横轴上，看不出细节。
一个正常的收敛曲线应该是：前期快速下降，中期平缓波动，后期基本稳定。如果曲线后期还在明显持续下降，说明迭代次数不够，要加大maxIter。如果曲线很快就完全不动了，说明种群可能早熟，所有粒子已经挤到了同一个局部最优点。
figure; plot(1:maxIter, trace, 'LineWidth', 1.5); xlabel('迭代次数'); ylabel('最优网损/MW'); title('粒子群收敛曲线'); grid on;5.3 用runopf兜底，验证粒子群结果是否可信
做完粒子群优化，我会顺手用 Matpower 自带的runopf再解一次同样的问题。这个对照有双重价值。
第一，runopf的内点法收敛精度高，它的网损结果可以作为理论参考值。粒子群如果能跑到和内点法接近的网损，说明算法实现没有问题；如果差得远，大概率是编码范围或罚函数出了问题。
第二，runopf返回的最优运行点也可以拿来和粒子群最优解做交叉验证。把两组解分别回写到mpc.gen里跑潮流，对比它们的约束满足情况。如果粒子群给的解在电压约束、无功约束上全部满足，而网损又接近内点法，那这个优化结果就可以放心落地方案了。
这里要提醒一点，runopf是确定性求解，粒子群是随机求解，两者结果有差异很正常，不必追求完全一致。粒子群的工程意义在于处理内点法搞不定的非标准模型，而不是替代它。
6. 实操中绕不开的坑：随机种子、无功越限与潮流不收敛
6.1 换一台电脑结果就不一样？随机性管理
粒子群是随机算法，第一次跑出网损 10.5MW，第二次跑出 11.8MW，这是正常的。但问题是你不希望"提交结果时还能复现当时的实验"。我刚开始跑的时候吃过这个亏：换一个随机种子，结果就变，评审时完全答不上来自己的数据是怎么来的。
解决办法很简单，优化循环开始前固定随机流：
rng(7);固定种子之后，同一套参数每次跑的结果完全一致，可以复现。我还会把最终结果连同rng种子的状态一并保存：
save(sprintf('pso_case30_%s.mat', datestr(now, 'yyyymmdd_HHMM')), 'gbest', 'gbestVal', 'trace', 'loss0');另外，不要只跑一次就信数据。工程实践里我会用 3 到 5 个不同种子各跑一遍，取最优结果。这样既保证可复现性，又保留随机算法的全局搜索优势。
6.2 发电机无功越限与PV-PQ切换问题
这是我在这个项目里踩过最深的一个坑。粒子群把机端电压调得很高时，发电机可能需要发出很大的无功，但每台发电机都有Qmax限制。Matpower 在潮流计算过程中，如果发电机无功达到上限，会自动把该节点从 PV 节点切换成 PQ 节点，并让无功固定在边界值上。这意味着你设定的"机端电压"实际上并没有达到，优化结果的物理含义已经变了。
最直接的判断方法是：优化结束后，检查res.gen里的无功出力是不是落在Qmin到Qmax区间内。一旦有越限，要么罚函数起作用，要么说明电压上限设得太宽松，应该把电压决策范围从 0.94 到 1.06 收紧为 0.95 到 1.05。
我在适应度函数里对无功越限做了 3000 的固定惩罚，效果好于连续罚函数。因为无功越限是"硬性不合理"，固定大额惩罚能让种群快速淘汰这种个体。
6.3 当粒子让潮流算不出来时怎么办
粒子群里总有一些个体特别激进，比如让某台发电机出力特别小，而某个大负荷节点的功率缺额全部压在一条重载线路上，潮流可能直接发散。遇到这种情况，代码会抛异常，如果不处理，整个for循环就断了，迭代矩阵还少一行，最后程序崩得非常难看。
我的方案就是前面代码里展示的try-catch加success判断。返回一个超大惩罚值1e6，让这个粒子在竞争中自然淘汰。这个操作可以说是粒子群和潮流计算耦合项目里最重要的工程习惯。没有它，你会花大量时间在和粒子群"斗智斗勇"上，而不是真正关注优化结果。
6.4 罚函数权重与早熟问题
罚函数权重过大，粒子会很快放弃探索那些"暂时越限但可能引导到好解"的区域，整个种群迅速收缩到一小块安全区域，陷入局部最优。我遇到过一次把电压罚函数权重设为 5000 的情况，结果迭代曲线在 20 代就平稳了，最终网损反而比基态只降了几个百分点。把权重降到 2000 之后，收敛曲线明显更健康，迭代到 150 代才稳定，网损也降得更低。
罚函数调参的经验是：先放宽，后收紧。第一轮优化尽量让粒子群多探索，把罚函数权重设小一点；第二轮再把关键约束的权重加上去，在可行域范围内精细搜索。这个"由松到紧"的策略，可以显著降低早熟概率。
6.5 把优化结果回写到MPC再做闭式校验
最后一步我强烈建议做：拿优化出的gbest，重新构造一份mpc，然后打印完整的潮流结果，包括每个节点的电压幅值、每台发电机的实际出力和无功、每条支路的负载率。不要只看一个网损数字就结束。这份闭式校验数据才是你写报告、写论文、向别人说明优化效果时真正有价值的材料。
mpc_opt = mpcbase; genBus = mpc_opt.gen(:, 1); slackIdx = find(mpc_opt.bus(:, 2) == 3); isSlack = ismember(genBus, slackIdx); mpc_opt.gen(~isSlack, 2) = gbest(1:5); mpc_opt.gen(:, 6) = gbest(6:11); res_opt = runpf(mpc_opt); loss_opt = sum(res_opt.branch(:, 14) + res_opt.branch(:, 16)); fprintf('优化后网损: %.4f MW, 降损率: %.2f%%\n', loss_opt, (loss0 - loss_opt) / loss0 * 100);这个习惯让我绕过了很多次"看起来优化了，实际上潮流根本不收敛"的尴尬，建议你也养成。
最后再分享一点我个人的工作习惯：每次调完参数，我会顺手把随机种子、参数配置、最终网损和迭代曲线全部存成一个.mat文件，文件名里带上日期和参数关键字。隔几周再回头看，仍然能准确还原当时做了什么事情。这不是什么高深技巧，但对做优化仿真的人来说，绝对能节省大量返工时间。
