电气互联系统有功-无功协同优化Matlab实现与调试实战

最近在做“碳中和”背景下的电气互联系统优化项目时，我翻遍了手头的资料，发现大部分讨论都停留在纯电网的有功经济调度，或者单独的无功电压控制，真正把“有功-无功协同优化”和“电气互联”这两个维度揉在一起、还给出完整Matlab代码实现的资料，真的不多。这个项目前前后后跑了快两个月，从最初的模型设计到代码落地，中间踩了不少坑，也积累了一些很实用的调试经验。这篇博文就把我搭建电气互联系统有功-无功协同优化模型的完整过程整理出来，从数学模型怎么搭、目标函数怎么选、约束怎么处理，到Matlab代码怎么组织、求解器怎么调，再到我实际遇到的问题和排查思路，一次性讲清楚。

这篇文章适合正在做电力系统优化调度方向的研究生、从事电网运行方式计算工作的工程师，以及准备用Matlab复现相关算法的初学者。如果你对无功优化、新能源消纳、多能互补系统建模感兴趣，这篇应该能帮你省下大量摸索时间。

1. 为什么“碳中和”让无功优化变得不再简单

1.1 新能源高渗透带来的新变量

以前做传统电网的无功优化，情况相对简单：火电、水电等同步机型机组数量可控，各节点电压特性和无功出力范围已经通过多年的运行经验摸得很透，目标就是在线路损耗和电压质量之间找一个平衡点。但“碳中和”目标下，风电、光伏大规模接入，电力系统的运行工况发生了根本性变化。

新能源机组出力有天然的随机性和波动性。光伏夜间出力归零，风电一阵一阵地变化，这些都会直接改变电网各节点的有功注入，进而影响无功潮流分布和节点电压水平。更麻烦的是，逆变器型电源并没有传统同步发电机那样的转动惯量和励磁调节能力，它的无功支撑能力受逆变器容量限制，且控制策略不同，呈现出来的无功外特性完全不一样。分布式光伏大量接入配电网之后，原本只在输电网层面考虑的电压越限问题，现在已经开始出现在10kV甚至380V台区，电压控制的时间尺度和空间尺度都大大延展了。

这还不算完。电气互联系统的概念在于，电网和天然气网、热力网之间出现了越来越多的耦合环节。燃气轮机既可以发电又能供热，电转气（P2G）设备能把多余的可再生能源电力转化成天然气储存起来，天然气网的压缩机也需要从电网取电。这些跨网耦合设备让有功和无功的联系从过去的“松耦合”变为“强耦合”，单纯用传统的解耦思路去优化，已经很难满足安全经济运行的要求。

1.2 有功和无功为什么必须放在一起看

传统的优化调度通常分两步走：先做有功调度，确定各机组的有功出力计划；再做无功优化，在已知有功计划的基础上调整无功源，改善电压质量和降低网损。这种两级优化的逻辑在系统运行工况相对稳定的年代没有太大问题，因为有功计划基本确定，无功优化的跷跷板空间本来就有限。

但新能源大规模接入后，这种“先有功、后无功”的解耦策略开始暴露出明显的弊端。首先是时间尺度的错位，新能源有功出力的快速波动意味着每隔十五分钟甚至五分钟，系统有功平衡状态就发生一次变化，而无功优化的结果往往基于之前的顶峰潮流断面，用在新的运行点上很可能已经不再最优。其次是安全风险的增加，如果只考虑有功而忽略无功电压的约束，可能出现优化结果在角度上“经济最优”，但实际上节点电压已经越限，根本不可执行。反过来，只做无功优化时又必须把有功出力看作固定量，没法通过调节有功来缓解电压问题。

举个我在实际算例里遇到的情况。某一时段风电大发，输电网末端节点电压偏高，如果只做无功优化，需要投大量的电抗器或让附近机组大量吸收无功，但是参与调节的机组可能已经顶到无功出力下限了，无功优化怎么算都找不到可行解。而如果把有功和无功放在一起协同优化，适当的降低受端附近风电场的出力（当然这是最后手段），让潮流分布更均匀，电压问题就能在不增加无功补偿容量的前提下得到缓解。这种“用有功配合解决无功问题”的思路，就是有功-无功协同优化的直观价值所在。

