电气互联系统有功-无功协同优化：建模、二阶锥松弛与Yalmip实现

上个月帮一位师弟调算例，他手里有现成的有功经济调度模型，加了无功优化之后，解算时间从不到1秒涨到了七八分钟，而且电压曲线算出来明显不对。我再一翻他的约束：发电机无功上限给的0.3 p.u.，变压器分接头根本没建进去，天然气网和电网的耦合变量只在目标函数里象征性乘了个系数。这哪是无功优化，是给有功调度挠了个痒痒。

说实话，类似的情况在“碳中和”背景下的综合能源系统研究里特别常见。“电气互联系统有功-无功协同优化”这个题目，论文里人人都写，代码能跑通的却不多。核心问题不是大家不会写目标函数，而是对“有功-无功为什么要协同”“电气互联系统的耦合界面在哪里”“非凸约束怎么在Matlab里落地”这几件事缺一条完整的逻辑链。这篇文章我就把自己跑通的一个模型从头拆一遍：问题背景、数学建模、二阶锥凸化、Yalmip代码实现、松弛紧度检验和调参避坑，全给你捋清楚。适合刚接触综合能源系统优化的研究生，也适合已经跑通纯电网潮流优化、想往电-气耦合方向扩展的朋友参考。

1. 碳中和目标下，无功优化为什么不再是“附属品”

1.1 高比例新能源带来的电压支撑缺口

传统电力系统里，同步发电机是绝对的电压支撑主力。励磁系统调无功，本质上是就地平衡无功功率、维持节点电压。碳中和目标之下，风光机组大规模接入，情况开始变得微妙。

光伏逆变器和双馈风机虽然理论上能发无功，但实际运行中它们的无功能力高度依赖有功出力。以光伏为例，中午大发有功的时候，逆变器容量接近饱和，无功裕度反而最小；傍晚出力降下来了，无功裕度回来了，但系统又往往处于负荷爬坡阶段，最需要无功支撑的恰恰是另一个节点。风机也是这样，很多双馈机型在低电压穿越期间会主动吸收无功，导致故障后电压恢复更加困难。

更麻烦的是，随着部分在运煤电机组逐步进入灵活性改造或退役序列，系统里的同步转动惯量和无功备用都在下降。可以这样理解：无功功率的“银行”里，储户少了，可提款的人没少。于是出现一个现象——有功调度计划做得再漂亮，电压曲线照样越限，线损也压不下来。

1.2 从“有功经济调度”到“有功-无功协同优化”的转变逻辑

传统经济调度只回答一个问题：在满足负荷和网络约束的前提下，各机组有功出力怎么分配最便宜。它把无功当成“事后校验”，先算有功，再校核电压，不行就切负荷或者调分接头。

但在新能源高渗透系统里，这个先后顺序行不通了。原因在于有功和无功在潮流方程里是强耦合的。比如，一条长距离输电线路输送的有功越大，线路两端无功损耗越多，末端电压跌落越厉害。此时要让系统继续安全运行，可能得让近端的燃气机组多发一点有功、同时把无功出力控制在更宽的范围，或者调整变压器分接头改变无功分布。这些决策互相牵连，分开优化必然得到次优解。

碳中和又给这个优化加了一条新约束：碳排放。燃气机组的碳排放强度比煤电低一半还多，但要让它顶上来，就得考虑气网的供气能力。于是，电网侧的“有功-无功协同”进一步扩展为“电网-气网协同”。这就是电气互联系统有功-无功协同优化的完整动机：既要有功经济、电压合格，又要系统整体碳排放可控，还要气网供气可行。

2. 电气互联系统到底在互联什么

2.1 先澄清一个容易误解的名词

“电气互联系统”里的“电气”，指的是电力系统和天然气系统，不是“电气工程”那个电气。更严谨的写法是“电-气互联系统”，即 Integrated Electricity and Gas System。英文里用 and 连接，中文口头简称成“电气”，一不留神就歧义了。和审稿人或者导师交流时，第一次出现最好写全称，避免被误读。

这两个网络的耦合，靠的是两类关键设备：

• 燃气轮机：电网侧看是发电机（输出有功和无功），气网侧看是天然气负荷（消耗燃料）。燃气机组的运行方式直接决定两个网格之间的功率流和气流。
• 电转气（Power-to-Gas，P2G）：电网侧看是负荷（消耗电能），气网侧看是气源（注入天然气）。新能源大发时段富余的电，通过电解水制氢、甲烷化合成天然气，回注气网。

工业和建筑里的另一类耦合设备是热电联产和燃气锅炉，但那是电-热互联的范畴，电气互联优化里往往当边界条件处理，不放进决策模型。

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} ]

无功调节手段按响应速度和控制方式分几类，在模型里处理方式完全不同：

我的建议是：第一版模型先只做发电机无功上下限和节点电压约束，把基座跑通。确认无问题后，再加变压器分接头和电容器投切。一次性把所有离散变量都扔进去，解算时间会爆炸，而且遇到不可行问题时根本没法定位是哪个约束惹的祸。

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确认所有求解器都被正确识别。新版 Matlab（2022 到 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 调算子与结果输出

模型组装完毕，调用optimize：

ops = 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 自带约束检查 end

checkset是 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 高频踩坑与排查指引

我把这两年调试电气互联模型遇到的典型问题整理成一张表，照着症状查原因比较快。

这些坑我基本全踩过，其中气压初始值为 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)这类参数控制分支方向，大算例能省下一半以上的求解时间。