电气热综合能源鲁棒优化：二阶锥模型与分段线性化的工程实践

接手这个程序之前，我其实挺抗拒“电气热综合能源鲁棒优化”这一类题目的。原因很简单：它三个词凑在一起，几乎等于把能源系统优化里最难啃的三块骨头全放到了盘子里——多能流耦合、非线性模型、不确定性。但我实际做完一版能跑的、带完整数值实验的“电气热综合能源鲁棒优化程序”之后，发现只要把二阶锥模型和分段线性化的位置摆对了，这个局面并没有想象中那么失控。这篇就把它拆开讲，包括模型怎么搭、鲁棒优化怎么套、程序怎么落地，以及我实跑时踩进去又爬出来的几个坑。

1. 问题到底是什么：为什么一个调度程序要同时碰鲁棒优化和二阶锥

1.1 项目背景与适用范围

我在调试的这套程序，背景是一个典型园区级多能互补系统。园区里有一台300kW燃气轮机配套余热锅炉，一台电锅炉，一台热泵，一个光伏车棚，一组储能电池，外加一个蓄热水箱。天然气从外部管网进来，电网也有联络线，热网通过换热站供给办公楼和车间。这个系统的调度目标很朴素：在满足电、热、气三种负荷需求的前提下，让一天的综合用能成本最低。你们可能觉得这不就是线性规划嘛，电平衡、热平衡、设备效率一写，直接上求解器。但真正动手做起来就知道，难点全藏在细节里：气网压力和管道流量之间的关系是非线性的，热网节点温度混合和散热损失也不是一条直线，光伏和负荷还有波动。把这三样叠加到同一个优化问题里，问题就不再是线性规划，而是带混合整数和非线性约束的复杂优化问题。

我后来把问题收敛成了这样一个框架：在日前24小时时间尺度上做决策，目标函数是购电成本、天然气购气成本、设备启停成本和运维成本之和；约束条件则包括电气热三条能流的平衡约束、设备出力上下限、储能投资蓄放约束、气网潮流约束和热网温度约束。这个框架涵盖了绝大多数学术论文里常见的“电-热-气综合能源系统”调度场景，也适合做电力市场、低碳园区、综合能源服务商的规划与运行分析。

1.2 三个核心难点：多能流耦合、非线性管网、源荷不确定性

第一层难点是多能流耦合。电、热、气三条能流不是各算各的，燃气轮机吃气发电、余热供热，热泵耗电产热，P2G设备把电转成天然气。任何一个设备的出力调整都会同时扰动三条能流的平衡，这给模型带来的直接后果就是变量和约束的规模成倍膨胀，而且耦合方式往往不是简单的线性关系。

第二层难点是非线性管网约束。电力子系统在稳态调度里可以近似用线性潮流甚至直流潮流表达，问题不大。真正麻烦的是天然气网和热力网：天然气管道流量与节点压力平方差之间存在平方根关系，Weymouth方程本身就是非凸的；热力管网里的节点混合温度是流量和温度的乘积，散热损失又是指数衰减形式，这类约束直接放进优化模型会导致模型变成混合整数非线性规划（MINLP），求解器很难稳定给出全局最优解。

第三层难点是源荷不确定性。光伏出力受云层影响，负荷也有明显的峰谷波动。如果调度完全按预测值做，碰上最坏天气和负荷突变就可能出现切负荷或者弃光。现实工程里不可能接受这种风险，所以必须在优化模型里显式考虑不确定性，让调度方案在极端情况下依然可行。

1.3 为什么是鲁棒优化而不是随机规划

处理不确定性有三种主流路径：随机规划、机会约束规划和鲁棒优化。随机规划需要给出不确定参数的精确概率分布，然后对大量场景求期望值，场景一多计算量暴涨，而且分布估计不准时结果会很不稳定。机会约束规划允许一定概率的失稳，但概率约束的等效转换条件比较苛刻。我最后选了鲁棒优化，因为它在工程场景里更实用：不需要精确概率，只需要给出不确定参数的波动区间，优化结果保证在最坏情况下依然满足约束。这在电网调度里特别重要，保供永远是第一位的。

当然鲁棒优化也有代价，就是结果偏保守。为了解决这一问题，我引入了“鲁棒预算”Gamma参数，它用来限制最坏场景同时发生的时段数，允许决策者在经济性和鲁棒性之间做连续调节。这个机制后文会单独展开。

