从非凸到凸：综合能源系统二阶锥松弛建模与MISOCP求解

把一套含电、气、热三类能源的综合能源优化程序从“能跑”调到“跑得稳”，我前后折腾了小半年。最典型的教训是：同样的园区数据，第一版用非线性求解器直接算潮流方程，初值稍微给偏一点，CHP出力的结果就能差出15%；第二版把潮流方程线性化，速度快了，但最优解代入实际潮流一验算，线路过载都发现不了；第三版终于走到正路上来——用MATLAB加Yalmip建模，把电网DistFlow方程里的电压-电流二次项做二阶锥松弛，天然气Weymouth方程做凸松弛，热网在质调节假设下线性化，整套模型交给Cplex做混合整数二阶锥规划（MISOCP），不猜初值、收敛稳定、全局最优有保障。这篇文章就把第三版的路子从头到尾拆开：物理建模、凸松弛原理、Yalmip代码、求解器调参、松弛紧性校验，按你拿到一个问题后实际操作顺序来写，适合正在做园区综合能源、微电网调度和能源互联网方向课题的同学，也适合想从传统非线性求解转凸优化的工程师。

1. 综合能源优化问题的骨架：三条能源网络的耦合与调度边界

1.1 优化目标不是单一的电费，而是整个能源站的总账

很多第一次接触综合能源优化的人，会把目标函数直接写成“购电费用最小”。这其实只对了一半。含电气热的系统里，天然气的购买费用、CHP和锅炉的启停损耗、设备的运行维护成本，甚至弃风弃光惩罚，全部都要进目标函数。不然很容易出现一种奇怪结果：电费降下来了，气费涨得更多，总账反而难看。

以最典型的日前经济调度为例，目标函数长这样：

Objective = sum(price_e .* P_buy * dt) + ... % 购电费用 sum(price_g .* F_gas * dt) + ... % 购气费用 sum(C_start .* u_start) + ... % 启停成本 sum(C_om .* P_device); % 运维成本

其中dt是各时段长度，单位小时；price_e是分时电价，price_g是气价。注意这里有个隐蔽细节：如果电价单位是“元/MWh”，电功率是MW，时间长度是h，乘出来直接就是元；但如果电价表给的是“元/kWh”，很多人直接代进去，目标函数数值一下子就差了1000倍，求解器大概率直接报数值问题。

1.2 CHP、锅炉、电转气：物理设备如何把三条网“缝”起来

电气热三个网络不是各算各的，它们靠耦合设备连在一起。最常见的三样：

• CHP（热电联产）：天然气进去，电和热同时出来。电出力与热出力之间有一个可行的配比范围，不是随便组合。
• 燃气锅炉：天然气进去，只出热。结构简单，效率相对固定。
• 电锅炉/热泵：电进去，出热。在弃风弃光时段特别好用，相当于把多余的电转化成热储存或直接供热。

如果系统里还有P2G（电转气），那就是第四条耦合链路：电变成天然气，回灌到气网。做这类课题时，我把这些设备想象成“配货中转站”——电、气、热三条传送带各自有容量限制，CHP就是那个从气传送带取货、同时往电和热传送带放货的机器，但放货比例有约束。二阶锥规划要干的事，就是在所有传送带容量和转运规则都满足的前提下，找到总成本最低的配货方案。

1.3 为什么只做能量平衡不够，还要带上网络约束

这是很多简化模型最容易被诟病的地方。如果只写“电功率平衡 + 热功率平衡 + 气源平衡”，相当于默认所有线路和管道都是无限容量、没有压降、没有电压限制。这种模型求出来的解，拿到真实物理网络里很可能是不可行的：CHP在某时段满出力，局部节点电压可能越上限；热网末端用户可能因为管道压降不够而供热不足。

所以完整的综合能源优化必须把网络约束写进优化模型。电网要有电压约束、支路容量约束；气网要有节点压力约束、管道流量与压力关系；热网要有供回水温度和节点热功率约束。但这恰恰引出了一个问题：这些物理约束几乎都是非线性的，直接扔给求解器，就会回到我最开始说的初值敏感、局部最优的泥潭。这也是二阶锥规划在这类问题里成为主流选择的核心原因。

