基于Yalmip和二阶锥松弛的电气热综合能源日前优化调度实战

去年帮一个团队做园区级综合能源日前优化调度，程序框架就是标题里那套组合：Matlab 写主程序，Yalmip 建模，Cplex/Gurobi 求解，模型层面最关键是二阶锥（SOCP）松弛。整套程序跑通不算难，难的是把电气热三条网络的约束揉进同一个优化问题里，还能让求解器稳定给出可信解。这篇文章把我做这块的经验完整写出来，从为什么用二阶锥、三网模型怎么拆、Yalmip代码怎么组织，到调试时最容易踩的数值坑，全部过一遍。适合正在写综合能源优化程序、打算用 Yalmip 复现论文但被非凸约束卡住的人参考。

1. 电气热综合能源程序到底在优化什么：从调度对象到数学问题

1.1 三条网络各自的“主角”和变量

所谓电气热综合能源优化，本质上是在一个时间尺度（通常是日前 24 小时）下，同时决策电网、天然气网、热网里的运行变量，目标是让系统总运行成本最低，同时满足三条网络的物理约束和所有设备的出力范围。

电网侧的变量主要是各节点注入有功/无功功率、支路潮流、节点电压幅值。气网侧的变量是气源出流、节点气压、管道流量。热网侧则是热源出力、节点供回水温度、管道热水流量。这三者的时间耦合来自储能设备和设备的爬坡约束，空间耦合则来自燃气轮机、P2G（电转天然气）、热电联产、电锅炉这类跨网设备。

很多刚接触这块的人会犯同一个错误：把三条网当成三个独立优化问题分开算，再人工拼起来。这样做在强耦合场景下根本得不到可行解，因为电网的负荷可能是气网通过燃气轮机供给的，热网的出力又反过来约束气网的消耗量。必须在一个模型里联立求解。

1.2 耦合设备才是综合能源的题眼

综合能源系统区别于单一电网优化的核心，是那些跨网设备。我常用的一组典型耦合关系如下表所示，这些设备的模型大多是线性或分段线性的，但它们的变量会同时出现在多个网络的节点平衡方程里，这是建模时最需要注意的索引关系。

每次写耦合约束前，我会先在纸上画一张“变量-节点-网络”对照表，用 Excel 列清楚每个设备输入变量属于哪条网络、从哪个节点取能、输出注入哪条网络、到哪个节点。否则写代码时特别容易把索引串位，程序报错排半天才发现是燃气轮机的电气节点编号写错了。

1.3 一个典型日前调度的目标函数长什么样

目标函数一般取系统总运行成本最小化：

• 向上级电网购电的成本（分时电价，按节点和时段累计）
• 从气网购气的成本（按气源计价）
• 设备启停成本（对燃气轮机、热电联产这类可启停机组）
• 可再生能源弃用惩罚（光伏/风电优先消纳的软约束）

我做的程序里还加了一个很小的电压/气压越限惩罚项，目的不是真的鼓励越限，而是给求解器一个“软台阶”，避免因为硬约束边界导致病态收敛。这种做法在工程上很常见，但要注意惩罚系数不能取得过大，否则会改变最优解的性质。

2. 为什么一定要用二阶锥：非凸约束的“坑”和松弛的逻辑

2.1 非线性潮流方程的麻烦

如果是纯直流电网模型，问题是一个线性规划，好解得很。但电气热综合能源系统里，电网需要体现电压和无功，天然气管道流量和气压是非线性关系，这些都让模型变成非凸。

非凸问题带来两个直接麻烦：第一，商用求解器不能保证找到全局最优，内点法很容易陷在局部最优解里；第二，求解时间和初值选择高度相关，换一组负荷数据可能就不收敛。对科研复现来说，这种“时好时坏”的程序是最要命的。我不止一次见过同学的程序跑出一组离谱的机组出力，原因就是非线性求解器从一个差的初始点掉进了局部陷阱。

二阶锥松弛的思路是：把非凸可行域“撑大”成一个凸的锥体，代价是放宽了某些等式约束为不等式，但换来整个问题变成凸优化。凸优化的好处是，Cplex/Gurobi 这类求解器可以又快又稳地拿到全局最优解，而且可以从数学上保证最优性。