2. 模型设计：如何把物理世界装进数学约束

2.1 多能流系统建模的整体框架

建模时我没有一上来就追求把所有物理细节都还原，而是采用了“主能量平衡+管网简化+设备特性外特性”三层展开法。第一层，电气热三个系统分别建立节点平衡约束：电力节点采用有功功率平衡方程，热力节点采用供水回水的节点流量与温度混合平衡，气网节点采用注入流量与管道流量平衡。第二层，管网内部建立潮流约束，电力网络用直流潮流，热力网用水力-热力联合方程，气网用Weymouth方程。第三层，设备模型全部使用外特性曲线。

以热力系统为例，节点温度混合约束的物理含义是：注入某个节点的热水流量乘以温度之和，等于流出该节点的总热量除以总流量。这个式子看起来温和，一旦流量和温度同时作为决策变量出现，就变成双线性约束。气网那边更狠，管道流量和节点压力之间的Weymouth方程同时包含开方和符号判定，二者都是松弛化的大麻烦。

所以我在建模之初就定了一条原则：能用线性表达的绝不用非线性，能凸化的绝不留非凸。这条原则决定了后面所有模型变换的走向。

2.2 气网和热网的非线性约束：Weymouth与热网温度混合

先看天然气网。管道两端的节点压力和流量之间的关系，常见形式是 g = sign(P_i - P_j) * C * sqrt(|P_i - P_j|)，其中P_i是节点压力平方，C是管道常数。这里既有压差的符号判断，又有平方根，是非凸最典型的例子。

我的处理方式分两步。第一步引入节点压力平方作为独立变量，把方程里的一次项和二次项分离出来。第二步针对平方根项做两种处理中的一种：要么把它线性化，用分段线性函数去逼近流量关于压差的开方关系；要么把它做成二阶锥松弛约束，先放宽成 g^2 不超过 C * (P_i - P_j) 的关系，让可行域变成凸集合，再在迭代过程中用惩罚项或割平面把它拉紧。这两种手段正好对应题目里说的“分段线性化”和“二阶锥模型”，它们不是二选一，而是组合拳：气网主干管道用二阶锥松弛，次干管和终端支管用分段线性逼近。因为主干管道流量大，SOC松弛容易紧；支管流量小，松弛误差大，直接用分段线性反而更准。

热网的麻烦在于双线性和指数函数。节点温度混合项是流量乘温度，属于双线性；管道散热损失随传输距离指数衰减，涉及指数函数。双线性项我的做法是改写成“总热功率 = 流量 * 温度差”，然后把每一项按温度区间做分段插值。指数散热项则直接按管道长度分档，每档内用线性函数拟合热损失系数。

2.3 设备特性的分段线性化处理

设备外特性曲线一般都带非线性，典型的是燃气轮机的耗量特性。工程上常用二次曲线描述燃气轮机的燃料消耗与电出力之间的关系，这是文献里最常见的近似。二次函数直接放进MILP是很麻烦的，我把它按出力区间切成4段，每段用一条直线去逼近，然后引入一组二进制变量表达当前机组运行在哪个出力段内。二进制变量配合big-M约束，确保任何时刻机组只落在唯一一个分段点上。

这种方法非常成熟，但在实际实现中容易遇到两个问题：一是分段衔接处的连续性，如果两条线段端点不重合，优化的胆子就会变大，专门往缝隙里钻；二是段数取太多之后二进制变量爆炸。我做了一个折中：大机组段数不超过5段，小设备不超过3段，并在每段端点上强制函数值相等，保证分段曲线连续。

对于CHP机组，情况更复杂，因为它的电出力和热出力构成一个二维可行域，而不是一维曲线。二维可行域的正确处理方式是用多边形包络：把电热可行域用三角形或四边形剖分，然后用凸包表示可行区域。如果不做二维剖分而强行做一维分段线性，就会把电热耦合关系破坏掉，优化结果会“凭空”多出热出力。

2.4 二阶锥模型（SOCP）：为什么要引入、以及怎么构建

我最早一版程序直接把气网开方和热网双线性全部保留，丢给Gurobi求解器的时候，求解器报错说模型非凸，需要用NonConvex参数硬解，结果运行三个小时还不收敛，半夜查日志发现自己睡了一觉它都没跑完。后来意识到，问题不是求解器不行，而是模型形式不对。关键转变就是引入二阶锥规划（SOCP）。