2. 协同优化模型怎么搭：目标函数、变量与约束的完整设计

2.1 目标函数怎么选：我为什么用多目标加权

有功-无功协同优化的目标函数选择，直接决定了整个优化的走向。我第一版模型只把系统总有功网损作为目标，跑出来的结果确实网损下来了，但是机组二氧化碳排放量反而上升了——因为为了降网损，优化算法倾向于让更多本地机组发电而减少远方大机组的出力，导致低效率高排放的小机组承担了更多负荷。这和“碳中和”的初衷完全相反。

于是我把目标函数改成多目标加权形式：

[ \min F = w_1 \cdot C_{loss} + w_2 \cdot C_{carbon} + w_3 \cdot C_{curtail} ]

其中：

• (C_{loss}) 是系统总有功网损折算后的经济价值，单位是元/小时；
• (C_{carbon}) 是系统运行产生的碳排放罚金，通过各发电单元的有功出力和单位电量碳排放因子相乘后累加得到；
• (C_{curtail}) 是弃风弃光的惩罚项，用弃电量和对应上网电价的乘积表达；
• (w_1)、(w_2)、(w_3) 是权重系数，我取的是 0.5、0.3、0.2，但需要根据不同算例重新标定。

这里有个细节值得说。通常碳排放的计算要区分能源类型，燃煤机组的排放因子要比燃气轮机高出一倍以上，而风电、光伏基本不排碳。在Matlab里我用一个列向量 co2Factor 对应每台发电设备的单位电量排放量，目标函数里直接做向量点乘，这样程序写起来很简洁，后续改排放系数也方便。需要注意不同量纲的归一化问题，网损可能只有几十兆瓦，但碳排放罚金可能是几百甚至上千的数值，不归一化的话权重系数基本失去意义。我采用的方法是先做一次基础潮流计算，用基准工况下的网损值和碳排放值分别做缩放基准，把三个分目标都归一化到0到1的量级，再乘以权重。

2.2 决策变量设计：不只是无功那么简单

既然是有功-无功协同优化，决策变量就必须兼顾两个维度。我设计的决策变量向量包含四类：

这里有个工程实现上的麻烦：变压器分接头是离散量，实际运行时是整数挡位，但优化算法中直接处理整数变量会大大增加求解难度。我在第一版代码里直接把它当作连续变量处理，优化结束后再就近取整，回代潮流验证约束。实践证明，在分接头档位不多的情况下，这种“连续求解加就近取整”的方法误差不大，工程上是可接受的。

2.3 约束体系：电气互联系统比纯电网复杂在哪

约束条件是电气互联系统优化模型里最需要花心思的部分。纯电网无功优化的约束无非是潮流平衡、节点电压上下限、发电机无功上下限、变压器变比范围、补偿装置容量这些。而电气互联系统在这个基础上，又增加了天然气网侧的一系列约束。

电力侧的等式约束是标准潮流方程，采用极坐标形式的牛顿-拉夫逊方法，节点有功和无功注入平衡方程分别为：

[ P_i = V_i \sum_{j=1}^{N} V_j (G_{ij} \cos\theta_{ij} + B_{ij} \sin\theta_{ij}) ]

[ Q_i = V_i \sum_{j=1}^{N} V_j (G_{ij} \sin\theta_{ij} - B_{ij} \cos\theta_{ij}) ]

这部分我直接用Matpower的runpf做潮流内核，通过修改MPC结构体里的发电机出力和变压器参数来注入优化变量，省去了自己编写潮流迭代和雅可比矩阵的麻烦。

天然气侧约束主要包括节点流量平衡方程、管道流量方程、气源出力上下限、节点气压上下限。天然气管网的稳态流量方程我用的是经典的Weymouth方程：

[ F_k = C_k \cdot sgn(\pi_i - \pi_j) \cdot \sqrt{|\pi_i^2 - \pi_j^2|} ]

