上个月帮一位师弟调算例他手里有现成的有功经济调度模型加了无功优化之后解算时间从不到1秒涨到了七八分钟而且电压曲线算出来明显不对。我再一翻他的约束发电机无功上限给的0.3 p.u.变压器分接头根本没建进去天然气网和电网的耦合变量只在目标函数里象征性乘了个系数。这哪是无功优化是给有功调度挠了个痒痒。说实话类似的情况在“碳中和”背景下的综合能源系统研究里特别常见。“电气互联系统有功-无功协同优化”这个题目论文里人人都写代码能跑通的却不多。核心问题不是大家不会写目标函数而是对“有功-无功为什么要协同”“电气互联系统的耦合界面在哪里”“非凸约束怎么在Matlab里落地”这几件事缺一条完整的逻辑链。这篇文章我就把自己跑通的一个模型从头拆一遍问题背景、数学建模、二阶锥凸化、Yalmip代码实现、松弛紧度检验和调参避坑全给你捋清楚。适合刚接触综合能源系统优化的研究生也适合已经跑通纯电网潮流优化、想往电-气耦合方向扩展的朋友参考。1. 碳中和目标下无功优化为什么不再是“附属品”1.1 高比例新能源带来的电压支撑缺口传统电力系统里同步发电机是绝对的电压支撑主力。励磁系统调无功本质上是就地平衡无功功率、维持节点电压。碳中和目标之下风光机组大规模接入情况开始变得微妙。光伏逆变器和双馈风机虽然理论上能发无功但实际运行中它们的无功能力高度依赖有功出力。以光伏为例中午大发有功的时候逆变器容量接近饱和无功裕度反而最小傍晚出力降下来了无功裕度回来了但系统又往往处于负荷爬坡阶段最需要无功支撑的恰恰是另一个节点。风机也是这样很多双馈机型在低电压穿越期间会主动吸收无功导致故障后电压恢复更加困难。更麻烦的是随着部分在运煤电机组逐步进入灵活性改造或退役序列系统里的同步转动惯量和无功备用都在下降。可以这样理解无功功率的“银行”里储户少了可提款的人没少。于是出现一个现象——有功调度计划做得再漂亮电压曲线照样越限线损也压不下来。1.2 从“有功经济调度”到“有功-无功协同优化”的转变逻辑传统经济调度只回答一个问题在满足负荷和网络约束的前提下各机组有功出力怎么分配最便宜。它把无功当成“事后校验”先算有功再校核电压不行就切负荷或者调分接头。但在新能源高渗透系统里这个先后顺序行不通了。原因在于有功和无功在潮流方程里是强耦合的。比如一条长距离输电线路输送的有功越大线路两端无功损耗越多末端电压跌落越厉害。此时要让系统继续安全运行可能得让近端的燃气机组多发一点有功、同时把无功出力控制在更宽的范围或者调整变压器分接头改变无功分布。这些决策互相牵连分开优化必然得到次优解。碳中和又给这个优化加了一条新约束碳排放。燃气机组的碳排放强度比煤电低一半还多但要让它顶上来就得考虑气网的供气能力。于是电网侧的“有功-无功协同”进一步扩展为“电网-气网协同”。这就是电气互联系统有功-无功协同优化的完整动机既要有功经济、电压合格又要系统整体碳排放可控还要气网供气可行。2. 电气互联系统到底在互联什么2.1 先澄清一个容易误解的名词“电气互联系统”里的“电气”指的是电力系统和天然气系统不是“电气工程”那个电气。更严谨的写法是“电-气互联系统”即 Integrated Electricity and Gas System。英文里用 and 连接中文口头简称成“电气”一不留神就歧义了。和审稿人或者导师交流时第一次出现最好写全称避免被误读。这两个网络的耦合靠的是两类关键设备燃气轮机电网侧看是发电机输出有功和无功气网侧看是天然气负荷消耗燃料。燃气机组的运行方式直接决定两个网格之间的功率流和气流。电转气Power-to-GasP2G电网侧看是负荷消耗电能气网侧看是气源注入天然气。新能源大发时段富余的电通过电解水制氢、甲烷化合成天然气回注气网。工业和建筑里的另一类耦合设备是热电联产和燃气锅炉但那是电-热互联的范畴电气互联优化里往往当边界条件处理不放进决策模型。2.2 为什么耦合进来的是天然气而不是别的很多人会问要支撑新能源用储能不是更直接吗储能确实行但它和天然气网的物理特性不同。天然气网的独特之处在于管道本身具有储气能力也就是“线储”。一条长距离输气管道首端进气量和末端出气量可以短期不一致靠管存缓冲。这个特性让气网在应对燃气机组快速调节时有一定柔性同时也让它变成一个有记忆性的网络——上一时刻管道里的存气会影响下一时刻的气压和流量约束。所以电气互联系统的优化不能只看静态断面。虽然很多学术文章为了可解性把模型简化成单时段我接下来给的代码骨架也是单时段的但心里要清楚严格意义上这是一个多时段耦合问题。如果要做日内滚动优化就得在模型里加入管道管存动态方程和储气库状态方程。从运行成本角度看气源购气成本与煤电煤耗成本在目标函数里是同等级的决策量。碳中和场景下煤电出力受碳排放约束挤压燃气机组出力占比上升气网侧的约束会越来越紧。如果不把气网建进去燃气机组的“天然气供应”就是无穷大变量系统等于少了一大块真实约束算出来的最优解大概率不可行。3. 协同优化模型的数学骨架3.1 目标函数把碳成本“晒”进运行成本里这个模型的目标函数我建议写成四部分相加煤电煤耗成本、气源购气成本、碳排放惩罚成本以及可选的弃新能源惩罚成本。数学形式如下[ \min ; \sum_{i \in \Omega_G} \left( a_i P_{G,i}^2 b_i P_{G,i} c_i \right) \sum_{s \in \Omega_S} c_{gas,s} F_{S,s} \lambda_{CO_2} \left( \sum_{i \in \Omega_G} e_i P_{G,i} \right) ]其中 ( P_{G,i} ) 是机组有功出力( F_{S,s} ) 是气源产气量( e_i ) 是机组碳排放强度tCO₂/MWh( \lambda_{CO_2} ) 是碳价元/tCO₂。这个形式的妙处在于碳成本不是作为约束强行卡某个排放上限而是折算成经济信号进入目标函数。机组调度会在“煤电便宜但碳排高”和“气电贵但碳排低”之间自动权衡。实际调参时我会对 ( \lambda_{CO_2} ) 做 20 到 200 元/t 的敏感性扫描看煤电出力占比和总排放的变化拐点。这个拐点就是系统的“经济碳价”比单纯用约束限额更有工程参考意义。3.2 电网侧约束DistFlow 模型与无功调节手段电网部分我推荐用支路潮流模型也就是常说的 DistFlow/Branch Flow 模型。它比传统极坐标潮流方程更适合做凸优化因为它用支路有功 ( P_{ij} )、支路无功 ( Q_{ij} )、支路电流平方 ( l_{ij} )、节点电压平方 ( v_i ) 作为变量结构规整。支路潮流方程的核心有三组支路功率平衡首端功率等于末端功率加支路损耗。支路压降方程( v_i - v_j 2(r_{ij} P_{ij} x_{ij} Q_{ij}) - (r_{ij}^2 x_{ij}^2) l_{ij} )。电流-功率-电压关系( P_{ij}^2 Q_{ij}^2 v_i l_{ij} )。第三组是二次等式也是模型非凸性的来源之一。后面凸化处理就是针对它。节点功率平衡方程把支路潮流和机组注入、负荷连接起来[ \sum_{j \in \delta(i)} P_{ij} P_{G,i} - P_{D,i}, \quad \sum_{j \in \delta(i)} Q_{ij} Q_{G,i} - Q_{D,i} ]无功调节手段按响应速度和控制方式分几类在模型里处理方式完全不同手段变量类型响应速度模型处理发电机无功出力连续变量毫秒级上下限约束变压器分接头离散变量分钟级整数或二进制变量并联电容器/电抗离散变量开关操作二进制/整数变量SVC/SVG连续变量毫秒级连续变量容量约束我的建议是第一版模型先只做发电机无功上下限和节点电压约束把基座跑通。确认无问题后再加变压器分接头和电容器投切。一次性把所有离散变量都扔进去解算时间会爆炸而且遇到不可行问题时根本没法定位是哪个约束惹的祸。3.3 气网侧约束Weymouth 方程与节点气压天然气网的建模比电网“粗糙”一些核心变量是节点气压 ( \pi_i ) 和管道流量 ( f_{ij} )。稳态气网的流量-压差关系由 Weymouth 方程描述[ f_{ij}^2 C_{ij}^2 \left( \pi_i^2 - \pi_j^2 \right) ]流量方向由气压差决定气压高的节点向气压低的节点流动。节点气流平衡方程为[ \sum_{j \in \delta(i)} f_{ij} F_{S,i} - F_{D,i} - F_{GT,i} ]其中 ( F_{D,i} ) 是气网负荷( F_{GT,i} ) 是燃气机组耗气量。气源有上下限节点气压也有上下限。长距离输气管道还要考虑压缩机站压缩机升压需要消耗一部分天然气或电力。模型第一版可以忽略压缩机或者把它简化成固定升压比等主问题稳定了再精细化。Weymouth 方程是大麻烦它包含气压平方项和流量平方项是非凸二次等式。后面会单独讲怎么处理。3.4 气电耦合约束燃气轮机和 P2G 的数学表达燃气机组的耗气量近似为有功出力的线性函数[ F_{GT,i} \alpha_i P_{GT,i} \beta_i ]注意单位换算( P_{GT,i} ) 单位是 MW( F_{GT,i} ) 通常用 m³/h 或 MJ/s必须通过天然气热值统一换算。常见热值取 35.6 MJ/m³ 上下具体按你用的数据手册来。P2G 设备的表达式反过来[ F_{P2G,i} \eta_{P2G} \cdot P_{P2G,i} / H_{NG} ]( \eta_{P2G} ) 是电转气综合效率当前工程水平在 50% 到 65%别用太乐观的数值。P2G 的投运给电网侧增加了一个可调负荷给气网侧增加了一个可调气源这是电网和气网之间唯一的双向耦合通道非常关键。4. 求解策略把非凸约束“驯服”成可解形式4.1 二阶锥松弛DistFlow 的经典凸化上一节提到支路潮流模型里的 ( P_{ij}^2 Q_{ij}^2 v_i l_{ij} ) 是非凸等式。解决办法是把它松弛成不等式[ P_{ij}^2 Q_{ij}^2 \leq v_i l_{ij} ]这个不等式可以等价写成标准二阶锥形式[ \left| \begin{bmatrix} 2P_{ij} \ 2Q_{ij} \ l_{ij} - v_i \end{bmatrix} \right|2 \leq l{ij} v_i ]在 Yalmip 里这一条写起来非常干净Constraints [Constraints, cone([2*Pij(k); 2*Qij(k); l(k) - V(i)], l(k) V(i))];松弛之后模型从非凸二次约束优化变成二阶锥规划SOCP如果还有变压器分接头、电容器投切这类整数变量整体就是混合整数二阶锥规划MISOCP。Gurobi、Mosek、CPLEX 都能直接解。但“松弛”意味着原问题可行域被放大了。如果最优解恰好落在被放大的那部分区域求出来的解对应的是虚假的物理状态需要专门检验松弛紧度。这个我在第 6 节详细讲。4.2 天然气网 Weymouth 方程的近似处理Weymouth 方程的凸化比 DistFlow 棘手。文献里做法五花八门我实际项目里用过两种说下适用场景。第一种是固定流向分段线性化。先根据历史运行方式或一次直流潮流估计固定每条管道的流向然后把“流量平方-气压平方差”的关系曲线分段线性化。这个方法实现简单、求解快缺点是流向固定后失去了最优性适合管道数量少、流向明确的场景。第二种是双向流量二进制变量建模。每条管道引入方向二进制变量配合 Big-M 约束处理双向流动。模型更精确但整数变量数量翻倍求解难度上升。对于 30 条以上管道的算例我一般先试第一种如果发现最优解里存在方向反转的管道再回去升级成第二种。具体做增量线性化时Yalmip 里可以用 SOS2 变量但 Gurobi 对 SOS2 的求解效率一般。我更推荐自己写增量线性化incremental linearization用普通二进制变量做分段选择。代码稍长但可控性强。4.3 求解器选型与参数设置模型是 MISOCP 的话Matlab 自带的intlinprog帮不上忙必须外接商用求解器。我这些年用得最多的是 Gurobi其次是 Mosek。两者对二阶锥约束的支持都很成熟许可证在高校里也容易申请。调用方式很简单Yalmip 负责建模求解器负责解ops sdpsettings(verbose, 2, solver, gurobi); ops.gurobi.mipgap 1e-4; % MIP 间隙阈值 ops.gurobi.timelimit 3600; % 单次求解时间上限 ops.gurobi.numericsfocus 1; % 数值稳定性优先对于几十个节点的电气互联系统模型规模不算大重点是整数变量数量。如果 MIPGap 设得太严格比如 1e-6大算例可能几个小时出不来。我的经验是 1e-4 足够工程使用学术论文里可以设到 1e-5再严意义不大。5. Matlab 代码实现从数据表到求解出图5.1 环境与数据准备代码实现的环境是 Matlab Yalmip 求解器。Yalmip 装好之后在命令行跑yalmiptest确认所有求解器都被正确识别。新版 Matlab2022 到 2026 那几版对 Yalmip 的支持没有本质差异只要注意下载对应版本的 Yalmip 压缩包即可。电网数据我建议直接用 Matpower 格式。loadcase(case39)能一次性拿到节点表、发电机表、支路表而且支路表里电阻、电抗、对地电纳都是标幺值省去大量手工录入工作。气网数据没有现成标准格式一般自己建结构体字段包括节点编号、管道首末端、管道直径、长度、压缩因子、气源上限、气负荷。5.2 电网建模核心代码下面这段是核心建模骨架我尽量把关键约束完整给出可以直接抄进你的主程序。数据来自 Matpower 的 mpc 结构体。mpc loadcase(case39); nb size(mpc.bus, 1); nl size(mpc.branch, 1); ng size(mpc.gen, 1); fb mpc.branch(:, 1); % 支路首端节点 tb mpc.branch(:, 2); % 支路末端节点 r mpc.branch(:, 3); x mpc.branch(:, 4); % 负荷 Pd mpc.bus(:, 3); % 有功负荷 Qd mpc.bus(:, 4); % 无功负荷 % 发电机节点映射 gen_bus mpc.gen(:, 1); gen_Pmax mpc.gen(:, 9); gen_Qmax mpc.gen(:, 5); gen_Qmin mpc.gen(:, 6); % 定义变量 V sdpvar(nb, 1); % 节点电压平方 l sdpvar(nl, 1); % 支路电流平方 Pij sdpvar(nl, 1); % 支路有功 Qij sdpvar(nl, 1); % 支路无功 Pg sdpvar(ng, 1); % 机组有功 Qg sdpvar(ng, 1); % 机组无功 Constraints []; % 支路潮流约束 for k 1:nl i fb(k); j tb(k); % 二阶锥松弛 Constraints [Constraints, ... cone([2*Pij(k); 2*Qij(k); l(k) - V(i)], l(k) V(i))]; % 支路压降方程 Constraints [Constraints, ... V(i) - V(j) 2*(r(k)*Pij(k) x(k)*Qij(k)) - (r(k)^2 x(k)^2)*l(k)]; end % 节点功率平衡 for i 1:nb inj_idx find(fb i); out_idx find(tb i); g_idx find(gen_bus i); if isempty(g_idx) Constraints [Constraints, ... sum(Pij(inj_idx)) - sum(Pij(out_idx)) -Pd(i)]; Constraints [Constraints, ... sum(Qij(inj_idx)) - sum(Qij(out_idx)) -Qd(i)]; else Constraints [Constraints, ... sum(Pij(inj_idx)) - sum(Pij(out_idx)) Pg(g_idx) - Pd(i)]; Constraints [Constraints, ... sum(Qij(inj_idx)) - sum(Qij(out_idx)) Qg(g_idx) - Qd(i)]; end end % 机组约束 Constraints [Constraints, 0 Pg gen_Pmax]; Constraints [Constraints, gen_Qmin Qg gen_Qmax]; % 电压约束 Constraints [Constraints, 0.95^2 V 1.05^2]; % 平衡节点电压锚定 ref find(mpc.bus(:, 2) 3); if ~isempty(ref) Constraints [Constraints, V(ref(1)) 1.0^2]; end注意电压约束里的平方我定义 ( V ) 是电压幅值平方所以上下限也要平方。很多人第一次写这个模型直接写0.95 V 1.05电压一下飞到 1.05 的平方根附近算出来结果怎么都对不上。5.3 气网建模与耦合约束气网建模比电网稍“脏”一些因为 Weymouth 方程的处理方式直接影响模型类型。下面给出固定流向分段线性化版本的骨架。% 气网参数结构体 % gas.bus: [编号 基准气压 气压下限 气压上限 气负荷] % gas.pipe: [首端节点 末端节点 管道系数C 流量上限] nbg size(gas.bus, 1); npl size(gas.pipe, 1); ngs size(gas.source, 1); Pi sdpvar(nbg, 1); % 节点气压平方 fs sdpvar(npl, 1); % 管道流量 Gs sdpvar(ngs, 1); % 气源产气量 % 管道约束 for k 1:npl i gas.pipe(k, 1); j gas.pipe(k, 2); % 固定流向首端气压高于末端 Constraints [Constraints, Pi(i) Pi(j)]; % 流量平方采用分段线性化这里简化为直接约束 Constraints [Constraints, fs(k)^2 gas.C(k)^2 * (Pi(i) - Pi(j))]; % 示意实际需分段 end % 节点气流平衡 for i 1:nbg in_pipe find(gas.pipe(:, 2) i); out_pipe find(gas.pipe(:, 1) i); src_idx find(gas.source(:, 1) i); Constraints [Constraints, ... sum(fs(in_pipe)) - sum(fs(out_pipe)) Gs(src_idx) - gas.bus(i, 5)]; end % 气压和气源上下限 Constraints [Constraints, gas.bus(:, 3).^2 Pi gas.bus(:, 4).^2]; Constraints [Constraints, gas.source_lb Gs gas.source_ub];我必须强调上面代码里fs(k)^2 gas.C(k)^2 * (Pi(i) - Pi(j))这一行是示意不是最终可用的公式。真正的分段线性化要用 SOS2 或增量二进制变量把流量平方在可行区间内近似成折线。你要是图省事直接写这个等式丢给求解器Gurobi 会把它当非凸二次约束要么拒绝求解要么直接弹错。耦合约束反而最简单% 燃气机组耗气量 Constraints [Constraints, F_gt alpha .* Pg_gt beta]; % P2G注入气网 Constraints [Constraints, F_p2g eta .* Pp2g / HNG];这两个线性约束把电网变量和气网变量串联起来模型闭环。5.4 调算子与结果输出模型组装完毕调用optimizeops sdpsettings(verbose, 2, solver, gurobi); ops.gurobi.mipgap 1e-4; ops.gurobi.timelimit 1800; diagnostics optimize(Constraints, Objective, ops); if diagnostics.problem 0 value(V) value(Pg) value(Pi) else disp(求解失败检查约束和求解器日志); checkset(Constraints) % Yalmip 自带约束检查 endcheckset是 Yalmip 最实用的调试工具之一它会把每条约束的残差打印出来快速定位是哪条约束过不了。6. 模型调参、松弛紧度检验与避坑地图6.1 锥松弛到底紧不紧三步检验法整个模型能不能用关键就看二阶锥松弛紧不紧。所谓“紧”就是松弛后的不等式解出来刚好取等号。如果不取等号说明模型放大了可行域算出来的“最优解”在原始问题上根本不存在。我的检验方法是三步走。第一步对每条支路计算松弛间隙gap abs(V(i).*l(k) - (Pij(k).^2 Qij(k).^2)) ... ./ max(1e-8, Pij(k).^2 Qij(k).^2);第二步统计最大间隙和平均间隙。工程上最大间隙小于 1e-3 基本可以接受小于 1e-4 就相当健康了。第三步如果间隙太大在目标函数里给每条支路的松弛量加一个很小的惩罚项比如Objective Objective 1e-4 * sum(V(fb).*l - (Pij.^2 Qij.^2));罚系数从 1e-4 开始试探逐步加大直到松弛间隙达标。注意系数不能太大否则它会反过来主导目标函数革命又变成帮凶。6.2 碳价、罚系数和目标函数权重怎么调碳价的设置直接影响模型行为我的标准操作是先跑三组对照第一组碳价设 0得到纯经济调度基准解。第二组碳价设 50 元/t看机组组合出力和碳排放总量变化。第三组碳价设 150 元/t看气源产气量、燃气机组耗气量是否触碰气网约束上限。这三组对比能快速定位气网的“瓶脖子”。如果第二组里气网约束已经饱满了说明气网参数比电网先一步限制了降碳空间需要回头检查是否该扩展管道路径或增加 P2G 容量。罚系数方面弃新能源惩罚要设得比碳价更有分量。否则优化结果可能是“为了减碳宁愿少用风光”这就不合理了。可再生能源的边际成本接近零弃掉它既费钱又费碳罚系数至少要比碳价高一个量级。6.3 高频踩坑与排查指引我把这两年调试电气互联模型遇到的典型问题整理成一张表照着症状查原因比较快。现象可能原因处理办法求解器返回 Infeasible气网气压上下限过严或燃气机组耗气量超过了气源上限放宽气压下限或把燃气机组最大有功调低再试求解器报 NaN气压变量初始值为 0Weymouth 方程里除以气压平方给气压设置一个宽松的下界比如 0.1 倍的基准气压电压曲线边缘化普遍贴着上限电压目标函数缺失只有成本项在目标函数里加电压偏移惩罚或用软约束解算时间骤增几个钟头跑不完变压器分接头、电容器整数变量太多先用连续变量松弛找初值再固定一部分离散变量松弛间隙不达标反复振荡罚系数过小锥约束取等号没有激励逐步加大松弛惩罚系数每调一次观察间隙变化气网管道流量和气压方向不一致固定流向与实际最优流向相反放开流动方向限制加方向二进制变量这些坑我基本全踩过其中气压初始值为 0 那个最隐蔽。气网节点气压如果从 0 开始Weymouth 方程里 ( \pi_i^2 ) 一阶导数为 0求解器在数值上非常容易出问题。解决办法就是给气压变量一个务实的下界别从零开始探索。6.4 整体调下来的几句实话这套模型我前后调了两周最大的体会是“别想着一步到位”。我自己的路径是先只跑电网 DistFlow 无功优化确认锥松弛紧度达标再把天然气网接进来但燃气机组只当固定气负荷最后才把 P2G 和双向耦合加进去。每一步都有可靠的中间结果做锚点出了问题知道该怪谁。第二个体会是数据单位比数学公式更容易坑人。功率用 MW 还是 p.u.气压用 kPa 还是 MPa天然气流量用 m³/h 还是 MJ/s稍不留神就差了三个数量级。我现在的习惯是在数据录入阶段就统一标幺化基准功率取 100 MVA基准气压取管道设计压力所有参数全部转成 p.u. 后再进模型。这样做的好处不仅是数值稳定后续调参时量级也直觉。最后一个想分享的细节是MISOCP 模型的初值质量对求解速度影响极大。先用连续松弛模型解一遍把解出来的整数变量值赋给原模型作为初始解可以让整数搜索的路径大幅缩短。Yalmip 里可以用assign手动赋初值配合sdpsettings(gurobi.branchdir, 1)这类参数控制分支方向大算例能省下一半以上的求解时间。