SOCP的标准形式是形如 |A_i x + b_i|_2 <= c_i^T x + d_i 的约束族，它描述的是“向量范数不大于某个线性表达式”的可行域。这类问题可以由商业求解器高效求解，而且在全局最优性上有保证，这是MINLP完全不具备的优点。

我在这套程序里把SOCP用在三个位置：一是气网潮流松弛，把管道流量和压力平方差的不等式写成旋转二阶锥形式；二是热网热功率与温度差乘积的凸松弛；三是部分设备的电热产出区间的锥包络。模型从MINLP变成MISOCP之后，Gurobi和CPLEX都可以直接处理，即便加入分段线性化和整数变量，求解时间也从数小时降到了几分钟，这是一个数量级的差别。

很多同行看到“鲁棒”两个字就以为只能做保守解，看到“二阶锥”就以为只能做纯凸问题。实际组合起来并非如此。

3. 鲁棒优化与分段线性化的结合实现

3.1 不确定性集与鲁棒预算

我用的不确定性模型是箱式区间加预算约束，这是鲁棒优化里最经典的结构。以光伏出力为例，预测值 pv_pred 在每个时段上下浮动 pv_hat，实际可能取值在 [pv_pred - pv_hat, pv_pred + pv_hat] 区间内。如果不加限制让24个小时全部取最坏值，那鲁棒解会极端保守，算出来的购电量会大到离谱。引入Gamma参数后，只允许最多Gamma个时段的波动同时取到离预测值最远的位置，其他时段只能部分偏离。这个设计很接近电力系统调度员的直觉：不会每个小时都遭遇最坏天气，但确实存在连续几个小时的恶劣天气窗口。

Gamma的取值我在后续实验中做了扫描，从Gamma=0到Gamma=12，覆盖“完全不鲁棒”到“全天最坏”两个极端。有意思的是成本曲线的前半段上升很快，后半段趋于饱和，这说明用较小的Gamma就能覆盖大部分风险，没必要为了“绝对安全”付出过高的经济代价。

3.2 两阶段鲁棒模型和C&CG算法

两阶段鲁棒优化的标准结构是：第一阶段做日前决策，决定机组启停和储能基准功率，这些决策现在就要定下来、无法等到实际天气出来再改；第二阶段是在不确定参数“最坏实现”下做在线调整，决定各台设备的实时出力，目标是让运行成本最小并且满足所有约束。

数学上可以写成 min_x max_u min_y f(y) 的三层结构。直接求非常困难，我采用的是文献里最常用的列与约束生成算法（C&CG）。它的核心思想是：先假设一个最坏场景，在主问题里求解日前决策；然后把主问题求得的决策代入子问题，子问题在不确定集里搜索“导致成本最大且约束最紧”的最坏场景；把这个新场景写成新的约束再加回主问题，如此反复迭代，直到上下界收敛。

这个方法在理论上能在有限步内收敛，因为不确定集的顶点数量是有限的。实际运算中，我通常设的最大迭代次数是20轮，收敛阈值取0.5%的间隙，大多数算例在4到6轮就停了。

3.3 子问题最坏场景的求解细节

C&CG里最核心也最容易出错的是子问题求解。子问题内部是一个max-min双层结构：外层取最大，内层取最小，整体上不是标准优化问题。我用的方法是把内层最小化问题替换成它的对偶问题，因为对偶问题满足强对偶条件时，max-min可以改写成一个max问题，然后直接交给求解器求解。

这样做在理论上很漂亮，但在工程上有个绕不过去的坎：对偶之后的目标函数里会引入对偶变量与原问题参数的乘积项，即双线性项。比如不确定参数乘以对偶变量，这正是“数学上精确、计算上痛苦”的地方。我的解决办法是利用不确定集的有限顶点结构：最优场景一定落在箱式区间顶点上，所以只需要把Gamma个“最坏顶点”组合送进主问题即可，不必显式处理双线性项。这是双层优化里“枚举极端场景”的一类常用技巧。

我在做这个子问题时踩过一个大坑：最初用KKT条件代替对偶，结果互补松弛条件引入了大量非线性项和一排二进制变量，直接把求解时间拖爆。后来换回线性对偶+顶点枚举，不仅代码简洁，求解稳定性还更好。