2.2 DistFlow的凸松弛推导

辐射状配电网潮流最常用的模型是 Branch Flow Model，也就是 DistFlow。定义每条支路 i-j 的潮流为 P_ij、Q_ij，电流幅值平方为 l_ij，节点电压幅值平方为 v_i，那么支路方程是：

v_j = v_i - 2(r_ij P_ij + x_ij Q_ij) + (r_ij^2 + x_ij^2) l_ij

l_ij = (P_ij^2 + Q_ij^2) / v_i

第二个等式是非凸的，因为 P_ij^2 + Q_ij^2 除以 v_i 不是凸函数。二阶锥松弛把它不等式化：

l_ij >= (P_ij^2 + Q_ij^2) / v_i

这个不等式等价于标准二阶锥约束：

code复制|| [2*P_ij; 2*Q_ij; v_i - l_ij] ||_2 <= v_i + l_ij

这个形式在 Yalmip 里可以直接用 cone() 写，非常方便。松弛的意义在于：原来可行域是一条曲线上的点，现在变成整个锥体内部的点。只要最优解恰好落在锥的边界上，松弛前后的解就是完全一致的，这个问题下面专门讲。

2.3 天然气Weymouth方程怎么转锥

天然气管道流量和气压的关系常用 Weymouth 稳态方程描述：

f_ij^2 = C_ij^2 (Π_i - Π_j)

其中 Π 是气压平方，C 是管道常数。这个式子右边的括号代表两端气压平方差，决定了气流方向。它同样是二次非凸等式。

做松弛时把它放宽为不等式：

f_ij^2 <= C_ij^2 (Π_i - Π_j)

这样形成一个旋转二阶锥。实际写代码时为了用标准二阶锥 cone() 形式，可以用一个常见的过参数化技巧：令 w = C_ij^2 ΔΠ，把约束写成以下等价形式：

code复制|| [2*f_ij; w - 1] ||_2 <= w + 1

因为展开后是 4 f_ij^2 + (w-1)^2 <= (w+1)^2，化简后正好得到 f_ij^2 <= w，也就是 f_ij^2 <= C_ij^2 ΔΠ。这个“加减 1 的过参数化”是SOCP建模里的常规操作，自己推导一遍就明白，别直接抄网上的代码，否则换单位后很容易出错。

2.4 松弛的“紧性”意味着什么

二阶锥松弛解出来以后，必须验证松弛是不是“紧”的。所谓紧，就是求出的解让原来那个不等式约束取到了等号。比如电网松弛后的 l_ij 是否真的等于 (P_ij^2+Q_ij^2)/v_i，气网松弛后的 f_ij^2 是否真的等于 C_ij^2 ΔΠ。

大多数论文里的辐射状电网算例，DistFlow 的 SOCP 松弛在很多工况下都是紧的。气网的 Weymouth 松弛紧性则没有电网那么稳定，尤其在管道流量很小、气压差很小的时候，解可能直接落在锥的内部。遇到这种情况，就需要用惩罚项或迭代法把解“拉回”边界，我在第 5 节会展开讲具体做法。

3. 三网联立的建模拆解与耦合约束

3.1 电网侧：DistFlow与节点电压约束

电网侧我用的是辐射状配电网，DistFlow 模型加节点注入平衡约束，节点功率平衡方程写成：

P_i = P_ij - sum(P_jk) - r_ij l_ij

Q_i = Q_ij - sum(Q_jk) - x_ij l_ij

这里 P_i、Q_i 是该节点所有注入源的净注入，包括上级电网购电、燃气轮机发电、光伏出力，扣除本地负荷。节点电压幅值平方有上下限，通常取 0.95^2 到 1.05^2。支路有功/无功潮流也要设上限，防止线路过载。

写代码时有个小技巧：把线路参数 r、x、容量上限都整理成矩阵，用分支索引直接切片，比一条条写约束省事得多，而且不容易漏。

3.2 气网侧：节点平衡、管道方程与压缩机

