前一阵帮课题组调试一个电气综合能源系统的日前优化调度模型，说实话，第一次从零开始建这个模型的时候我心里是有点发怵的。原因倒不是电力和天然气网络本身复杂，而是我怎么都绕不开那个让人头疼的“非线性”。一开始我直接按最常见的方式去写：天然气管道用Weymouth方程，电力潮流用交流潮流方程，再加上机组出力约束，整个模型直接变成一个混合整数非线性规划（MINLP）。结果不言而喻，求解器一进去就开始分支定界，半个小时过去还在原地打转，连个可行解都算不上。
后来我下决心把模型往二阶锥优化（SOCP）方向改造。整个过程花了不少时间，也踩了不少坑，但最终的效果确实值得：模型在MATLAB+YALMIP环境下用Cplex求解，分钟级别就能出结果，而且我还能用松弛间隙去验证解是否真的满足原始方程。这篇就把我在这个过程中的建模思路、数学推导、代码实现和排错经验完整地写出来，给正在做电气综合能源系统优化调度的同学一个可以直接参考的实践路径。
1. 电气综合能源系统调度为什么绕不开二阶锥优化
1.1 我最初把它当成MINLP，然后被求解时间教训了一顿
先说清楚我在算的问题：一个典型的电气综合能源系统，电网侧有火电机组、负荷，气网侧有气源、气负荷，两者通过燃气轮机耦合。调度的目标是在满足负荷需求的前提下，让火电购电成本加上天然气采购成本的总和最小。看起来好像很简单，但真正把每个物理过程“翻译”成数学约束时，问题就变得复杂了。
电网上我一开始写的是交流潮流方程，节点电压幅值和相角互相嵌套，有功和无功也解耦不了。气网那边更麻烦，管道气流满足一个叫Weymouth的方程：管道流量与两端节点压力平方差的平方根成正比。这个方程天生是非线性的，而且不是那种好处理的“凸非线性”，是一上来就破坏整个模型凸性的那种。再加上燃气轮机把电、气两侧硬生生连在一起，模型里既有0/1变量（机组启停），又有一堆非线性项。
我把这套东西扔给求解器后，遇到了很多做这个方向的人都会遇到的窘境：可行性难找，找着了也证明不了最优；就算证明最优，一天时间也未必够用。我在那段时间里甚至怀疑是不是自己的公式写错了。后来回头一查，发现根本不是自己写错，而是问题本身的数学模型就不是一个“现代求解器友好”的模型。
1.2 模型里面真正让人头疼的非线性是哪些
要让问题变得可解，我最后做的选择是：对非关键部分做适度近似，对关键非线性做凸松弛，让整个模型变成一个二阶锥优化问题。先别急着觉得“近似”不靠谱，这套做法在很多电气综合能源系统调度论文里已经很常见了，核心是保证松弛在最优解处是紧的。
具体看，电力网络我直接用了直流潮流（DC Power Flow），也就是忽略无功、忽略电压幅值变化、假设相角差足够小。这样做的好处是潮流方程退化成线性方程，电网这一侧的约束全部线性化。直流潮流对输电网层面的经济调度来说精度是可以接受的，工程上都是这么用。
气网这一侧不能直接线性化，因为Weymouth方程的非线性是管道气流最基本的物理规律，硬性线性化要么误差大，要么需要引入一大堆分段线性变量，麻烦不说，效果也未必好。我采用的方案是把Weymouth方程改写成一个二次等式约束，然后把它松弛成不等式的旋转锥约束。关键就在于“松弛”之后，模型的可行域从非凸变成了凸，二阶锥优化问题就诞生了。
这样改造之后，整体模型就变成了一个包含线性约束、二阶锥约束和连续变量的凸优化问题。没有整数变量、没有非凸约束，Cplex、Gurobi这类求解器在求解这类问题上已经非常成熟，速度和稳定性都远非MINLP能比。
2. Weymouth方程的二阶锥松弛：数学推导与适用边界
2.1 从气流方程到旋转锥约束
先看Weymouth方程最常见的写法：
F_pq = C_pq * sign(p_p - p_q) * sqrt(p_p² - p_q²)
这里F_pq是管道流量，p_p和p_q是管道两端节点的压力，C_pq是管道常数。这个方程的意义很直观：压力差越大，气流越快，但两者是开方关系，不是线性关系。麻烦的地方在于流量符号由压力差方向决定，绝对值再乘C，整条式子既非凸也非凹。
处理第一步，定义变量Π_p = p_p²，把压力平方当作独立变量。这样Weymouth方程写成：
F_pq² = C_pq² * (Π_p - Π_q)
其中约定Π_p ≥ Π_q，也就是气流方向从p流向q。这里注意，等式约束本身依然是二次的，而且是非凸的。接下来做关键一步：把等式约束松弛成不等式：
F_pq² ≤ C_pq² * (Π_p - Π_q)
这一步的意义在于把可行域从一条曲线放宽到一整块区域，而这块区域恰好是凸的。为了理解它是凸的，可以把它写成标准的旋转锥形式。假设ΔΠ = Π_p - Π_q ≥ 0，那么：
4F_pq² + (ΔΠ - 1)² ≤ (ΔΠ + 1)²
展开后左边是4F_pq² + ΔΠ² - 2ΔΠ + 1，右边是ΔΠ² + 2ΔΠ + 1，整理得到4F_pq² ≤ 4ΔΠ，也就是F_pq² ≤ ΔΠ。所以原来的不等式就等价于一个标准二阶锥约束：
|| [2F_pq; ΔΠ - 1] || ≤ ΔΠ + 1
这个写法非常重要，因为很多求解器对“二次式小于等于线性式”这种形式不一定能自动识别为锥结构，但如果按照上面这种标准二阶锥形式写，Cplex和Gurobi都能直接处理，求解效率完全不一样。
2.2 松弛“紧”到底是什么意思，怎么判断
很多人第一次接触松弛问题会担心一件事：把等式放宽成不等式之后，求出来的最优解还满足原来的Weymouth方程吗？答案取决于松弛在最优解处是不是“紧”的，也就是最优解处那个不等式到底取等号还是严格小于号。
为什么通常会取等号？我的理解可以这样想：系统调度希望成本最小化，而天然气从低压端往高压端“凭空流”是不可能的，物理上管道流量必须靠压力差驱动。在目标函数最小化且存在供需平衡约束的情况下，最优解往往不会浪费压力差。直观说就是，能少买气就少买气，流量和压力差会不自觉地被推到Weymouth方程对应的边界上。
但“往往”不等于“一定”，所以我强烈建议每跑完一次模型都做松弛间隙校验：逐条管道计算 value(F_pq²) - value(C_pq² * (Π_p - Π_q)) ，理论上该值不应大于零（因为约束是≤），如果它近似等于0（比如1e-4以内），说明松弛是紧的，解对原始方程有意义。我在求解部分会给出这段校验代码，整个流程里它应该和结果输出放在同等重要的位置。
2.3 为什么电力网络也能顺手享受SOCP的好处
虽然我自己的实现里电网用的是直流潮流线性约束，但如果你想做得更精细，SOCP在电网侧同样有很强的作用，尤其是配电网或者某些需要计及电压幅值的场景。
配电网潮流里最常用的是DistFlow模型，节点电压幅值、有功/无功注入、线路功率之间存在耦合关系。把电压平方项单独设成一个变量，并且把原本的非凸二次等式约束松弛成不等式，得到的就是一个典型二阶锥松弛。这个松弛在很多辐射状配电网场景下被证明是紧的，这也是近年来一大批配电网优化论文能跑得动的原因。
所以你会发现，二阶锥优化在电气综合能源系统里的意义是全方位的：电网侧的交流潮流可以凸松弛，气网侧的Weymouth也可以凸松弛，两边都变成了同一种能用成熟求解器高效处理的模型。这是我后来觉得这个方向特别值得投入时间研究的主要原因之一，模型表达能力和求解效率被同时保证了。
3. MATLAB+YALMIP建模实战：从变量定义到求解器配置
3.1 为什么用YALMIP，而不是手写内点法
我记得第一次看SOCP代码的时候，心里也闪过一个念头：自己写一个内点法求解器是不是更酷？这个念头在我实际算了一个小时之后就彻底消失了。手写内点法要考虑的问题实在太多：锥约束的投影、迭代步长、收敛判据、数值稳定性，这些都是专门的数值优化领域工作，不是普通调度研究者该重复造轮子的事情。
YALMIP是一个MATLAB里的建模语言，它的价值在于把优化模型的“数学表达”和“求解器调用”分开。你用符号变量定义目标函数和约束，然后YALMIP自动判断模型类型，把问题转成求解器需要的格式。尤其是二阶锥约束，YALMIP可以直接用cone或者rotate这样的指令构造，最后传给Cplex或Gurobi求解，整个过程清爽很多。
唯一要记住的是，YALMIP本身不是一个求解器，它只是一个“翻译器”。真正算的时候还是靠底层求解器，所以装好Cplex或者Gurobi这一步不能省。如果没有商用求解器授权，用开源的Sedumi或者SDPT3也行，但速度会慢一些，我建议还是申请一份学术版Cplex。
3.2 数据准备和单位约定
建模之前先把单位统一好，否则后面到处都是坑。我的做法是电网功率全部用MW，所有负荷、机组出力、燃气轮机出力都是MW。气网所有流量统一用能量单位，不直接使用体积流量，具体换算方法是体积流量乘以天然气热值，得到对应的能量流（MW）。这样做的好处是燃气轮机耦合方程可以直接写成：P_GT = η_GT * S_GT，不用再乘热值密度之类的系数，代码逻辑清爽很多。
我用的案例系统是这样设置的：电力网络6节点，三台火电机组，节点4/5/6带负荷；天然气网络6节点，两个气源，一个燃气轮机挂在气网节点4，同时向电网节点3提供电出力。管道常数C_pipe我做了标幺化处理，压力变量直接用压力平方的标幺值，单位问题后面的坑专门讲。
下面是天然气系统的主要数据定义：
% 天然气网络：6节点，6条管道（双向，建模时预先固定方向） ng = 6; pipe_edges = [1 2; 1 4; 2 3; 4 5; 5 6; 3 6]; Cpipe = [0.8; 0.8; 0.7; 0.9; 0.7; 0.8]; % 管道常数（已标幺化） Pi_min = [9; 9; 9; 7; 7; 7]; % 压力平方下限 Pi_max = [16; 16; 16; 16; 16; 16]; % 压力平方上限 gas_load = [0; 0; 0; 5; 4; 4]; % 各节点气负荷（能量单位MW）注意这里gas_load是我随手填的算例数据，你换了真实系统之后要替换成自己算例的数值。
3.3 核心约束的代码写法
天然气侧最关键的就是Weymouth二阶锥约束。我在代码里没有直接用F_pipe² <= C² * (Pi_i - Pi_j)，而是写成标准二阶锥形式，目的就是确保求解器能识别成锥约束，而不是当作一个普通的二次约束，避免求解器误判导致性能下降甚至不可解。
Pip = sdpvar(ng, 1); % 节点压力平方变量 Fpipe = sdpvar(length(pipe_edges), 1); % 管道流量变量，非负 Sgt = sdpvar(1, 1); % 燃气轮机耗气量 Constraints = []; % 压力上下限 + 非负约束 Constraints = [Constraints, Pi_min <= Pip <= Pi_max]; Constraints = [Constraints, Fpipe >= 0, Sgt >= 0]; % Weymouth二阶锥松弛 for k = 1:length(pipe_edges) i = pipe_edges(k, 1); j = pipe_edges(k, 2); delta_pi = Pip(i) - Pip(j); % 标准旋转锥： || [2F; delta_pi - 1] || <= delta_pi + 1 Constraints = [Constraints, ... cone([2*Fpipe(k); delta_pi - 1], delta_pi + 1), ... delta_pi >= 0]; end接下来是节点气平衡约束。每个节点上，气源注入量加上管道流入减去管道流出，应该等于该节点的气负荷加上燃气轮机耗气量。我这里的管道边方向是预先固定的，也就是默认气流方向从from端到to端，实际算例中如果发现最优解里某条管道流量方向反了，就需要回过来调整建模方向。
% 燃气轮机所在气网节点 node_gt = 4; % 节点气平衡 for n = 1:ng in_flow = sum(Fpipe(pipe_edges(:, 2) == n)); out_flow = sum(Fpipe(pipe_edges(:, 1) == n)); balance = Ssource(n) + in_flow - out_flow - gas_load(n); if n == node_gt balance = balance - Sgt; end Constraints = [Constraints, balance == 0]; end这里我用了逻辑索引的写法，MATLAB里比较流畅，不会出现复杂的find循环。
电力侧的直流潮流约束相对简单。线路有功潮流等于电纳乘以两端相角差，节点上净注入等于净流出。注意我这里做了一个约定：节点注入为正，线路潮流从from端流出时为正。具体符号约定你可以按自己习惯调整，但所有节点必须保持一致。
nb = 6; theta = sdpvar(nb, 1); Pfire = sdpvar(3, 1); % 火电出力，分别在节点1/2/5 P_gt = sdpvar(1, 1); % 燃气轮机电出力，接在节点3 Pline = sdpvar(length(line_edges), 1); line_edges = [1 2; 1 5; 2 3; 2 6; 3 6; 4 5; 5 6]; Bline = [-20; -20; -20; -20; -20; -20; -20]; % 线路电纳 for k = 1:length(line_edges) i = line_edges(k, 1); j = line_edges(k, 2); Constraints = [Constraints, Pline(k) == Bline(k) * (theta(i) - theta(j))]; end Constraints = [Constraints, theta(1) == 0]; % 参考节点 load_e = [0; 0; 0; 40; 50; 60]; % 各节点电负荷 for n = 1:nb inj = -load_e(n); % 负荷作为负注入 if n == 1, inj = inj + Pfire(1); end if n == 2, inj = inj + Pfire(2); end if n == 5, inj = inj + Pfire(3); end if n == 3, inj = inj + P_gt; end net_out = sum(Pline(line_edges(:, 1) == n)) - sum(Pline(line_edges(:, 2) == n)); Constraints = [Constraints, inj - net_out == 0]; end燃气轮机的耦合约束写在最后，因为Sgt和P_gt分属不同网络，现在就通过效率关系把它们联系起来：
eta_gt = 0.45; Constraints = [Constraints, P_gt == eta_gt * Sgt];到这里，模型的约束部分就完整了。
3.4 目标函数和求解设置
目标函数是火电发电成本加天然气采购成本。火电成本用二次函数，天然气的成本我用线性函数，这在实际算例里已经是常见选择。
a = [0.02; 0.025; 0.03]; b = [30; 32; 35]; c = [100; 110; 120]; fire_cost = sum(a .* Pfire.^2 + b .* Pfire + c); gas_cost_coef = [4.2; 4.0; 4.3]; % 气源单位能量成本 gas_cost = gas_cost_coef' * Ssource; Objective = fire_cost + gas_cost;求解器设置我直接把Cplex作为首选，如果没装可以换成gurobi：
options = sdpsettings('solver', 'cplex', 'verbose', 2, 'showprogress', 1); optimize(Constraints, Objective, options);如果一切顺利，几分钟内就能得到最优解。
4. 求解完成之后，第一件事应该看什么
4.1 松弛间隙：你的解是不是真的满足原始方程
跑完optimize之后，很多人第一反应是直接value变量、画图、写报告，但我建议第一件事先做松弛间隙校验。原因很简单：二阶锥松弛只是在数学上等价于原约束的“放宽版”，如果不检查最优解处那些锥约束是否取等，你可能拿到一个在原问题里根本不成立的解还浑然不觉。
校验代码很简单：
relax_gap = zeros(length(pipe_edges), 1); for k = 1:length(pipe_edges) i = pipe_edges(k, 1); j = pipe_edges(k, 2); Fval = value(Fpipe(k)); Pival_i = value(Pip(i)); Pival_j = value(Pip(j)); lhs = Fval^2; rhs = Cpipe(k)^2 * (Pival_i - Pival_j); relax_gap(k) = lhs - rhs; end disp(relax_gap);因为我建模时写的是F² <= C²*(Π_i-Π_j)，所以lhs - rhs理论上只能小于等于0。如果某个数值大于0，说明你的模型定义或者求解器识别环节出了问题。如果等于0或者非常接近0（比如1e-4上下），说明最优解处Weymouth等式成立，松弛是紧的，这个解才是真正符合管道物理规律的解。我在实际算例里跑出来的结果，松弛间隙全部在1e-5到1e-6量级，基本可以断定松弛有效。
4.2 从调度结果倒推模型有没有写错
松弛间隙过了，代表数学上没问题，但物理上不一定是合理的。我习惯再检查几组输出：
第一，燃气轮机出力和耗气量是否满足效率关系。这个其实不用看，因为它是等式约束，强制满足。真正需要看的是燃气轮机落在了可行经济区间，还是顶到了上限边界。如果顶到上限，说明这个算例里气价相对电价的优势很大，系统倾向于用燃气轮机多发。
第二，节点压力是否都在上下限内，且管道两端的压力差是否与流量方向一致。比如节点i压力平方大于节点j，那么F_pipe(i,j)应该非负，如果算出来Fpipe为负，说明我预设的管道方向与实际最优流动方向相反，这个在模型计算里不会报错，但说明初始方向假设不对，需要调整后再审视结果。
第三，电网线路潮流有没有越限。直流潮流约束本身只有等式约束，我没有加线路容量上下限，如果你想加，直接在Pline上加上限和下限即可，但要注意线路潮流有正负两个方向，上限和下限要对称设置，或者根据具体线路允许的单向流动来设置。
4.3 用敏感度分析验证耦合建模是否正确
还有一个我常用的验证手段：把气源价格系数稍微调高10%，看燃气轮机出力是否下降、火电出力是否上升。如果耦合模型正确，这种敏感性必然是反向的，否则说明模型里耦合关系写错或者目标函数单位不一致，此时排查起来就快捷多了。
敏感度分析本质上不需要重新推导，你只需要把gas_cost_coef改掉，重新optimize一次，记录P_gt的变化即可。这个方法虽然简单，但真的能暴露很多隐藏问题。我印象最深的一次，就是因为气网单位换算不对，导致敏感性分析结果完全反直觉，我才发现问题不在求解器，而在物理量纲上。
5. 我在实际跑代码中遇到的五个典型坑
5.1 求解器“并不认识”你的锥约束
我第一次用YALMIP写Weymouth约束时，图省事直接写成了Fpipe^2 <= Cpipe^2 * (Pip(i)-Pip(j))。语法没报错，Cplex也跑起来了，但求解速度慢得离谱，而且结果间隙很大。后来我把约束显式写成cone形式后，求解时间一下子从几十分钟降到几分钟。原因就是YALMIP对“二次≤线性”这种形式有时会当作非凸二次约束或通用约束处理，求解器没办法直接利用锥结构。
所以我强烈建议：凡是能用cone或者rotate写出来的SOCP约束，就不要偷懒写原始二次不等式形式。这个教训在我后来做其他SOCP问题时反复被验证。
5.2 压力平方的量纲害我浪费了一整天
天然气节点压力在工程单位里经常是MPa量级，比如上游压力4MPa，下游3.2MPa，平方之后一个16、一个10.24，看起来尚可。但如果你的压力单位是kPa，一个4000kPa的节点，压力平方就是1.6e7，这个数值直接放进约束里，和流量、负荷那些几十几百的数据放在一起，数值尺度会差出好几个数量级，轻则求解精度下降，重则直接报数值奇异。
解决办法很粗暴：所有压力做标幺化。先定一个基准压力p_base，然后把压力平方的变量定义成(p/p_base)²。这样变量范围就变成了0.8到1.2这种适中尺度，整个模型数值稳定性立刻改善。我在代码示例里填的Pi_min、Pi_max都是标幺后的值，就是这个意思。
5.3 双向管道的方向问题
我记得第一次让管道允许双向流动时，模型直接跑出了负流量，然后压力差方向的约束就冲突了。因为cone形式本身默认ΔΠ≥0，如果你让流量变量可正可负，就得额外引入方向选择变量和0/1变量，模型就从纯连续SOCP变成了混合整数SOCP，求解难度直接上一个台阶。
在工程实践中，天然气管道的方向在一个调度时段内通常是比较明确的，要么由上游气源到下游负荷，要么在环状管网里动态变化。我的建议是，第一版就先按“预知方向”建模，也就是说你根据节点压力等级预先假设每条管道从高压端流向低压端。这样模型简单，速度快，跑完之后再去检查最优解里的流量方向有没有和预设冲突。如果真有一条管道方向反复横跳，再去考虑用MISOCP处理那条管道。
5.4 燃气轮机耦合关系的单位换算
这个问题我前面已经提过，但值得再单独强调一下。电网用MW，气网如果用m³/h，那么耦合方程就得写P_gt = η_gt * GCV * S_gt，还得把气体的热值、密度、标准状态参数全部乘进来，稍不留神系数就少乘一个。我后来的做法是彻底统一到能量单位，也就是把管道流量、气源供气、气负荷包括燃气轮机耗气全部换算成MW再进模型，这样耦合方程就是单纯的线性效率关系。
具体的换算逻辑是：S_gt_mw = S_gt_flow * gas_heat_value / 3600，其中gas_heat_value取标准天然气的低位热值，单位要匹配。这些换算看着不起眼，但一旦算错，整个调度结果全都会变味。
5.5 Cplex报infeasible时的排查顺序
如果你把模型搭好之后Cplex直接报infeasible，先不要急着怀疑求解器或者YALMIP。我自己的排查顺序是从简单到复杂：
第一步，检查各个节点平衡约束是不是多了或少了一项。电力系统里最容易漏的是参考节点相角为0的约束，少了它潮流方程不唯一；天然气系统里最容易漏的是燃气轮机耗气那一项，少了它气源供给总是不平衡。
第二步，检查上下界是不是矛盾。比如压力平方范围设置太窄，或者管道常数Cpipe太小，导致可行域根本不存在。前一段我提到压力平方量纲问题，如果上下界是标幺值而流量没有标幺，两者差距过大，也会报无解。
第三步，把目标函数临时改成常数（比如0），只求解可行性问题。如果这样还是infeasible，那基本可以断定约束之间存在不等式冲突，拿一组具体数值去逐条对算一遍，很快就能找到问题。
这套排查顺序我后来成了习惯，比盲目加罚函数项有效得多。
6. 从静态SOCP向多时段、不确定性和电转气场景的扩展
6.1 多时段调度的变量追加
前面这套模型都是单一时段的调度，实际系统中调度员更关心未来24小时甚至更长时间序列的优化调度。扩展思路其实是线性增长：给每个时段复制一组决策变量，在时段之间加上储能约束或爬坡约束，再给目标函数加一个对时间求和的外层。
具体到SOCP，你只需要把Pip、Fpipe、Pfire这些变量从向量变成矩阵，每个时段对应一列。不同时段之间的耦合主要体现在机组爬坡约束和气网储气罐的动态上。这些耦合约束本身是线性的，所以整体依然是一个大型SOCP，只要矩阵稀疏性处理得当，Cplex跑起来依然很舒服。
6.2 天然气储能与管网动态
如果想把模型做得更贴近实际，天然气网络的动态过程不能完全忽略。管道本身存在“线包”效应，也就是管道可以储存一部分天然气，短期压力变化会导致暂态流量和稳态流量不完全一致。这种模型在数学上一般用偏微分方程描述，直接放进优化里很难。
实际工程中常见做法是把管道分段，用有限差分把偏微分方程离散成线性或双线性约束，有的文献干脆用储气罐节点来近似管道的储气能力。这套处理方法在SOCP框架下是可以衔接的，只要储气的容量和充放气速率约束保持线性，或者至多引入一些锥约束，模型依然可控。
6.3 P2G与碳排放约束
除了传统的“气到电”方向，现在很多综合能源系统还会考虑电转气（P2G）技术，也就是在新能源大发的时候，把多余电力通过电解水制氢，再掺混或甲烷化注入天然气系统。这个过程实际上给电、气两个网络增加了一条反向的耦合通道。
P2G的核心约束是进气流量与消耗电力之间的关系约束，以及产气注入气网的压力边界条件。前一个约束通常是线性效率关系，后一个在计及注入压力时需要满足气网节点的压力界，同样可以写成锥约束形式。另外，如果要考虑碳排放上限，只需要在目标函数或者约束上增加一个线性的碳排放流平衡，整个SOCP骨架不用动。
我自己对这套代码最大的体会是：把模型从MINLP改成SOCP之后，真正的精力和时间反而用到物理建模和结果验证上了，这其实才是做研究该有的状态。如果你正在被电气综合能源系统优化调度的求解速度折磨，我建议真的可以停下来花一两天时间把Weymouth方程的二阶锥松弛吃透，然后把代码重写一遍。等你看到Cplex清清爽爽在几分钟内返回结果，再去验证一下松弛间隙都接近零的时候，你会觉得之前所有绕路都是值得的。