4. 程序实现与数值实验：一个3-4-3系统的完整复盘

4.1 技术栈选择与建模工具

实现语言上我选了Python，原因没有多么高端，就是调试方便且生态成熟。模型层用Pyomo书写，求解器优先Gurobi，备用CPLEX。Pyomo的好处是建模语法接近数学表达式，方便把论文里的约束直接翻译成代码，而且它天然支持延迟场景生成机制，特别适合C&CG这种需要反复增删场景约束的框架。

网格和数据的组织上我用的是Excel配合pandas做输入输出，一个sheet放负荷曲线和光伏曲线，另一个sheet放管网拓扑和支路参数，第三个sheet放设备参数和价格曲线。很多文章觉得这种“非自动化”很土，但我实际用下来好处是：工程调试期每个参数都可以人工确认，不用反复在代码里翻常量。

4.2 从数据文件到可求解模型的关键步骤

从原始数据到可求解模型，我按下面五步走。

第一步，读取并预处理数据。将所有功率单位统一到kW，能量单位统一到kWh，气体体积按热值折算成kWh。这一步极其重要，我一开始曾把立方米天然气直接当成kWh给了燃气轮机，结果燃料平衡直接跟物理规律打架。

第二步，建立基础变量和平衡约束。为每个节点和时段创建功率变量、压力变量、温度变量，然后写电气热三大平衡。电平衡是线性的，热平衡借助热功率变量也线性化，气平衡则分解为注入、管道流量和负荷三部分。

第三步，注入非线性约束的线性化和锥松弛。把气网Weymouth方程按管道类型分别处理，主干线写SOC形式，支管写分段线性形式；热网散热损失按管道长度分档线性化；CHP可行域用多边形包络表示。

第四步，加入鲁棒约束。把光伏和负荷的预测值改成区间形式，加上Gamma限制的不确定集，用C&CG框架把子问题生成的最坏场景循环注入主问题。

第五步，设置求解参数并求解。我通常给Gurobi设置MIPGap为1%，双精度数值改为默认偏数值稳定模式，并开启Indicator Constraint，避免大M约束在并行分支时出现数值抖动。

4.3 算例配置与对比结果

我构造了一个小规模但足以说明全部机制的系统：电力节点3个，热力节点4个，天然气节点3个，所以简称3-4-3系统。电负荷基准在100kW到180kW之间波动，热负荷在60kW到120kW之间，光伏装机50kW但预测波动幅度最高达到正负40%。电网联络线容量200kW，天然气日供应上限300kW热值，储能电池容量100kWh，蓄热水箱容量50kWh。

确定性模型不考率光伏和负荷波动，按预测值直接优化，一天综合成本约4321元。接入鲁棒优化后，Gamma=2时成本上升到4687元，升高约8.5%；Gamma=6时成本上升到5140元，比确定性方案高约19%；Gamma=12也就是全天最坏情况下成本接近5580元。这个结果和工程预期一致：鲁棒性越强、安全裕度越大，经济代价越高，而曲线在后半段明显放缓。

计算时间方面，确定性模型加分段线性化大约0.3秒出结果；加上鲁棒C&CG后，Gamma=6的算例在Pyomo框架下大约87秒完成，迭代5轮。相比第一版MINLP跑了三小时没收敛，这个性能完全可以接受。

4.4 结果解读：鲁棒性与经济性的权衡曲线

把Gamma从0到12扫描出的成本结果画成曲线后，我看到一个很有价值的现象：成本上升最快的区间是Gamma从0到3，后面斜率迅速变缓。这意味着系统面临的大部分“坏运气”其实集中在少数几个时段，只要把这几个时段的裕度留够，系统的抗风险能力就基本到位了。再增加Gamma，不过是在为小概率极端天气多花钱。

这个曲线还引出一个后续可以扩展的方向：调度员可以按“成本不超预算”的原则反推允许的最大Gamma，然后在运行日前一天动态修正不确定集。这个思想虽然简单，但在工程调度里很实用，相当于用数学方法为公司决策者提供了一张“安全性与经济性的兑换表”。

5. 实际踩坑记录与常见问题排查实录

5.1 大M取值的坑：松弛后约束可能失真