气网节点平衡方程的核心是：

气源注入 + 流入流量 - 流出流量 + P2G产气 = 气负荷 + 燃气轮机耗气 + 燃气锅炉耗气

每个气节点的气压平方 Π 有上下界。管道流量方向在程序中需要处理：如果管网拓扑允许双向输气，最好引入两个非负变量 f_ij^+ 和 f_ij^-，让净流量 f_ij = f_ij^+ - f_ij^-，再用大 M 约束保证至少一个方向为 0。

压缩机站是另一个建模点。稳态模型中常常用压缩比约束：

Π_j <= CR_max * Π_i

或者用更细的分段线性模型。如果只是学模型和验证算法，先用固定压缩比或干脆忽略压缩机，把精力放在 SOCP 松弛本身，等基本框架稳定后再加压缩机也不迟。

3.3 热网侧：节点法模型与流量问题

热网比电网和气网多了一个时间常数和温度变量。常用的节点法模型把热网分成供水管网和回水管网两层，每个热负荷节点的热功率方程为：

Φ_i = c_p * m_i * (T_s,i - T_r,i)

问题是，流量 m_i 和温度 T_s、T_r 乘积是双线性项，不能直接放进 SOCP。这里有两条路：

第一条，采用质调节的简化假设，即各节点流量在优化周期内固定不变，那么流量是已知参数，热功率和温差就是线性关系。这个模型对大多数规划研究够用，也是很多论文的实际做法。

第二条，如果一定要优化流量，就需要把流量的变化离散化，用混合整数线性化处理双线性项，或者做迭代：先给一组流量初值，解 SOCP，然后把得到的温度代回去修正流量，重复迭代直到收敛。第二种方法工程上更实用，但要注意迭代过程可能不收敛，需要设置最大迭代次数和松弛阈值。

3.4 耦合设备模型：燃气轮机、P2G与热泵

耦合设备建模的原则很简单：把设备的输入变量写到对应网络的节点平衡方程里，把输出变量写到另一条网络的节点平衡方程里，中间用效率或可行域约束关联。例如燃气轮机：

F_gas,gt = P_gt / η_gt

这个等式把气网消耗和电网出力关联起来。热电联产机组如果考虑电热可行域，可以用顶点表格的凸组合形式变成一组线性约束，Yalmip 里直接写成矩阵不等式即可。

P2G 和电锅炉更简单，前者是电转气：F_gas,p2g = η_p2g * P_p2g，后者是电转热：Φ_eb = η_eb * P_eb。热泵加一个 COP 系数就行。这些约束全部线性，不会破坏凸性。

4. Yalmip代码骨架：从变量声明到二阶锥约束的落地写法

4.1 变量体系与标幺值约定

代码组织的第一件事是标幺值。我不建议在 Yalmip 里直接用真实单位混着写，尤其是气压高达几十万帕、而功率只有几兆瓦的时候，数值矩阵的条件数会很难看，求解器收敛会很吃力。通常的做法是电功率全部用 MW，电压幅值用标幺值，气压用 MPa 量级，热量用 MWth，流量用 kg/s。这样各个变量的数量级都落在 0.01 到 100 之间。

变量定义示例：

matlab复制nb = 33;            % 电网节点数
nl = 32;            % 电网支路数
ng = 6;             % 气网节点数
nh = 6;             % 热网节点数
T  = 24;            % 调度时段

Pg   = sdpvar(nb, T);      % 各时段购电有功
Qg   = sdpvar(nb, T);      % 各时段购电无功
V2   = sdpvar(nb, T);      % 电压幅值平方
I2   = sdpvar(nl, T);      % 支路电流幅值平方
Pij  = sdpvar(nl, T);      % 支路有功
Qij  = sdpvar(nl, T);      % 支路无功
Pi2  = sdpvar(ng, T);      % 气网节点气压平方
Fg   = sdpvar(nl_gas, T);  % 气网管道流量

每行变量都加上注释，告诉未来的自己这个变量属于哪条网络、哪个节点、什么单位。这句话听着像老生常谈，但程序规模一大，没有注释的变量表会让人崩溃。