其中 (\pi_i) 和 (\pi_j) 是管道两端节点的压力，(C_k) 是管道常数。这个方程的非线性很强，直接放进优化模型里会让问题迅速变成高度非凸，实际计算中很容易陷入局部最优。

电气互联系统的耦合约束是模型的核心难点。燃气轮机的燃料消耗来自天然气网，排放出力的同时消耗天然气；电转气设备则相反，消耗电力生产天然气，两套网络的潮流方向相互关联。我在模型里把电转气设备和燃气轮机作为两个双向耦合装置，用一组线性耦合约束来表述它们对电力网和天然气网的影响，用系数矩阵把天然气流量和电力功率关联起来。

2.4 求解思路：为什么不直接用内点法

模型搭好之后我考虑过多个求解路径。商业求解器像Gurobi处理混合整数线性规划很强，但这个模型本质上是一个混合整数非线性规划（MINLP），直接丢给求解器很难收敛。内点法在纯无功优化里表现不错，但加上天然气网络的非线性方程之后，雅可比矩阵的结构变得极其复杂，实现难度上升了一个量级。

最后我采用了粒子群算法（PSO）作为外层优化器，配合内层Matpower潮流计算进行约束校验和罚函数处理。选择这条路的主要原因有三点：第一，PSO不要求目标函数和约束函数可导，对Weymouth方程这类非线性约束比较宽容；第二，实现代码量可控，核心循环加上适应度评估不到两百行；第三，对多峰非凸问题的全局搜索能力比基于梯度的算法强。

当然，粒子群也有它自己的问题，早熟收敛、参数敏感、计算量大，这些坑我在第四节会细说。这里先给出我最终采用的求解框架，整体流程就是“外层粒子群生成决策变量 → 内层潮流计算评估安全约束 → 罚函数处理越限 → 迭代更新粒子位置”。

3. Matlab代码实现：从数学模型到可运行程序

3.1 数据准备与系统建模

我用了修改后的IEEE 30节点系统作为电力侧测试算例。之所以选这个系统，一是它规模适中，既有输电网典型特征又不至于计算量过大；二是Matpower自带case30数据文件，读取修改非常方便，适合做算法验证。

天然气网数据需要自己构建。我搭了一个简化版六节点天然气系统，包含两个气源节点、四个负荷节点和五条天然气管道，其中一个燃气轮机节点和电网节点17关联，一个电转气设备和电网节点21关联。把天然气网数据写成结构体数组，字段包括节点编号、节点类型（气源/负荷）、气压上下限、管道编号、管道长度和直径、管道常数等。

Matlab里数据初始化的代码组织方式是这样的：

% 电力系统数据 mpc = loadcase('case30.m'); % 读入IEEE 30节点数据 % 天然气系统参数 ngNode = struct('id', 1:6, ... 'type', {'slack','load','load','load','load','load'}, ... 'pressureMin', [4.2, 0.8, 0.8, 0.8, 0.8, 0.8], ... 'pressureMax', [5.8, 4.5, 4.5, 4.5, 4.5, 4.5]);

这里有个小经验，天然气压力我用了标幺值，基准压力取6.0MPa，实际运算中所有压力都不必带上单位数量级，计算方便很多。

3.2 粒子群主程序和目标函数封装

粒子群的主程序逻辑很清晰，核心就是速度位置更新和适应度评估。需要注意的一个技术细节是，不同决策变量应该有独立的越界处理方式。发电机有功出力和无功出力是连续变量，直接截断在上下界之间就行；变压器变比虽然是离散的，但我在优化过程中先按连续处理，最后再就近取整；而补偿装置的无功出力需要保证在容量范围内。

我贴一下主循环的核心代码片段：

for iter = 1:maxIter for i = 1:nPop % 速度更新 v(i,:) = w * v(i,:) + c1 * rand * (pbest(i,:) - pop(i,:)) + ... c2 * rand * (gbest - pop(i,:)); % 位置更新与边界处理 pop(i,:) = pop(i,:) + v(i,:); pop(i,:) = min(max(pop(i,:), lb), ub); % 调用适应度函数 fitnessNew = fitnessFunc(pop(i,:), mpc, ngData); if fitnessNew < fitnessPbest(i) pbest(i,:) = pop(i,:); fitnessPbest(i) = fitnessNew; if fitnessNew < fitnessGbest gbest = pop(i,:); fitnessGbest = fitnessNew; end end end w = wMax - (wMax - wMin) * iter / maxIter; % 惯性权重线性递减 end