2. 为什么偏偏是二阶锥：凸松弛的数学依据与适用边界

2.1 DistFlow的非凸等式与“松弛为锥”的完整推导

配电网最常用的潮流模型是DistFlow，适用于辐射状网络。对一条支路(i, j)，用P_ij、Q_ij表示支路首端有功、无功功率，v_i是节点i电压幅值的平方，l_ij是支路电流幅值的平方，r_ij、x_ij是支路电阻、电抗，那么：

• 电压降落方程：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

问题出在最后这个等式。它不是凸的，直接放进优化模型会得到一个非凸可行域，求解时对初值极其敏感。二阶锥松弛的做法，就是把这个等式“放松”成不等式：

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

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

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

不信可以自己推一下：两边平方，展开以后就是P_ij^2 + Q_ij^2 <= v_i * l_ij，正好是上面的式子。

这里的直觉可以这样理解：原约束是一条弯曲的非凸曲线，松弛成锥以后可行域变大了。如果求出来的最优解恰好落在锥的“边界”上，也就是不等式取等号，那么这个解和原非凸问题完全一致，我们说松弛是“紧”的。文献里的精确性定理告诉我们，对辐射状配电网、目标函数合理、没有特殊逆向潮流场景的绝大多数情况，SOCP松弛都是紧的。工程项目的实践也基本符合这个结论。

2.2 天然气Weymouth方程和热网约束的凸化处理

天然气管道稳态流量和节点压力之间，最常用的是Weymouth方程：

f_pipe^2 = C^2 * (alpha_i - alpha_j)

其中alpha_i是节点i压力的平方，C是管道常数。这是一个非凸二次等式。工程做法同样是松弛成不等式：

f_pipe^2 <= C^2 * (alpha_i - alpha_j)

再借助一个漂亮的变换写成标准二阶锥：

d = alpha(i) - alpha(j); Constraints = [Constraints, d >= 0]; Constraints = [Constraints, cone([2*f_pipe(k); C2 - d], C2 + d)];

C2是C^2。这个写法的好处是，无论Yalmip还是Cplex/Gurobi，都把它当成真正的二阶锥约束来处理，而不是回调非线性求解器。需要说明的是，气网流向如果可能反转，上述写法就不够用，一般要引入二进制变量刻画方向，模型复杂度会明显上升。实际项目里通常可以根据气源和负荷位置预判流向，或者只在少数管道上保留方向变量。

热网是三者里最麻烦的。热功率本身是流量和温差的乘积：H = c_p * m * (T_s - T_r)。如果同时把流量m、供水温度T_s、回水温度T_r都当变量，这就是双线性项，直接突破SOCP框架。工程上常用两条路绕开：一是质调节，把流量固定在设计工况值，只优化温度；二是量调节，温度固定，只优化流量。两种做法都能让热网约束保持在线性体系内。如果确实需要两者同时可变，我建议用“固定流量解一次SOCP，再把解出来的流量代回去更新温降系数，重新求解”的外层迭代，通常两三次就能收敛。这一点后面建模部分会细说。

2.3 SOCP与LP、NLP方案的取舍

不同建模方式没有绝对优劣，只有适合和不适合。我自己的使用感受可以整理成下面这个表：

所以对含电气热的综合能源优化，我的默认选择是：连续变量走SOCP，需要机组启停就加0-1变量变成MISOCP。只有在SOCP实在没法表达某些强非线性设备特性时，才考虑对那一小块做分段线性化或外层迭代。

3. 电-气-热联合建模：从物理方程到代数约束的完整推演

3.1 电网侧建模与标幺化

电网建模第一步是统一单位。我习惯把功率统一成MW，电压统一成标幺值pu。电压幅值的平方v一般限制在0.95到1.05的平方之间，也就是0.9025到1.1025，这个范围比直接用电压百分数要直观得多。

支路数据和DistFlow约束的代码骨架是这样的：