4.2 二阶锥约束在Yalmip里的几种写法

DistFlow 约束可以这样写：

matlab复制for k = 1:nl
    i = branch(k, 1);
    j = branch(k, 2);
    r = branch(k, 3);
    x = branch(k, 4);
    for t = 1:T
        % 电压降方程
        Constraints = [Constraints,
            V2(j,t) == V2(i,t) - 2*(r*Pij(k,t) + x*Qij(k,t)) + (r^2+x^2)*I2(k,t)];

        % 二阶锥松弛
        Constraints = [Constraints,
            cone([2*Pij(k,t); 2*Qij(k,t); V2(i,t) - I2(k,t)], V2(i,t) + I2(k,t))];
    end
end

气网 Weymouth 约束转标准锥：

matlab复制for k = 1:nl_gas
    i = gas_branch(k, 1);
    j = gas_branch(k, 2);
    C = gas_branch(k, 3);
    for t = 1:T
        dp = C * (Pi2(i,t) - Pi2(j,t));
        % 过参数化二阶锥
        Constraints = [Constraints,
            cone([2*Fg(k,t); dp - 1], dp + 1)];
    end
end

如果管道方向固定，可以不加双向变量，直接按拓扑方向写。如果方向不固定，千万记得引入方向变量，并在锥约束外面用大 M 约束保证 f_ij^+ 和 f_ij^- 不同时为正。

4.3 求解器调用与典型运行规模

目标函数和约束组装好后，调用求解器的写法非常短：

matlab复制ops = sdpsettings('verbose', 2, 'solver', 'cplex', 'showprogress', 1);
ops.cplex.mip.tolerances.mipgap = 1e-4;
ops.cplex.sim.tolerances.optimality = 1e-6;

result = optimize(Constraints, Objective, ops);

我实际测试的典型规模是 33 节点电网 + 6 节点气网 + 6 节点热网 + 24 时段，变量数大概在八千到一万五千之间，Cplex 求解 SOCP 一般在几十秒到几分钟内可以结束。如果规模更大，建议先用 4 小时或者 8 小时时段测试模型正确性，再跑全天。

求解器选择上，我比较过几种，简单结论如下表：

如果模型里带有整数变量（机组启停、分段线性化），实际上是一个混合整数二阶锥规划（MISOCP），Cplex 和 Gurobi 都能处理，只是求解时间会明显增加。我的经验是尽可能把整数变量数量控制在几百个以内，不然调试会非常磨人。

5. 求解器调试实录：数值病态、M值与卡死解的真实战场

5.1 数值尺度是最容易被低估的坑

我开始做这个程序的时候，气网压力单位直接用了 Pa，结果 Cplex 求解时一直报数值警告，解的精度也一塌糊涂，有时候明明是可行的负荷数据，却提示无解。后来把所有压力统一改成 MPa，把功率统一改成 MW，问题立刻好转。

数值尺度不只是单位问题，还包括约束系数的大小对比。比如一个约束里电压平方是 1 左右，而管道常数 C 可能是几千，这种系数差距会让求解器内部的预处理很辛苦。我的做法是：每个网络内部先标幺化，再通过耦合设备的效率公式统一能量单位，最后在目标函数里用真实价格系数还原成本量纲。

5.2 大M处理双向管道时的手感

双向气体管道建模要用大 M 法：

f_ij = f_ij^+ - f_ij^-
0 <= f_ij^+ <= M * b_ij
0 <= f_ij^- <= M * (1 - b_ij)

其中 b_ij 是二进制变量。M 取多少是个经验问题。取太小可能把实际可行的流量上限截断；取太大会让线性松弛变得很差，求解变慢甚至数值不稳定。我的经验是 M 取管道物理最大流量的 1.2 倍，而不是随便输入一个大数。

如果你发现同样的模型在 Cplex 里时快时慢，很大概率是 M 值取大了。这时候把 M 缩小一点，往往立竿见影。

5.3 松弛结果不紧怎么补救

如果验证发现 SOCP 松弛不紧，常规做法有几种。