这里 (w) 是惯性权重，初始值取0.9，结束时取0.4，这样前期注重全局搜索，后期加强局部精修。

目标函数封装是核心中的核心。我在 fitnessFunc 里面先解码决策变量，然后调用 runpf 计算潮流，再根据潮流结果计算网损、碳排放、弃电量和约束越限惩罚。

function f = fitnessFunc(x, mpc0, ngData) % 从决策向量中还原各变量 nGen = length(mpc0.gen(:,1)); Pg = x(1:nGen); Qg = x(nGen+1:2*nGen); tap = x(2*nGen+1:2*nGen+nTap); qc = x(2*nGen+nTap+1:end); % 修改MPC结构体里的发电机参数 mpc0.gen(:, 2) = Pg; % PG mpc0.gen(:, 3) = Qg; % QG % 调用潮流计算 results = runpf(mpc0); % 网损、碳排放、弃风弃光惩罚计算 P_loss = sum(results.branch(:, 14) + results.branch(:, 15)); % 线路有功损耗 C_carbon = results.gen(:, 2)' * co2Factor; C_curtail = initialWind - Pg(windIdx); % 注意只有windIdx对应风电机组 % 电压越限惩罚项，通过罚函数处理约束 voltagePunish = sum(max(0, results.bus(:,8) - vMax).^2) + ... sum(max(0, vMin - results.bus(:,8)).^2); f = w1 * (P_loss / P_loss_base) + ... w2 * (C_carbon / C_carbon_base) + ... w3 * (C_curtail / C_curtail_base) + ... lambda * voltagePunish; end

这段代码就是整个优化过程的核心。我实际跑的时候，每代粒子数和进化代数分别设置在60和300，一次完整优化大约需要二十分钟，在可接受范围内。

3.3 天然气网约束的处理方式

电气互联系统中天然气网约束不能像纯电网罚函数那样大而化之地处理，因为天然气压力若越限，不仅意味着不经济，还可能是物理上不可行的。我的处理办法是把天然气优化变量（气源出力和电转气功率）也纳入粒子群决策向量，在适应度函数中用Weymouth方程求天然气潮流，再对气压越限做惩罚。

天然气潮流求解没有像纯电网那样成熟的Matlab工具箱，我手写了一个节点压力迭代函数。思路是给定气源节点压力，估算管道流量初值，迭代修正内部节点压力直到满足所有节点流量平衡。写成Matlab函数大概几十行。

实际体会是，天然气子问题因为规模小，收敛很快，反倒是和电力网耦合的迭代过程更需要小心。燃气轮机消耗天然气量是随其有功出力变化的，电转气设备产生天然气量也和其消纳的电力相关，这两组交互变量如果不先给定初值，潮流计算容易振荡。我的经验是把燃气轮机的天然气消耗和电转气的天然气产出也用罚函数先放宽，经过几次迭代后限制逐步收紧，这样收敛性和稳定性好很多。

3.4 结果可视化与对比分析

优化结束后有两类图需要画：一是优化前后的电压剖面图，二是优化过程的收敛曲线图。电压剖面图直接对比基准情况、仅无功优化、有功无功协同优化三种工况下的各节点电压标幺值。收敛曲线则用于判断粒子群是否真正收敛，以及有没有早熟迹象。

figure; plot(1:nBus, V_base, 'k-o', 'LineWidth', 1.5); hold on; plot(1:nBus, V_qonly, 'b-s', 'LineWidth', 1.5); plot(1:nBus, V_pq, 'r-^', 'LineWidth', 1.5); xlabel('节点编号'); ylabel('电压标幺值'); legend('基准情况', '仅无功优化', '有功-无功协同优化'); grid on;