nb = 5; nl = 4; branch = [1 2; 2 3; 2 4; 4 5]; % from, to r = [0.01; 0.02; 0.015; 0.01]; x = [0.03; 0.05; 0.04; 0.03]; v = sdpvar(nb, 1); l = sdpvar(nl, 1); Pij = sdpvar(nl, 1); Qij = sdpvar(nl, 1); Constraints = []; for k = 1:nl i = branch(k,1); j = branch(k,2); Constraints = [Constraints, ... v(j) == v(i) - 2*(r(k)*Pij(k) + x(k)*Qij(k)) + ... (r(k)^2 + x(k)^2)*l(k)]; Constraints = [Constraints, ... cone([2*Pij(k); 2*Qij(k); v(i)-l(k)], v(i)+l(k))]; end Constraints = [Constraints, v >= 0.95^2, v <= 1.05^2, l >= 0];

节点功率平衡我建议用关联矩阵一次性写，避免手写循环把方向搞混。核心逻辑是：流入节点的功率加上该节点注入（购电、分布式电源、CHP电出力），等于流出节点功率加该节点负荷。这里不展开完整代码，但提醒一句：branch矩阵里的方向定义必须和功率平衡公式保持一致，否则算出来的潮流符号会全反。

3.2 气网侧建模与压降方程

气网变量主要是节点压力平方alpha、管道流量f_pipe、气源注入量F_source和节点气负荷F_load。节点气负荷来自CHP的气耗和燃气锅炉的气耗。约束就那么几类：

• 气源点注入有上限，节点压力有上下限；
• 管道Weymouth约束用上一章的SOCP写法治；
• 压缩机如果存在，耗气量常简化成流量的一定比例，或者干脆用分段线性近似；
• 压缩机的升压比约束一般也是线性的。

这样做完，整个气网约束保持在SOCP体系内，不会引入非线性求解需求。

3.3 热网侧建模：双线性项的两个务实降阶思路

热网建模比电网、气网更“工程”一些。完整的动态热网模型涉及管道热延迟，稳态优化一般不考虑延迟，只考虑节点热功率平衡和温度混合关系。节点热功率约束是：

H_node = c_p * m * (T_s - T_r)

c_p取4.182 kJ/(kg·℃)，流量m单位kg/s，温差单位℃，算出来热功率单位是kW，再除以1000才是MW。这是单位最容易翻车的地方之一。

如果流量和温度都是变量，这就不再是SOCP。我在工程里最常用的做法是质调节：流量按设计工况固定下来，那么热功率约束变成线性约束，T_s、T_r是决策变量，H_node也是变量。温度还有上下限约束，比如供水温度不超过120℃，回水温度不低于40℃。管道温降的指数项，我会做成分段线性近似；节点温度混合约束本质是“流入该节点所有管道热量之和等于流出热量之和”，在流量固定以后同样线性化。

如果你不想固定流量，那就用外层迭代。第一次按设计流量解SOCP，拿出各支路流量后重新计算温降系数和混合系数，再解第二次。实测下来两三次迭代足够稳定，而且每轮都是凸问题，不会出现NLP那种中途发散的情况。

3.4 耦合设备可行域与目标函数

CHP是最需要认真建模的设备，它的电出力和热出力之间存在一个二维可行域。精确模型是复杂的非线性区域，工程上常用多边形近似，用一组线性不等式描述：

Constraints = [Constraints, P_chp >= Pmin, P_chp <= Pmax]; Constraints = [Constraints, H_chp >= Hmin, H_chp <= Hmax]; Constraints = [Constraints, P_chp + alpha_chp.*H_chp <= beta_chp];

alpha_chp、beta_chp来自可行域多边形顶点的拟合。加上启停变量后，变成：

u_chp = binvar(n_chp, 1); Constraints = [Constraints, Pmin*u_chp <= P_chp <= Pmax*u_chp]; Constraints = [Constraints, Hmin*u_chp <= H_chp <= Hmax*u_chp];

燃气锅炉的模型就简单很多：H_gb = eta_gb * F_gb，效率常数。电锅炉是H_eb = eta_eb * P_eb。P2G是F_p2g = eta_p2g * P_p2g。这些设备方程全部线性。

