做电力系统研究的朋友对“无功优化”这四个字应该都不陌生。以前我们做无功优化思路很清晰给定有功调度结果再去调整无功补偿设备、变压器分接头让电压合格、线损最小。这套思路在传统电网里跑了很多年已经很成熟。但近几年有一个趋势越来越明显——随着“碳中和”目标落地风电、光伏、燃气轮机大规模并网电网和天然气网开始深度耦合有功调度和无功电压控制已经没法再当成两个独立的优化问题来解了。换句话说在一个电气互联系统里你在某台燃气轮机上多发了有功立刻会引起无功支撑能力的变化进而影响电压、线路损耗甚至碳排放指标。本文围绕的这种“牵一发动全身”的耦合关系构建了电气互联系统下的有功-无功协同优化模型并用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里跑不通、调不稳一切都是空谈。把非线性项做合理的线性化把约束按模块组织把每个变量的物理意义都搞清楚大部分实现层面的问题都可以迎刃而解。如果你也正在做类似的课题希望这篇文章能帮你少走几步弯路。