对比结果方面，我跑出来的典型数据是：仅做无功优化时网损下降约9%左右，电压合格率有所提升，但碳排放基本没有变化，甚至还有轻微上升；而有功-无功协同优化后，网损下降约15%，电压偏移量显著减小，碳排放量和弃风弃光量也都有可观下降。这说明协同优化相比传统两步法的优势确实是实打实的，不是模型里的虚假收益。

4. 实操中踩过的坑与排查技巧

4.1 潮流不收敛：第一步就要排查电网参数

这个坑几乎所有人都会踩。我第一次把优化变量注入到MPC结构体里，然后调用runpf，结果一连几次算例都报错，提示潮流不收敛，我一度以为是粒子群参数不对，后来才发现是发电机有功出力的边界设置出了问题。

Matpower的case30数据里，所有发电机的有功上下界是有特定范围的，必须保证你生成的随机粒子初始值都落在 (P_{\min} \leq P \leq P_{\max}) 之内。如果粒子初始值里有超出边界的，潮流第一次迭代就可能发散。排查方法是单独写一个小脚本，检查所有粒子的初始值是否严格满足爬坡约束和上下限约束，先确保每一组随机初始解都能计算潮流，再谈优化。

另外，变压器变比的修改也不是直接改MPC结构体里branch那一列那么简单的。Matpower的branch矩阵中变压器支路的第9列和第10列对应变比和移相角，但是变压器变比超过上下界会导致雅可比矩阵奇异，表现为潮流计算在变压器支路附近反复振荡。我的处理办法是先固定变压器变比跑一次潮流，确认系统本身是收敛的，再放开变量。如果系统本身不收敛，问题绝对不在优化算法，而在数据设定上。

4.2 粒子群早熟收敛：一个防不胜防的问题

粒子群最大的问题就是早熟收敛，尤其是这种变量维度高、非线性强的模型。我第一轮实验跑出来的优化结果，多次独立运行后最优解相差很大，说明算法已经陷入了局部最优。调参方向我给出几条经验。

惯性权重 (w) 的取值很关键。之前我用固定的0.5，发现后期粒子很难精细化搜索，总是稳定在一个次优解附近。改成线性递减从0.9降到0.4之后，优化质量明显改善，但代价是收敛速度变慢。学习因子 (c_1) 和 (c_2) 都取1.5左右是比较常规的配置，但如果发现粒子种群多样性不够，可以适当增大 (c_1) 而减小 (c_2)，让粒子更多地向自身历史最优学习而不是直接冲向全局最优。

种群规模和迭代代数的初始值不要拍脑袋。我在公用的大规模算例上做过测试，粒子数从30增加到60可以明显提高解的质量，但超过80之后基本没有提升，纯增加计算时间。迭代代数300代以内一般足够，如果300代后目标函数还在持续下降，说明初始的粒子群位置生成不够好，要考虑用Sobol序列或者拉丁超立方采样生成初始种群，而不是纯随机均匀分布。

4.3 罚函数系数怎么调：经验值分享

罚函数法是启发式算法处理约束的主流手段，但罚系数设置不当会让优化结果完全失真。我最早使用固定罚系数，数值调大了之后，粒子的搜索压力全集中在满足约束上，反而忽略了目标函数本身的优化；调小了之后又可能出现电压越限、无功越限但目标函数值很好看的“假最优”。

我给出一套比较稳妥的动态罚函数策略：迭代初期罚系数取较小值，允许粒子在可行域边界附近探索；迭代后期逐步增大，迫使最终结果严格满足约束。具体实现是在粒子群迭代循环中每代更新罚系数：

lambda = lambda0 * (iter / maxIter)^2;

初始 (\lambda_0) 我取的是1，但这个值要参考目标函数各分量的量级动态设置。有一种通用做法是先用不含约束惩罚的目标函数跑十几代，观测目标函数最大的量级，然后把罚系数设定为目标函数典型值的1.5倍左右。这么做的好处是罚项既不会被淹没，又不会喧宾夺主。