目标函数方面，除了前面说的购电购气成本，还要考虑设备启停和运维。这里提一个常见问题：目标函数里如果没有网损项或惩罚项，松弛锥可能“偷懒”，导致SOCP解虽然可行但不紧。解决办法后面专门讲，建模阶段先记住：目标函数最好加一个极小系数的网损惩罚项，比如1e-5 * sum(r.*l)，能显著改善松弛紧性。

4. MATLAB+Yalmip程序实现：变量、约束与求解接口

4.1 数据组织和决策变量声明

写这类程序，我强烈建议把所有网络参数装进struct或者独立脚本，不要在代码里到处写魔法数字。下面是一个单时段的示例，多时段只需把变量扩展成24列，约束按t=1:T循环即可。

% 统一单位：电功率MW，电压pu，温度℃，压力MPa % 电网参数 grid.nb = 5; grid.branch = [1 2; 2 3; 2 4; 4 5]; grid.r = [0.01; 0.02; 0.015; 0.01]; grid.x = [0.03; 0.05; 0.04; 0.03]; grid.P_load = [0; 0.2; 0.15; 0.25; 0.1]; % 决策变量 v = sdpvar(grid.nb, 1); l = sdpvar(size(grid.branch,1), 1); Pij = sdpvar(size(grid.branch,1), 1); Qij = sdpvar(size(grid.branch,1), 1); P_buy = sdpvar(1, 1); P_chp = sdpvar(1, 1); H_chp = sdpvar(1, 1); H_gb = sdpvar(1, 1); F_gas = sdpvar(1, 1);

变量名用Pij这种带方向的东西没问题，但要注意Yalmip里sdpvar默认就是全变量，不需要加'full'参数。有些老版本写法加了也不报错，但新版本可能把你当成矩阵对称去处理，坑过不少人。

4.2 电网DistFlow与SOCP约束组装

电网约束组装我在3.1已经给出了代码，这里补充节点功率平衡的写法。推荐用关联矩阵一次性构造：

% 关联矩阵 incidence(j,k)：支路k从节点j流出取1，流入节点j取-1 incidence = zeros(grid.nb, size(grid.branch,1)); for k = 1:size(grid.branch,1) incidence(grid.branch(k,1), k) = 1; % from incidence(grid.branch(k,2), k) = -1; % to end % 节点注入 = P_buy + P_chp（按实际接入位置映射） P_inj = zeros(grid.nb, 1); P_inj(1) = P_buy; % 节点1是上级电网连接点 P_inj(3) = P_chp; % CHP接在节点3 Constraints = [Constraints, incidence*Pij + P_inj == grid.P_load];

这里的符号约定是：incidence(j,k)=1表示支路功率从j流出，-1表示流入j。节点功率平衡写成“流出减流入加注入等于负荷”。实际项目里还要加无功平衡，方法和有功一致，把Qij和Q_load放进去。如果只做有功调度，无功可省略，但电压约束会失真，所以一般至少保留无功潮流。

4.3 气网、热网和设备约束的写法

气网的SOCP约束写法前面已经给了。热网在质调节假设下的约束写法如下：

% 质调节：流量固定为设计工况 m_fixed T_s = sdpvar(nh, 1); T_r = sdpvar(nh, 1); H_node = sdpvar(nh, 1); % 节点热功率 for n = 1:nh Constraints = [Constraints, ... H_node(n) == 4.182e-3 * m_fixed(n) * (T_s(n) - T_r(n))]; Constraints = [Constraints, T_s >= 70, T_s <= 120]; Constraints = [Constraints, T_r >= 40, T_r <= 70]; end

CHP可行域和启停约束、燃气锅炉和电锅炉的效率方程，按3.4的写法直接拼进Constraints。这里特别提醒：如果H_node是变量，并且它同时出现在热网节点平衡和CHP热出力约束里，不要重复声明两个不同变量去“相等”，直接在热网节点平衡里使用H_chp和H_gb本身，少一层中间变量就少一批数值问题。

% 热网节点热平衡：源节点注入等于负荷 Constraints = [Constraints, H_chp + H_gb == sum(H_load)];