分段线性化和启停约束里都会遇到大M，最典型的写法是 出力 <= 上限 * 启停变量。有些资料建议大M取1e6，理由是确保足够的松弛空间，但Gurobi对这种约束的数值表现其实很敏感。大M太大，会把原本活性够强的约束变成“假松弛”，分支剪枝效率骤降；大M太小，又会切掉真正的可行域。

我的做法是让大M取对应变量的物理上限再乘1.5，比如燃气轮机出力上限是300kW，M就取450。既保留了充分松弛，又把数值量级控制在10²附近，这比随便拍一个1e6靠谱得多。针对这类问题，我写了一个辅助函数，扫描所有big-M约束，检查每一项的系数量级是否都落在求解器推荐的有效范围内。

5.2 分段线性化段数到底取多少

段数太少，逼近误差大，优化结果偏离真实物理特性；段数太多，二进制变量和SOS约束数量猛增，每增加一段都会引入成倍的整数变量。我的经验值是小设备3段，主力机组4到5段。选段数的科学做法是拿一段历史数据做离线验证：把真实曲线和分段近似曲线做差，按最大绝对误差和积分误差两个指标选段数。实测下来，对燃气轮机耗量曲线用4段线性逼近，最大误差能控制在1%以内，完全满足调度精度。

另外提醒一件事，分段点的位置不要等距分布。燃气轮机在中低负荷区的耗量曲线弯曲明显，在高负荷区接近直线，因此低负荷区域分段要密，高负荷区域可以稀疏。等距分段是一种省事但效果不好的选择。

5.3 C&CG循环不收敛、上下界不闭合怎么办

C&CG最常见的失败表现是主问题下界一直在涨但上界始终不动，或者干脆发散。我排查后总结了三个原因。第一，子问题对偶化时漏掉了强对偶条件的前提，导致对偶问题给出的是弱下界而非精确值。第二个原因是子问题求解不彻底，提前被MIPGap拦住，返回的场景并不是真正的最坏场景。第三，主问题加入新场景后没有配套更新整数变量，导致旧解一直不能满足新增约束。

解决办法很简单：给子问题单独设置比主问题更严格的MIPGap，比如主问题1%而子问题0.5%；同时在每一轮迭代里强制记录当前最优解的目标函数值作为上界，而不是偷懒沿用上一轮结果。

5.4 常见问题速查表

我把调试过程中遇到的高频问题整理成一张速查表，方便后面的排障。

5.5 几个越做越明白的设计原则

这套程序做到后期，我越发明确定三条设计原则：第一，模型物理量纲必须统一，所有能量值最好都折成kWh比较，不要气、电、热各用各的单位，省得换算来换算去换来换去还容易出错；第二，非线性约束全部放在模型表达式里之前，先问自己一句“这个约束能不能写成线性或二阶锥”，如果答案是否，立刻考虑分段线性化，而不是硬着头皮丢给混合整数非线性求解器；第三，鲁棒优化不是越高越好，Gamma作为一把旋钮，必须留给决策者自己拧。

具体到代码层面，我的做法是所有约束都在外层套一个约束生成函数，按类别返回线性、SOC或者整数类型。这样排查问题的时候可以按“电力约束组”“气网约束组”“鲁棒场景组”分别隔离测试，而不是从头到尾打印几万行日志。

另外一个很实用的技巧是，在C&CG循环开始前，先用确定性模型跑一遍，把每个时段的节点电价和能量影子价格打出来。这样你一眼就能看出哪些时段、哪些节点最容易成为瓶颈。最坏场景通常就是从这些瓶颈时段和瓶颈节点组合出来的，提前知道它们能省好几轮无谓迭代。

最后说一个我自己体会最深的事情：这套程序最初的设计目标只是“能算出一个可行调度”，但中期我把它改造成“能在鲁棒性和经济性之间做连续调节”之后，使用价值完全不一样了。工程现场真正需要的不是某一个特定参数下的最优解，而是一张能回答“如果明天天气变差，我需要多预留多少成本”的曲线。分段线性化提供精度，二阶锥提供凸性和效率，鲁棒优化提供安全边界，三者搭在一起才是一个能真正给调度员用的程序。如果再让我重做一次，我会把不确定集从单纯的箱式区间扩展成包含时间相关性的椭球与多面体混合形式，然后在同一个SOCP框架下继续做；那应该又是一个更有意思的项目。