如果约束越限实在无法完全消除，建议检查一下约束本身有没有矛盾。比如某个负荷节点的电压下限设置在0.95，但该节点接入的分布式光伏容量又非常大，低负荷高光伏出力时电压自然升高到1.05以上。这种物理层面的矛盾不是罚系数能解决的，必须调整算例参数或者增加无功补偿装置容量。

4.4 电气互联耦合迭代不收敛的专项排查

电气互联系统相比纯电网多了一个天然气网的边界条件。我遇到过的典型情况是：燃气轮机在某一时刻有功出力设定很高，天然气网为了供应该气量，管道流量已经接近饱和，此时管道末端压力下降得很厉害，甚至低于气压下限。这种越限不是罚函数逐步收紧就能解决的，必须直接限制燃气轮机的出力上限。

我的做法是在粒子群初始化的边界设置阶段，先进行一次天然气网的最大供气能力分析，根据管道参数算出各燃气节点的最大进气流量，再通过热值折算成对应的燃气轮机最大有功出力，把这个值作为决策变量上界。这样粒子群在初始化时就不会产生超出天然气网物理能力的个体，后续迭代也更容易收敛。

另外，电转气和燃气轮机两个耦合设备在同一节点同时出现时要注意逻辑一致性。电转气在消耗功率的同时产出天然气，如果它产出的天然气又被同一节点的燃气轮机消耗掉，就会形成闭环补气，目标函数很可能出现一个“伪最优”，看似网损降低了很多，实际是两个耦合设备的内部循环在自嗨。排查方法很简单：检查耦合节点的天然气流量净注入值是否为0，如果系统在优化后出现了明显的闭环流动，就要在约束中增加净注入量为0的条件，或者把耦合设备分配到不同节点。

5. 给新手的几条实用建议

说几个个人体会比较深的事情。

第一，先从简单模型开始，不要一上来就追求完整的气网模型。我最初接手这个方向时，恨不得把所有设备都建进模型里，结果代码写了半个月，算例跑出来的结果根本无法解释。后来退回一步，用标准的IEEE 30节点加一个六节点天然气系统，先把电力潮流和天然气潮流的接口打通，确认优化算法在简单模型上表现稳定，再逐步增加设备类型和约束条件，过程顺畅了很多。

第二，一定要相信Matpower但也要验证Matpower。我有一次优化结果里电压波形特别奇怪，局部出现了尖峰，排查了很久才发现是我修改MPC结构体时把某个节点的无功负荷字段意外改了，导致潮流计算前提已经失真。后来我养成一个习惯，每次修改MPC结构体后先做一个快照保存，再用case30原始数据做对照潮流，确认电压幅值和相角的偏差在合理范围内，再继续优化步骤。

第三，权重系数不要照搬文献。很多论文里的权重设置是基于特定算例调出来的，换一个系统或者换一组负荷曲线可能完全失效。我的做法是先用等权重跑一遍，看三个分目标的数值范围，然后根据工程重要性适当调整权重，每调一次跑一遍并整理对比表。虽然多一点工作量，但最后的结论可靠得多。

第四，保存每一代种群的中间结果。我吃过一次大亏，一次跑了十几个小时的数据因为忘了保存中间种群，最后电脑断电全部丢失，只能从头再来。现在我在粒子群迭代里每隔50代就把种群和最优解保存一次，即使中途崩溃也能恢复，成本很低但收益极大。

第五，审慎对待“降碳”效果的量化和展示。模型中通过排放因子把碳排放转成目标函数的一部分，这本质上是工程层面的简化手段，实际系统的碳排放核算远比这个复杂。在撰写报告时，我会明确区分“模型目标里的碳排放代价”和“真实系统的碳排放计量”，这个边界如果模糊了，评审专家一眼就能看出来。

电气互联系统的有功-无功协同优化这个方向，难点不在某一个环节，而在于把电力网、天然气网和优化算法很好地咬合在一起。这篇写出来的模型和代码还是相对基础的版本，后面还可以在这个框架上扩展多时段滚动优化、考虑风光出力的不确定性、加入储能系统参与调节等。如果后面有时间，我会继续把这些扩展实现在这个框架上，到时候再和大家分享更深入的实战经验。