第一，加惩罚项：在目标函数里加上对“松弛量”的惩罚。对电网，惩罚 (l_ij - (P_ij^2+Q_ij^2)/v_i)；对气网，惩罚 (C_ij^2 ΔΠ - f_ij^2)。惩罚系数先给一个小值，解完检查，如果不紧再逐步增大。这个方法的局限是惩罚太大会把目标函数带偏。

第二，迭代固定法：先求解松弛问题，把部分状态变量固定，比如固定管道流向，再作为新的约束加入模型重新求解。这个过程相当于把大问题拆成两步，实际操作很多情况下能拿到紧解。

第三，换用精确非线性求解器做后验验证。我通常用 Ipopt 把 SOCP 解作为初值代入非线性模型，如果误差很小，就说明解是可信的。这一步对写论文很重要，因为审稿人一定会盯着松弛紧性问。

5.4 遇到无可行解的排查链路

程序报 infeasible 的时候，不要慌，按顺序排查：

1. 先去掉所有二阶锥松弛，只求解等式约束，看是不是模型本身矛盾。这一步跑不通就说明基础变量或拓扑有误。

2. 再检查节点功率平衡方程的符号。我犯过最多的是购电符号方向反了，导致电从负荷流向母线，模型当然无解。

3. 然后查耦合设备的索引。燃气轮机消耗的气网节点编号和发电的电网节点编号是不是对应上了；P2G 是否把产气写到了正确的气网节点。

4. 最后用 Yalmip 的 check 命令查看每条约束的残差，找最大误差的约束，基本就是问题所在。

这个链路我自己用了很多次，每次都能在十分钟内定位问题。强烈建议你在代码里给每条约束加注释和分组，这样 check 输出的时候能直接看到是哪一组约束出问题。

6. 从跑通到可信：程序验证方法与算例完整性

6.1 小算例验证：IEEE 33节点加三网耦合是起步配置

我建议起步算例不要大，我给团队的配置是 IEEE 33 节点配电网 + 比利时天然气 6 节点网 + 一个 6 节点热水网，24 个时段，三处耦合点：一个燃气轮机接在电网某节点和气网某节点之间，一个电锅炉接在电网节点供热负荷，一台热电联产供一部分热负荷。这个规模既能暴露模型问题，又不会因为算不动而拖慢调试。

小算例还有一个好处：可以手工估算几组结果。比如把气价设得极高，程序应该自动把所有燃气轮机出力降到最低，改用购电和电锅炉；如果程序没有出现这个趋势，说明某个耦合约束写错了，目标函数没有真正指导调度决策。

6.2 回代检验与非线性模型对比

SOCP 求解完之后，我习惯做三个验证：

1. 把最优解代回原非凸等式，计算最大相对误差，看是否在 1% 以内。

2. 把同一算例交给非线性求解器重解一次，对比各时段总成本和主要设备出力曲线，偏差不应超过几个百分点。

3. 改变负荷水平，观察最优解是否平滑变化。如果负荷只加了一点点，设备出力却大幅跳变，很可能模型里有坏的数值解或松弛不紧。

三个验证都通过，我才敢说这个程序的结果是可信的。做研究写论文时，这部分验证过程本身也是很有价值的内容，审稿人往往非常看重。

6.3 可以继续扩展的方向

框架稳定之后，往这几个方向扩展都比较顺手：加入储能设备（电池、储热罐、储气库），引入可再生能源出力的不确定性做鲁棒优化或随机优化，把单一成本目标扩展为碳排放和成本的多目标，或者在热网侧加入动态管道温度模型。这些扩展都不需要推翻二阶锥框架，而是在现有约束基础上继续加层。

我个人的体会是，二阶锥松弛在电气热综合能源里的角色很像一个稳定的底座：它让原本难解的非凸模型变成一个工程上可用的凸优化问题，又比纯线性模型多保留了电压、气压、温度的网络物理特征。很多研究结论，包括最优化调度策略、多能互补效益、设备容量配置，都是建立在这个底座之上的。如果底座本身数值不稳、松弛不紧，后面所有结论都不可靠。所以别急着堆复杂约束，先把 SOCP 这一层做扎实。