4.4 求解设置与结果提取

模型组装完成后，求解设置是我的固定套路：

ops = sdpsettings('solver', 'cplex', 'verbose', 2); ops.cplex.param.mip.tolerances.mipgap = 1e-3; ops.cplex.param.timelimit = 3600; sol = optimize(Constraints, Objective, ops); if sol.problem == 0 fprintf('求解成功，总成本 %.2f 元\n', value(Objective)); Pij_opt = value(Pij); v_opt = value(v); l_opt = value(l); else disp(sol.info); end

value()是Yalmip提取变量数值的标准接口。多时段程序里，P_buy是24维向量，value(P_buy)直接就得到24个时段的购电功率曲线，画图非常方便。

5. 求解器选型与参数调优：从“能算”到“算得快”

5.1 为什么默认Cplex/Gurobi而不是Sedumi

Yalmip背后支持很多求解器，但如果综合能源模型里出现了binvar，问题就变成MISOCP。SDPT3和Sedumi这类纯内点法求解器只能处理连续SOCP，不认二进制变量。Cplex和Gurobi对MISOCP的支持非常成熟，分支定界、割平面、锥检测一条龙，基本不用你操心。如果你的问题没有整数变量，那么Mosek也很快，但从生态和接口稳定性来看，我还是优先推荐Cplex或Gurobi。

需要提醒的是，Yalmip默认会自动挑一个它能探测到的求解器。如果装了解算器但Yalmip没识别出来，先运行yalmiptest检查接口是否配好。很多“Solver not found”其实是mex接口没安装，不是模型问题。

5.2 数值缩放是最容易被低估的一步

这条经验值钱。电气热综合能源系统里，电压平方可能只有0.95，热网温度可能是一百多，天然气压力如果用了帕斯卡单位就是几百万，这几个数量级差到天上去了。Cplex内点法做锥约束检测的时候，对数值尺度极其敏感。变量量级差距超过1e6，常见的表现就是报Numerical problems，或者明明有可行解却判成不可行。

我踩过最惨的一次，就是天然气节点压力用了Pa做单位，模型怎么调都不可行，后来统一改成MPa，同一个模型秒解。从那以后，我的固定习惯是：电功率全部MW，电压全部pu，温度全部℃，压力全部MPa，成本全部万元或元，任何变量尽量压在1e-3到1e3这个区间。

5.3 常用求解器参数与中等规模算例表现

Cplex的常用参数：

ops = sdpsettings('solver', 'cplex', 'verbose', 2); ops.cplex.param.mip.tolerances.mipgap = 1e-3; % MIP间隙 ops.cplex.param.mip.tolerances.integrality = 1e-5; % 整数容差 ops.cplex.param.threads = 8; % 并行线程 ops.cplex.param.timelimit = 3600; % 时间上限

Gurobi对应的是：

ops = sdpsettings('solver', 'gurobi', 'verbose', 2); ops.gurobi.MIPGap = 1e-3; ops.gurobi.TimeLimit = 3600; ops.gurobi.NumericFocus = 2;

NumericFocus是我个人很喜欢的一个参数，数值有问题时调到2或3，求解器会更谨慎地处理矩阵计算，代价是速度慢一点，但总比不可行强。中等规模算例（电网100多个节点、气网30个节点、热网20个节点、24时段，含80个左右0-1变量），在我自己的测试环境下，纯连续SOCP一般十几秒能出结果，加整数变量后往往要到几分钟。如果你遇到几小时跑不完的情况，先别改参数，回头检查是不是给太多设备加了0-1变量——很多设备在调度周期内根本没有启停需求，固定成1就行，能砍掉大量分支。

6. 松弛紧性检验与排错实战

6.1 残差检查：判断二阶锥松弛是否“紧”

SOCP解出来后，不能直接拿去汇报。你要先判断松弛是不是紧的。判断方法很简单：检查每个支路是不是满足原来的等式约束。

