前一阵帮课题组调试一个电气综合能源系统的日前优化调度模型说实话第一次从零开始建这个模型的时候我心里是有点发怵的。原因倒不是电力和天然气网络本身复杂而是我怎么都绕不开那个让人头疼的“非线性”。一开始我直接按最常见的方式去写天然气管道用Weymouth方程电力潮流用交流潮流方程再加上机组出力约束整个模型直接变成一个混合整数非线性规划MINLP。结果不言而喻求解器一进去就开始分支定界半个小时过去还在原地打转连个可行解都算不上。后来我下决心把模型往二阶锥优化SOCP方向改造。整个过程花了不少时间也踩了不少坑但最终的效果确实值得模型在MATLABYALMIP环境下用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. MATLABYALMIP建模实战从变量定义到求解器配置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作为首选如果没装可以换成gurobioptions 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清清爽爽在几分钟内返回结果再去验证一下松弛间隙都接近零的时候你会觉得之前所有绕路都是值得的。