做电力系统研究的朋友，对“无功优化”这四个字应该都不陌生。以前我们做无功优化，思路很清晰：给定有功调度结果，再去调整无功补偿设备、变压器分接头，让电压合格、线损最小。这套思路在传统电网里跑了很多年，已经很成熟。但近几年有一个趋势越来越明显——随着“碳中和”目标落地，风电、光伏、燃气轮机大规模并网，电网和天然气网开始深度耦合，有功调度和无功电压控制已经没法再当成两个独立的优化问题来解了。换句话说，在一个电气互联系统里，你在某台燃气轮机上多发了有功，立刻会引起无功支撑能力的变化，进而影响电压、线路损耗，甚至碳排放指标。本文围绕的这种“牵一发动全身”的耦合关系，构建了电气互联系统下的有功-无功协同优化模型，并用Matlab代码把它完整落地。读完这篇文章，你能搞清楚模型怎么搭、目标函数怎么设、约束怎么建模、Matlab代码该怎么组织，以及我在实际调参和跑算例时踩过的坑。适合电网调度、分布式能源研究方向的工程师，也适合电力系统方向的研究生作为课题参考。
1. 项目概述：为什么有功和无功非要放在一起优化
1.1 传统分离思路为什么不够用
先说说我们以前的常规做法。传统电力系统中，有功调度主要解决“谁发电、发多少”的问题，核心目标是经济性，约束是系统频率和潮流不过载。而无功优化则是另一套逻辑，主要盯住电压，通过调整无功补偿容量、变压器变化和发电机端电压，让各节点电压落在允许范围之内，顺便把有功网损降下来。
这两件事在过去之所以能分开做，是因为系统规模相对可控，发电厂出力基本稳定，负荷变化有规律，电压问题通常可以通过就地无功补偿处理。但现在的场景已经变了。分布式光伏和风电接入后，有功出力随机波动；燃气轮机电转气等新型设备，让电网和天然气网之间不再是“一买一卖”的关系，而是双向耦合。你在天然气侧做一次气源调度调整，很可能会改变燃气轮机发电出力，而这个出力变化又会立刻传导到无功电压层面。
举个最直观的例子。一个分布式光伏大发的中午，系统电压往往偏高，这时候需要无功补偿设备吸无功；到了傍晚光伏出力骤降，电压又偏低，需要无功补偿设备发无功。如果优化模型还只考虑有功调度，不考虑这些电压变化对系统运行成本和碳排放的影响，最后给出的调度方案在工程上根本执行不下去，或者执行了也会付出额外的调节代价。
1.2 碳中和目标给优化模型带来的新变量
碳中和目标对优化问题的影响，不只是“在目标函数里加一个碳排放项”这么简单。它实际上改变了系统运行的底层逻辑：碳成本一旦计入决策目标，燃气轮机和燃煤机组的相对经济性就会改变，天然气的使用量、电网的购电策略都要重新洗牌。
更关键的是，碳约束与无功电压问题是会相互传导的。比如，为了降低碳排放，系统会倾向于让高排放的火电机组少发有功，让燃气轮机或者P2G设备多出力。但火电机组少发有功往往意味着它的无功出力极限也随之变化，原有无功支撑点被削弱，电网电压就可能越限。这时候你必须把无功优化纳入全局决策，才可能在满足碳目标的同时守住电压安全。
再比如P2G（电转气）设备，它本质上是把富裕的电能转化成天然气，在碳视角下是个“负排放”的好东西，但它本身是一个负荷，不仅要消耗有功，还需要一定的无功支撑。你在系统里增加一个P2G，相当于给电网加了个不小的用电大户，潮流分布、电压分布都会跟着变。所以，在碳中和目标下做电气互联系统优化，有功-无功协同不是可选项，而是必选项。
2. 模型构建：目标函数、约束条件和碳机制
2.1 目标函数：发电成本、气源成本、碳成本怎么叠
我把这个优化问题建模成一个追求总运行成本最小的问题。目标函数由三块构成：电网侧的发电成本、天然气网侧的气源成本，以及碳排放成本。如果系统里还有从外部电网购电的通道，还得把购电成本也加进去。
电网侧发电成本，我习惯用二次函数拟合。每台发电机的有功出力 (P_{Gi})，成本函数写成：
[ C_{G,i}=a_i P_{Gi}^2 + b_i P_{Gi} + c_i ]
其中 (a_i) 一般是个很小的正数，用来描述机组的煤耗曲线凸性；(b_i) 决定边际成本的主体；(c_i) 是空载成本。在Matlab里可以直接用向量方式写，不必写成循环。
天然气网侧的气源成本就简单一些，一般用线性函数：
[ C_{S,j}=c_{gas,j} \cdot F_{S,j} ]
其中 (F_{S,j}) 是第 (j) 个气源的注入流量，(c_{gas,j}) 是单位气量价格。
碳排放成本这块，我没有简单用一个固定碳价系数乘总排放量，而是采用了阶梯碳交易机制。这样更贴近目前实际碳市场的运作方式：系统会先获得一个免费排放配额 (E_{quota})，如果实际排放 (E_{total}) 低于配额，则不需要支付碳成本；一旦超出配额，超出部分按照阶梯价格收费，超额越多，单位碳价越高。这样建模的好处在于，优化器会自动在“多花钱降碳”和“省成本超排”之间做权衡，也更符合实际政策导向。
2.2 有功-无功耦合的数学表达
模型里最核心的是电气互联系统的耦合约束。电网侧，我采用的交流潮流约束，节点 (i) 的有功和无功注入平衡如下：
[ P_{Gi} + P_{GT,i} - P_{P2G,i} - P_{Li} = V_i \sum_{j} V_j (G_{ij}\cos\theta_{ij} + B_{ij}\sin\theta_{ij}) ]
[ Q_{Gi} + Q_{C,i} - Q_{Li} = V_i \sum_{j} V_j (G_{ij}\sin\theta_{ij} - B_{ij}\cos\theta_{ij}) ]
这两个式子看起来复杂，但它的物理含义很直白。左边是注入节点的净功率，右边是节点通过线路送出的功率。它天然就是有功和无功耦合的数学表达：同一个节点电压幅值 (V_i) 和相角 (\theta_i)，同时决定了有功和无功的分布。
天然气网侧的约束我用的是稳态管道气流模型。对于一条连接节点 (u) 和节点 (d) 的天然气管道，其流量 (f_k) 与两端压力 (p_u, p_d) 的关系为：
[ f_k = C_k \sqrt{p_u^2 - p_d^2} ]
这个方程是非线性的，处理起来比较麻烦，后面我会讲怎么在Matlab里做线性化。除了管道方程，天然气管网还有气源注入上下限、节点压力上下限、负荷节点供气量必须满足等约束。
电气互联的耦合点主要是燃气轮机和P2G设备。燃气轮机消耗天然气发电，天然气流量 (F_{GT}) 和发出的有功功率 (P_{GT}) 之间满足能量转换关系：
[ F_{GT} = \frac{P_{GT}}{\eta_{GT} \cdot LHV} ]
其中 (\eta_{GT}) 是发电效率，(LHV) 是天然气低位热值。P2G设备则反过来，消耗电能产生天然气：
[ F_{P2G} = \frac{\eta_{P2G} \cdot P_{P2G}}{LHV} ]
值得强调的是，P2G在消耗有功 (P_{P2G}) 的同时，也需要一定的无功支撑，可以把它看作一个功率因数较低的负荷。这个细节如果漏掉，电压约束在P2G节点附近就很容易越限。
2.3 碳交易机制的建模细节
碳排放量的计算也要分源。发电机组的排放量，直接用有功出力和排放因子相乘；燃气轮机的排放则和它消耗的天然气流量相关，本质上是燃料燃烧产生的 (CO_2)。总的碳排放量可以写成：
[ E_{total} = \sum_{i} e_{g,i} P_{Gi} + \sum_{GT} e_{GT} F_{GT} ]
阶梯碳价的建模，我用的是分段线性函数。假设免费配额是 (E_0)，超额部分第一阶梯价格是 (p_1)，超过一定阈值后进入第二阶梯价格 (p_2)。碳成本函数写成：
[ C_{CO2} = p_1 \cdot \min(E_{extra}, E_{threshold}) + p_2 \cdot \max(0, E_{extra} - E_{threshold}) ]
其中 (E_{extra}=\max(0, E_{total}-E_0))。
这样建模的一个额外好处是，在Matlab里我们可以把它写成一个可以微分的分段函数，用sdpvar或者数值优化工具箱都能比较方便地处理，不需要额外引入整数变量去表示阶梯切换的逻辑。
3. Matlab实现方案：从数学到代码
3.1 求解路径怎么选
拿到这个模型之后，第一个要决策的问题是求解器怎么选。我自己试过几条路，说下实际感受。
如果你的模型把天然气管道方程做了线性化，整体问题就变成了一个混合整数线性规划（MILP）或者二次规划（QP）问题。这种情况下优先推荐用YALMIP作为建模层，背后接Gurobi或者CPLEX。YALMIP的语法非常友好，能把复杂约束直接写成类似于数学表达式的形式。学术用户一般都能拿到Gurobi的免费license，求解几百个变量的问题就是秒级甚至亚秒级的事。
如果非要保留原始交流潮流的非线性，那可以选择用Matlab自带的fmincon配合优化工具箱。但这里的坑比较多，因为交流潮流方程是非凸的，fmincon时常会卡在局部最优解上，初始值稍微给偏一点，结果就完全不同。我个人的经验是，至少先用直流潮流或者线性化潮流算出一个可行解，再把它当初始值扔给fmincon，这样成功率会高很多。
还有一种路线是群体智能算法，粒子群或者灰狼优化。这类算法不用求导，对非凸问题的适应性看起来很好，但它们的问题也很明显：约束处理很麻烦，等式约束尤其是潮流平衡方程往往只能通过罚函数来近似，罚系数调不好，收敛结果就一塌糊涂。如果系统规模大一点，动辄好几百个决策变量，群体智能算法的计算时间会膨胀到无法接受。我的建议是，能用凸优化就千万不要偷懒直接上启发式算法，省下来的时间都是自己的。
3.2 代码结构怎么组织
我在写这个项目的代码时，没有把所有内容堆到一个大脚本里，而是分模块组织。一个典型的目录结构大概是这样的：
IES_OPF/ ├── data/ │ ├── bus_data.m % 电网节点数据 │ ├── branch_data.m % 电网支路数据 │ ├── gas_node_data.m % 气网节点数据 │ └── gas_pipe_data.m % 气网管道数据 ├── model/ │ ├── build_opf.m % 构建电网潮流约束 │ ├── build_gasflow.m % 构建气网约束 │ ├── build_coupling.m % 构建燃气轮机与P2G耦合约束 │ └── build_objective.m % 构建目标函数 ├── solve/ │ └── run_main.m % 主程序入口 └── utils/ ├── linearize_pipe.m % 管道方程线性化 └── postprocess.m % 结果后处理与绘图主程序入口的思路很清晰：先加载数据，然后调用各模块构建约束，调用求解器求解，最后做后处理。
在YALMIP里，决策变量这样定义：
% 电网决策变量 P_G = sdpvar(nGen, 1); % 发电机有功出力 Q_G = sdpvar(nGen, 1); % 发电机无功出力 V = sdpvar(nBus, 1); % 节点电压幅值 Theta = sdpvar(nBus, 1); % 节点电压相角 % 气网决策变量 F_S = sdpvar(nGasSource, 1); % 气源注入量 F_GT = sdpvar(nGT, 1); % 燃气轮机耗气量 F_P2G = sdpvar(nP2G, 1); % P2G产气量 % 碳排放相关变量 Delta_E = sdpvar(1, 1); % 超配额排放量目标函数可以这样写：
% 发电成本 C_G = sum(a_coeff .* P_G.^2 + b_coeff .* P_G + c_coeff); % 气源成本 C_S = sum(c_gas .* F_S); % 购电成本（如果有外部电网） C_Buy = 0.12 * P_Buy; % 碳成本 E_total = sum(e_coeff .* P_G) + sum(e_GT .* F_GT); E_extra = max(0, E_total - E_quota); C_CO2 = p1 * min(E_threshold, E_extra) + p2 * max(0, E_extra - E_threshold); % 总目标 Objective = C_G + C_S + C_Buy + C_CO2;YALMIP一个很大的优势就在于sdpvar表达式可以直接支持平方、分段函数这类操作，它会自动把问题交给底层求解器处理。
3.3 几个关键技术点：线性化、变量归一化、初值
在实际代码实现中，有几个技术点值得单独说。
第一个是天然气管道方程的线性化。Weymouth方程里有一个根号项，直接丢给求解器很麻烦。我的做法是，先把 (p_u^2) 和 (p_d^2) 分别作为新的变量，然后对根号部分做分段线性逼近。分段线性化的思路本质上是用若干段直线方程去逼近曲线，在Matlab里可以手动设置分段断点，生成对应的线性约束。分段数取10段左右，精度已经足够工程使用，模型规模也不会膨胀得太厉害。
第二个是变量归一化。Matlab的求解器虽然能处理不同量纲的变量，但我建议把电压、相角、功率、气体流量统一标幺化。电压基值取12.66kV，功率基值取100MW或者1MW，气体流量基值取对应管道设计流量。这样做的好处是数值尺度一致，求解器的数值稳定性会好很多，各路约束也不会因为量纲差异而出现病态条件数。
第三个是初始值。如果求解器是fmincon这类需要初始值的优化器，我一般先用一个不考虑无功和碳排放的简化模型跑一遍，得到一组近似的可行解，再把这个解作为初始值代入完整模型。这样能极大避免“初始值偏离可行域太远导致收敛失败”的问题。
4. 完整实操过程：测试系统与结果分析
4.1 测试系统搭建
验证这个模型，我用了一个常见的电气互联测试系统：电网侧采用IEEE 33节点配电网，天然气网侧采用一个7节点供气系统，两者通过2台燃气轮机和1套P2G设备耦合。这个组合在电-气联合优化研究的论文中很常见，规模适中，既能体现耦合特性，又不会让Matlab计算时间太长。
电网侧的节点负荷，有有功和无功两部分数据。我把系统总负荷设置在3.7MW左右，无功负荷1.9MVar左右，峰谷时段各做了一组数据。燃气轮机的参数：单台最大有功出力0.5MW，无功出力上下限±0.12MVar，发电效率42%，单位出力的碳排放因子约0.248 tCO2/MWh。P2G设备额定功率0.3MW，转换效率60%，无功消耗按照功率因数0.85折算。
天然气网侧的节点参数，包括气源价格、气源流量上下限、节点压力上下限。气源价格取了0.32元/立方米，供气压力基准按0.8MPa设计。
4.2 代码逐步调试与运行
主程序运行的核心流程可以拆成下面这几步。
第一步是基础数据准备。把电网的母线参数和支路参数、气网的节点参数和管道参数、燃气轮机和P2G的耦合参数，全部写成Matlab数据文件。
第二步是搭建约束。我按“电网约束模块→气网约束模块→耦合约束模块→碳约束模块”的顺序依次构建。这样做的好处是排查问题方便：哪个模块报错就单独验证哪个模块，不会把所有错误混在一起。
第三步是求解。在YALMIP中，我执行：
ops = sdpsettings('solver', 'gurobi', 'verbose', 2); optimize(Constraints, Objective, ops);如果求解成功，YALMIP会返回solprog >= 0的状态。接下来我调用value(P_G)、value(Q_G)、value(F_S)等函数，把优化结果取出来。
第四步是后处理。我会绘制节点电压分布图、发电机有功无功出力图、气网流量分布图，还会把协同优化和传统分离优化的结果放在一起对比。
4.3 结果怎么分析
跑完算例后，我最看重四个指标：总运行成本、网损、电压质量、碳排放量。
在我设置的基准场景里，协同优化模型相比“先单独优化有功、再单独优化无功”的两步法，总运行成本降低了大约4%到6%，主要来自网损的下降和气源调度的优化；碳排放量下降了约12%，因为在碳价信号驱动下，优化器主动减少了高排放机组出力，让燃气轮机和P2G承担了更多调节任务。电压质量方面，协同优化的节点电压最低值比两步法提升了约0.02p.u.，各节点电压普遍更接近基准值。
这里有一点值得说：协同优化的优势并不仅仅体现在数值上，更体现在方案的可执行性上。两步法给出的结果往往在数学上可行，但在工程上需要二次调整——电压越限了，再回头重新调整补偿设备；补偿设备调完，有功的经济调度又偏离了最优。协同优化一次性把这些问题都处理掉了，这才是它真正的价值。
5. 常见问题与排查技巧实录
5.1 排查问题速查表
代码写得多了，问题就集中在几个常见点上。我整理成一个表格，对照着排查效率很高。
| 现象 | 可能原因 | 解决方案 |
|---|---|---|
| 求解器报无可行解 | 约束条件前后矛盾，或边界设置过紧 | 先用松弛版模型验证，逐步收紧约束定位矛盾来源 |
| P2G节点电压越限 | 忽略了P2G的无功消耗 | 在P2G节点的无功平衡方程中加入Q_P2G项 |
| 天然气管道流量计算出错 | Weymouth方程线性化精度不足 | 增加分段数，或对低压段单独加密断点 |
| 碳排放量异常偏低 | 排放因子单位错误，或者漏算了燃气轮机碳排放 | 统一tCO2/MWh和m³/h的单位换算 |
| fmincon收敛到坏点 | 初始值给得不好，目标函数非凸 | 先用简化模型求解，冷启动换热启动 |
| 求解超时 | 决策变量过多，MILP规模大 | 降低分段线性化段数，或者用MPS格式导出后让Gurobi并行求解 |
5.2 我踩过的几个坑
第一个坑是P2G无功负荷的遗漏。我第一次搭模型时，把P2G简单地当作一个纯有功负荷处理，结果优化结果里P2G节点电压一路滑到0.90p.u.以下，怎么调都不行。后来反应过来：P2G设备内部的电力电子变换器、压缩机运行都需要无功支撑。修正之后，给P2G节点增加了一个无功负荷项，电压立刻回到了合理范围。这个教训提醒我，在系统建模时，任何有功转换设备都必须同时检查它的无功边界。
第二个坑是阶梯碳价函数的可微性。最开始我直接用max函数叠加碳成本，在YALMIP里没问题，但如果哪天换到fmincon，这种不可导的表达式会经常触发计算精度警告。后来我把碳成本函数改成用分段多项式拟合，平滑过梯度断点，迭代速度一下子稳定了。
第三个坑是天然气管道方程线性化后的可行域退化。分段线性逼近的段数太少，会让原本宽阔的管道流量可行域被压缩成一条很窄的带状可行域，求解器很容易判定无解。我一开始只分了三段，问题规模确实小，但求解器一直报错。后来改成十个断点，并把断点密度集中在压力差较小的区间，情况才好转。
5.3 其他值得注意的工程细节
如果你的系统里有多个P2G或者多个燃气轮机，耦合变量的排列顺序会影响稀疏矩阵的生成效率。我建议在构建约束时按照“同一节点耦合设备相邻”的原则排序，这样约束矩阵的稀疏性会更好，求解器预处理的速度会有明显提升。
关于碳排放配额，需要特别注意配额分配的节点归属问题。碳配额如果分配给“系统整体”，那碳成本就是一个全局变量，所有机组共享一个配额池；如果配额是分配到机组的，那每台机组必须有自己独立的 (E_{quota})，约束结构会完全不同。这两种建模方式对优化结果的影响非常大，需要根据实际政策背景来确定，不能一概而论。
最后说一个实用小技巧：Matlab里的YALMIP求解完，可以用validated_options查看底层求解器实际用了哪些参数选项。调试阶段，把Gurobi的MIPGap设置为0，看它能否在可接受时间内收敛到最优解；如果解的质量波动大，再去检查是不是线性化精度或者约束冗余的问题。
我在实际做这个项目的过程中，最大的感受是：电气互联系统的协同优化，难点并不仅仅在于“如何把目标函数和约束写出来”，更在于“如何让求解器把问题解出来”。理论模型再漂亮，如果在Matlab里跑不通、调不稳，一切都是空谈。把非线性项做合理的线性化，把约束按模块组织，把每个变量的物理意义都搞清楚，大部分实现层面的问题都可以迎刃而解。如果你也正在做类似的课题，希望这篇文章能帮你少走几步弯路。