res_max = 0; for k = 1:nl i = branch(k,1); j = branch(k,2); res_k = abs(v_opt(i)*l_opt(k) - (Pij_opt(k)^2 + Qij_opt(k)^2)); res_max = max(res_max, res_k); end fprintf('最大锥残差: %.3e\n', res_max);

在标幺制下，残差如果小于1e-5，基本可以认为松弛是紧的，解可以直接用。如果残差到了1e-2甚至更大，说明某些支路被松弛得很松，最优解虽然SOCP可行，但放到原物理方程里根本不对。最有效的补救办法：在目标函数里加一个很小系数的网损惩罚项1e-5 * sum(r.*l)。原因很朴素，目标函数里没有网损项时，潮流变量对模型来说是“免费的”，求解器怎么松快怎么来；加上网损惩罚后，它就有动力把l压在真实物理值附近。

气网也要做同样检查，看f_pipe^2和C^2*(alpha_i - alpha_j)的偏差。我见过不少项目电网上来就检查，气网完全忽略，最后结果里气网管道流量和压力的关系对不上，整个调度方案直接被否掉。

6.2 Infeasible和Numerical problems的排查链路

遇到不可行，我有一套固定的排查顺序，比盲目调参有效率得多：

第一步，把约束按模块拆开：只留电网约束，求解；再分别只留气网、热网约束，确认每个模块单独都是可行的。这个步骤能快速定位问题是出在哪一圈。

第二步，检查单位。这是最频繁的幕后黑手。尤其是热功率方程里那个4.182的系数，用kJ/(kg·℃)还是kW·h/(kg·℃)，结果差3600倍。压力用Pa还是MPa，又差1e6倍。一次不可行排查，先花十分钟把所有单位统一，往往就解决了。

第三步，检查负荷有没有“无源节点”。比如某个热负荷节点边上没有热源、没有管道连接，热量进不去，模型当然不可行。电网里如果出现孤立负荷节点，同理。

第四步，如果模型是MISOCP且报Numerical problems，把整数约束临时改成连续变量，看SOCP能不能解。连续能解、加整数就出问题，多半是数值缩放不良，参考5.2去处理。

6.3 常见报错速查表

下面这张表是从我自己的踩坑记录里整理出来的，按出现频率排序：

这里再补充一个Yalmip特有的坑：如果把sdpvar变量直接赋值给普通double矩阵，MATLAB会报类型错。遇到这种情况不是模型错，是数据结构设计问题，把存储容器改成cell或者直接重声明变量就好。

6.4 结果验证：与非线性模型对照的实用技巧

SOCP结果出来以后，我习惯做一次“物理后验”而不是直接信求解器。做法是：把求出来的P、Q、v代入原始非凸DistFlow方程，看电压降落方程是不是成立，潮流等式是不是吻合。这不是重新求解非凸问题，只是验证，速度快也不涉及初值。气网和热网同理。

如果项目允许，我会再用非线性求解器跑一个不要求全局最优的可行解来对照，看成本和SOCP结果差多少。如果SOCP成本远低于NLP可行解，大概率是SOCP松弛不紧，或者某些设备可行域被过度放大了。反过来，如果两者成本接近，那说明松弛质量很好，结果可以放心用。

还有一个容易被忽略的工程检查：看热网温度结果是否落在物理合理区间。比如某个时段热负荷为0，但求解器给出的CHP仍在供热，那基本可以断定热网模型少加了CHP最小热出力约束，或者热负荷平衡那里漏了用能下限。这种“数学可行、工程怪异”的解，靠残差是看不出来的，必须结合物理直觉去审视出力曲线。

最后说点个人习惯。我现在接到电气热综合能源优化程序，第一件事不是写代码，而是先画三张网络图，把所有设备的耦合关系列成一张表，统一好单位，再动手。很多“程序跑不起来”的坑，其实都是物理关系没理顺、单位不一致导致的。另外一个非常实用的小技巧：先跑纯SOCP不带整数的连续模型，确认收敛和残差指标合格之后，再加启停二进制变量。这样能很快区分是建模错误还是求解器数值问题。如果你也卡在“同样的模型，人家收敛我不收敛”这种问题上，先从这一步开始排查，大概率有奇效。