1. 为什么P2G和碳捕集要放进同一个模型问题的起点与耦合逻辑做综合能源系统优化调度的人早晚都会遇到这么一个问题电热联供机组CHP一开热是满足了但电出力跟着被绑死风电光伏大发的时候电又多到没处去只能弃风弃光碳排放又随着煤/气消耗一路上涨。单独加一套电转气P2G设备能把多余的电变成天然气存起来但P2G制甲烷恰恰需要二氧化碳作为原料。单独加一套碳捕集设备能降低火电/CHP的碳排放但捕集本身要消耗大量电能和热能又反过来挤压系统的供电能力。两个设备单看都有价值放在一起却产生了天然的双向互补关系——碳捕集提供CO2P2G消耗CO2P2G消耗电能碳捕集也消耗电能但同时降低了系统碳排放。我最初做这个课题时想法很朴素在传统的热电联供综合能源系统里把电储能、热储能、气储能都加上再引入P2G设备和碳捕集设备然后看它们之间如何协同调度。但真正动手建模、写MATLAB代码、调结果之后才发现这个模型的核心难点根本不在“加设备”而在于三个层面一是两种新设备与CHP机组之间形成了复杂的能量-物质耦合闭环二是碳捕集设备的能耗特性会显著改变系统电热平衡关系三是这些设备的运行约束直接影响了日前调度的可行域处理不好就会出现一堆匪夷所思的劣质结果。这篇博文就把我实际搭建这个调度模型的过程、踩过的坑、MATLAB代码的组织方式以及一些值得注意的建模细节原原本本整理一遍。适合正在做综合能源系统优化调度、碳捕集与P2G耦合建模的同行参考也适合刚接触这类模型、想用MATLAB快速搭一个可运行版本的朋友做地图式导览。我不会把代码逐行粘贴因为那会变成一份没有灵魂的代码说明书我会重点讲清楚“每一部分代码背后的建模决策是什么为什么这么写换了场景要改哪里”。2. 模型骨架目标函数与约束条件的数学化拆解2.1 目标函数的费用构成别只盯着购能成本典型的日前优化调度目标函数一般是总运行成本最小化。我实际使用的目标函数包含四块购能成本购电、购气、碳交易成本、设备运维成本、弃风弃光惩罚成本。之所以单独列出碳交易成本是因为碳捕集设备的收益必须在经济性上体现出来如果不能把“捕了多少碳”换算成钱模型就会本能地“不想开”碳捕集设备。购电成本按分时电价计算这个比较简单——逐时电价向量点乘逐时购电量。购气成本类似按天然气价格计算注意CHP耗气和P2G产气都要在气平衡中体现。运维成本我习惯按出力比例计算单位是“元/MWh”而不是固定启停费用因为对于日前连续调度来说边际运维成本对调度结果的影响更平滑。弃风弃光惩罚成本是为了让模型在“多发电”与“经济性”之间取得合理折中惩罚系数通常设置为略高于上网电价否则模型会倾向于大量弃风。碳交易成本这里容易踩坑。很多人直接把碳排放量减去免费配额再乘以碳价但这个系统里存在两种碳流向CHP燃烧产生碳排放碳捕集设备从烟气中捕集一部分CO2后排放到大气中的才是净排放P2G设备如果制甲烷会消耗掉一部分捕集到的CO2。因此在目标函数里碳交易成本应当基于“净排放量”即净排放量 总碳排放量 - 碳捕集量 - 免费配额其中碳捕集量进入P2G制甲烷流程的那部分虽然物理上最终以甲烷燃烧的形式回到大气但在本模型的时间尺度内视为“已捕集暂存”。不同文献对这个问题的处理方式不太一样有的把P2G消耗的CO2作为二次排放处理我认为在日前调度这种短时间尺度下简单等效为“捕集即减排”是合理且符合直觉的。2.2 电、热、气三条平衡约束CHP是中间的枢纽节点热电联供系统的基本平衡关系并不复杂电能平衡、热能平衡、天然气平衡外加各设备的出力上下限和爬坡约束。但真正写起来CHP机组的建模方式是关键分水岭。我这里采用抽凝式CHP的简化可行域模型。所谓“可行域”就是在给定电出力Pe和热出力Hh的条件下机组能否稳定运行的区域。抽凝式机组的典型特征是可调热电比电出力提高时热出力可以在一个范围内变化而不是固定随动。这个可行域通常用几个线性不等式描述例如0 ≤ Hh ≤ H_max Pe_min(Hh) ≤ Pe ≤ Pe_max(Hh)其中电出力上下限是热出力的线性函数。写MATLAB约束时要把这些线性关系写成矩阵不等式放在同一个约束块里。如果用YALMIP直接对每个时段用循环append约束就行如果用纯intlinprog就要把所有时段的变量拼接成大向量然后用稀疏矩阵批量构建约束。前者代码可读性好后者求解效率高我后续会详细说这个取舍。电平衡约束购电 风电 光伏 CHP电出力 P2G消耗取负数 碳捕集电耗取负数 电负荷 电储能充电取负/放电取正。注意P2G和碳捕集在公式里不是“负荷”而是“可变电需求”它们的耗电量由调度决策内生决定这就让电平衡变得非线性耦合——但好消息是如果P2G效率是常数、碳捕集单位电耗是常数那么它们各自只是线性项问题仍然是线性规划。热平衡约束CHP热出力 热储能放热/充热 燃气锅炉供热如果有 热负荷。P2G在电转甲烷过程中产生的大量余热是否回收会直接影响热平衡。很多模型忽略这一项实际上P2G电解水过程的余热相当可观如果系统热负荷较大余热回收相当于免费热源。我建议建模时把这个回收系数加进去哪怕只是0.3左右也会明显改变调度结果。气平衡约束购气 P2G产气 储气罐放气 气负荷 CHP耗气 储气罐充气。这里有一个单位统一的问题——P2G的输入是MWh电输出是m³天然气而CHP耗气单位也是m³购气单位是m³但价格单位是元/m³。所有涉及能量转换的地方要么统一折算到MWh要么统一折算到m³不要混合使用。我当时在这个量纲问题上反复对了好几遍数据才意识到错在哪里。2.3 碳排放约束从简单系数到碳捕集参与平衡没有碳捕集时碳排放约束无非是“总排放 ≤ 配额 可购买额度”线性且简单。加入碳捕集后问题变得微妙捕集量不是独立决策变量它取决于CHP机组的烟气量和捕集设备的处理能力。如果要精确建模捕集量与CHP出力之间存在非线性关系——烟气量随燃料消耗量变化捕集能耗随CO2处理量变化。实际工程中更可行的路线是做分段线性化处理。我用的方案是以CHP电出力为横轴把碳捕集设备的捕获能力分三段描述每一段用不同的捕集率和不同单位能耗。在MATLAB中通过引入二进制变量和big-M约束来实现分段选择。这种做法避免了非线性求解器带来的收敛问题又能保留“低负荷时捕集率下降、单位能耗上升”的物理特性对于MILP求解器来说非常友好。3. 碳捕集与P2G的耦合细节从CO2流向看系统闭环3.1 碳捕集不是“即捕即走”它的能耗会反过来压制CHP出力碳捕集设备最反直觉的地方在于它捕集了碳却消耗了能量而这些能量恰好来自它服务的CHP机组。捕集系统的能耗分为两部分固定能耗开启溶剂循环泵、风机等和可变能耗再沸器再生溶剂的蒸汽消耗、压缩机的电耗。固定能耗的存在意味着碳捕集设备存在最小开启时间或最低处理量约束否则模型会为了减排目标频繁启停它这在物理上是不现实的。我建议把碳捕集设备建模为一个带有“最低运行负荷”的单元。假设设备设计捕集能力为C_maxt/h最低连续运行负荷为20%C_max低于这个值系统无法稳定运行。在日前调度里这意味着一旦碳捕集投入就必须保证捕集量不低于设计值的20%。加上这项约束后调度结果会呈现出更真实的模式碳捕集要么不开要么以一个可观的处理量运行而不是在低捕集率区间频繁抖动。另一个需要注意的问题是碳捕集消耗的热能来源。如果碳捕集再沸器需要用蒸汽而系统里有CHP机组那么捕集设备消耗的热能会直接挤占对外供热能力热平衡方程里要额外增加一项“碳捕集热耗”。如果忽略这一项相当于白送了碳捕集设备一份热能补贴结果会偏向于过度使用碳捕集。我做的模型中碳捕集能耗建模为电耗热耗双通道电耗进电平衡热耗进热平衡。3.2 P2G的电转甲烷路径CO2、氢气、天然气的三重耦合P2G技术路径分两步第一步电解水制氢第二步氢与CO2反应生成甲烷即甲烷化反应。从系统调度角度看第一步消耗电产出氢气和热量第二步消耗氢气与CO2产出甲烷和热量。如果只把P2G视为“电转气”中间商氢气不建模那么模型就退化为“一台电转甲烷机”这个简化对整体调度影响不大但会丢失一个重要约束——甲烷化反应需要连续的CO2供应而CO2供应量本身来自碳捕集设备二者之间存在严格的时序耦合关系。具体一点说设t时段碳捕集设备的捕集量为C(t)这部分CO2一部分直排被压缩存储另一部分输送给P2G。P2G在t时段能生产的甲烷量受到min(实际捕集到的CO2可用量氢气可用量)的约束。氢气可用量则来自电解槽产氢它有一个储氢缓冲罐所以氢气供应可以跨时段平滑。这种中间缓冲的存在使得P2G和碳捕集之间不是严格的同时刻耦合而是“当天累计量约束”写约束时要注意区分。我实际使用的是累计量限制一天内P2G消耗的总CO2不超过碳捕集设备当日送入存储的CO2总量加上初始存量。这比逐时约束更合理也更稳健因为没有必要把24小时物流都卡死在每个时段一一对应。逐时约束很容易导致模型无解或结果过于保守尤其是在碳捕集设备后半夜出力低、P2G却需要全天连续产气的场景。3.3 CO2储罐要不要建别为了好看加设备一个很容易被忽略的细节是CO2储罐。如果P2G和碳捕集在时序上不对齐大概率不对齐就必须有CO2缓存机制否则模型只能强行让两者同步运行经济性会大打折扣。我在第一版代码里没有CO2储罐结果求解器给出的碳捕集设备几乎全天满负荷只为了喂饱P2G看起来“绿色”实际上运行成本远超预期。加CO2储罐之后模型自由度增加了但也要付出代价储罐容量上限、初始CO2存量、期末存量约束都需要设置。我建议容量按“当日P2G总CO2需求量的30%~50%”来设置初始存量取容量的一半期末存量必须等于初始存量这样才能保证日滚动的连续性。这个约束在很多论文里不写但实际做连续多日调度时非常重要——如果期末CO2存量被模型悄悄清空第二天的调度就会陷入无米之炊的窘境。提示CO2储罐的建模与储气罐、储热罐几乎完全一样都是“容量上下限时段递推方程”。如果你已经有了储能设备的通用建模模板CO2储罐基本可以复用只是单位变成吨。4. MATLAB代码框架从数据结构到求解器调用的工程决策4.1 数据输入的模块化用结构体组织不要散落一堆变量写优化模型的MATLAB代码最忌讳的是把所有变量扔在工作区跑一个脚本全靠记忆。我的建议是把所有输入数据打包成一个结构体或类。例如caseData.time 1:24; caseData.load.e [...]; % 电负荷 caseData.load.h [...]; % 热负荷 caseData.load.g [...]; % 气负荷 caseData.renew.wind [...]; % 风电预测出力 caseData.renew.pv [...]; % 光伏预测出力 caseData.price.e [...]; % 分时电价 caseData.price.g [...]; % 气价 caseData.dev.chp struct(pmin,...,pmax,...,hmax,...); caseData.dev.p2g struct(eta_elec,0.75,eta_heat,0.2,co2_need,...); caseData.dev.ccs struct(capacity,...,power_cons,...,heat_cons,...); caseData.storage.co2 struct(cap,...,init,...);这样做的最大好处是换一个案例比如换典型日数据时只需要替换caseData这个结构体调度模型主体代码一个字都不用改。我后来给别人调试代码时发现很多问题都是因为数据散落在多段脚本里某一处改了另一处没改结果边际条件对不齐导致结果荒谬。4.2 YALMIP建模范式与纯intlinprog的选择可读性vs性能MATLAB里做MILP调度主流方案有两种。一种是YALMIP工具包配合Cplex/Gurobi等求解器用sdpvar定义变量用Constraints [Constraints; 表达式]不断追加约束代码跟数学公式几乎一一对应另一种是纯MATLAB的intlinprog把所有决策变量平铺成一个向量x用矩阵Aineq、bineq、Aeq、beq表达约束代码可读性差但不需要额外安装工具包。我实际推荐的做法分阶段建模初期用YALMIP快速尝试不同约束形式调试表达式、检查量纲平衡一旦模型结构和约束确定且求解顺利再考虑是否重写为纯intlinprog版本。如果你的场景规模不大比如24时段、几十台设备YALMIP版完全够用没必要为了性能牺牲可读性。如果场景扩大到8760小时或数百节点才需要考虑用矩阵形式重写并调用Cplex直接建模。YALMIP的关键代码骨架大致是这样PL sdpvar(1, 24); % 购电量 PG sdpvar(1, 24); % 购气量 Pe sdpvar(1, 24); % CHP电出力 Hh sdpvar(1, 24); % CHP热出力 Pp2g sdpvar(1, 24); % P2G电耗 Gp2g sdpvar(1, 24); % P2G产气量 Cccs sdpvar(1, 24); % 碳捕集量 Pccs sdpvar(1, 24); % 碳捕集电耗 % 储能变量... Constraints []; for t 1:24 Constraints [Constraints, Pe(t) PL(t) ... caseData.load.e(t) ...]; % 各类上下限约束 end Objective sum(...); ops sdpsettings(solver,cplex,verbose,2,showprogress,1); optimize(Constraints, Objective, ops);这一段看起来简单但真正常被忽略的是sdpsettings里的参数。我建议把verbose设为1或2这样求解器打印的日志能帮你在模型无解时快速定位是哪个约束导致不可行。还有mip tolerance即MIP相对间隙要设置到合理值默认值可能偏松导致结果在工程可接受误差边缘徘徊。我通常设为1e-4。4.3 求解之后的结果整理别只画几张图就完事求解完成后第一件事不是画图而是做“平衡校验”。把调度结果回代到三条平衡方程里逐时计算残差。残差超过1e-6的时段基本就是约束写错了或者变量搞混了。这种自查方式能节省大量调bug时间。我习惯用三个bar图呈现核心结果第一个图是电功率平衡堆叠图购电、风电、光伏、CHP电出力、P2G耗电、碳捕集耗电、电负荷第二个图是热平衡堆叠图CHP热出力、锅炉供热、热储放热、热负荷第三个图是CO2流向图排放、捕集、P2G消耗、储罐变化。三个图放在一张figure的不同subplot里一眼就能判断结果是否合理。我见过不少同行只输出目标值和几条关键曲线结果模型内部产生了明显的能量不平衡却无人察觉非常可惜。5. 调试中必须跨过的五道坎数值病态、非凸性与收敛问题5.1 不可行问题的定位方法一步一步“丢约束”找真凶MILP模型报不可行infeasible是最常见的挫败场景。我的排查方法很简单先把约束分组注释掉一部分看模型能否可行。比如先只保留电平衡和所有设备出力上下限跑一次如果能行再逐步加回热平衡、气平衡、储能约束加回来之后哪一步开始不可行问题就出在那部分。有一次我调试一个加入CO2储罐的版本模型一直不可行最后发现是期末CO2存量等于初始存量的约束与P2G全天消耗量、碳捕集最低运行负荷三者产生了矛盾——捕集设备出力虽然满足最低负荷但全天CO2产量填不满P2G需求和期末存量要求。解决方式是调整CO2储罐初始容量而不是放松约束。这类问题如果不对照“约束分块逐步回代”的流程很难快速锁定。另一个常见的不可行原因是孤立节点某时段储能设备SOC递推方程中的初始SOC和末尾SOC设置互相冲突。比如第1小时SOC初值0.5最后1小时要求SOC回到0.5但中间时段储能容量上限和充放电功率上限不足导致任何时候都满足不了“回到0.5”的要求。这时候要么降低期末要求要么增大储能额定容量。5.2 big-M参数不是越大越好M值失当会让求解器很痛苦MILP里用big-M处理逻辑约束非常普遍例如“P2G产气量只有在某种状态开启时才允许为正”这类约束。写起来很顺手但M取值过大会造成数值病态——Cplex或Gurobi在求解时会看到极大系数预求解器难以处理对偶间隙收敛困难最终要么超时要么结果偏差。我一般遵循一个原则M取该变量物理最大值的1.1~1.5倍。比如P2G电耗最大不超过5MW那M设为6~8即可没必要用1e6这种“安全值”。数值病态还会表现为同一次模型前后两次求解结果不一致或者在极小参数扰动下结果剧烈变化。如果怀疑M取值有问题可以用求解器自带的presolve日志看有没有“large coefficient”警告。5.3 单位不一致是最隐蔽的bug来源你必须统一能量单位综合能源系统里能量单位特别容易混电是MWh气是m³或kWh热是GJ碳是吨。每个变量单独看都没问题一旦放进同一个目标函数或同一个平衡方程就立刻出乱子。我的惯例是电、热、气全部折算成MWh气体按热值折算例如1 m³天然气热值约为10 kWh即0.01 MWh碳以吨为单位目标函数里CO2价格单位是元/t储能SOC用MWh表示不用百分比。在MATLAB里我建议在数据输入阶段就完成所有单位折算而不是在约束里反复乘系数。因为约束里系数越多越容易漏乘或错乘。折算完成后所有涉及单位转换的参数比如P2G产气热值只在数据定义处出现一次其他地方的公式看起来完全清爽。5.4 非线性环节的线性化分段逼近别太贪心碳捕集捕集率随负荷变化的曲线、电解槽效率随负载率变化的曲线本质都是非线性。有些同行为了“精确”建模把非线性保留了然后调用非线性求解器结果经常掉进局部最优或干脆无法收敛。我认为综合能源系统日前调度的精度不足以支撑“精确非线性”的投入产出比分段线性化是更务实的选择。分段线性化的关键是把非线性曲线分成3~5段每一段用线性关系近似。例如电解槽效率从0.65到0.85、随负载率变化我把它分成4段负载率0~25%、25%~50%、50%~75%、75%~100%每段用不同的线性效率常数。然后用sos2约束或引入二进制变量表示当前处于哪一段。YALMIP有直接的分段函数处理工具但手写约束也不难。注意分段过多会增加二进制变量拖慢求解。4段对大多数工程场景已经足够不必为了几个百分点的精度让求解时间翻好几倍。5.5 求解器超时与次优解这是一个工程决策而非纯粹数学问题即便用了线性化包含二进制变量的MILP在较大场景下也可能求解缓慢。我遇到过一个含50个二进制变量的24时段模型默认参数下Cplex要跑15分钟才收敛到可行解。这个速度对研究型项目还能忍但如果你想做场景对比分析比如跑20个典型日就会非常痛苦。应对方式有几种一是调大MIP相对间隙例如从默认1e-4放宽到5e-4或1e-3在工程上这通常只会带来1%以内的目标值偏差但求解时间可能缩短5~10倍二是给二进制变量设置合理的start point通过warm start引导求解器进入好的分支搜索路径三是减少不必要的整数变量数量比如把一些原本用二进制表示的“最小运行负荷”约束改为不启用先看对结果的影响再决定是否加回来。6. 让模型能落地的几个扩展方向6.1 从日前调度到滚动调度模型结构不需要大改日前调度模型往往是一天求解一次但实际系统中风电、光伏、负荷预测值会随时间修正因此更实用的做法是“滚动时域优化”即每1小时求解一次未来24小时的调度方案只执行当前时段的指令。好消息是如果你用的是YALMIP结构体数据组织方式从日前到滚动式改动非常小只需要在每次求解前更新caseData中的预测曲线并把储能/储气/CO2储罐的实际SOC或存量作为初始值传入。这里唯一要小心的是终止条件——最后一次滚动时段的最终SOC或存量要求要重新设定不能让模型“为了回初值而回初值”。6.2 考虑不确定性鲁棒优化与场景法的取舍如果觉得确定性模型太理想化下一个自然的方向就是风电/光伏出力的不确定性建模。我自己的经验是场景法scenario-based stochastic programming比分布鲁棒法更容易在MATLAB里落地只需要把风电场景集提前用历史数据生成比如k-means聚类出10个典型场景然后在目标函数里除以场景数量取平均。但场景法会迅速扩大变量和约束规模。处理好这个问题的核心在于场景缩减——把100个原始场景缩减为10个代表性场景。MATLAB自带的kmeans函数就能完成大部分工作。缩减后的场景如果仍然导致求解过慢还可以用多阶段决策分解Benders分解是最常见的不过那已经超出了绝大多数研究项目的需求。6.3 从单一园区到多园区互联新增的是联络线约束如果模型要扩展到多个综合能源系统区域互联物理结构只是增加了区域间联络线通常用一条功率交换线来建模。核心约束变为联络线传输功率不超过上限区域内各设备运行逻辑不变目标函数再增加跨区购售电成本项。拓展时最容易出错的其实是联络线功率方向约束的方向一致性A区输出到B区的功率在A区的电平衡方程里是“出力”在B区的电平衡方程里是“负荷”符号取反后必须确保两者同时被约束在同一个界面储能模型上。在MATLAB里把所有区域变量放在同一个求解向量中联络线功率用同一个变量表示并加正负号约束能最大程度避免这个坑。6.4 加入需求响应用户侧灵活性才是隐藏红利模型做到后面你会发现供给侧设备CHP、P2G、碳捕集、储能都调完了还是无法完全消纳全部风电特别是某些极端气象日。这时候需求响应就非常有价值一部分可转移负荷比如工厂电解铝、充电桩聚合可以平移1~2小时可削减负荷可以接受削峰指令。在调度模型里加入需求响应本质上是增加电负荷平衡方程的调整变量。可转移负荷要建模为“总量不变、时段可变”引入二进制变量和时段关联约束可削减负荷则相对简单只需在规定上下限范围内允许削减并支付补偿费用。这些扩展对原模型的影响是局部性的不会推翻CHP、P2G、碳捕集的耦合约束因此非常值得加。写在最后的一些实际体会这套模型我前后迭代了三个版本从第一版不考虑CO2储罐、直接用逐时耦合约束到第三版引入CO2储罐累计约束和分段线性化最大的体会是综合能源系统优化调度的价值恰恰藏在那些“看似多余”的细节里。碳捕集的最低运行负荷、P2G的余热回收、CO2储罐的初始状态——每一条都在约束模型但每一条背后都是真实物理系统对调度计划的无情限制。模型不是越复杂越好而是越贴近设备真实运行约束越好。如果你正在写类似的MATLAB代码我的建议是先把“无解调试”这一关想清楚。准备一个约束分块开关系统比如一个struct里放一些布尔量每个布尔量控制一组约束是否启用。这样调试时你就不需要不断注释代码而是通过修改这个struct快速定位问题。这个小小的习惯能替你省下大量“盲人摸象”式的排查时间。最后一个实操小技巧求解得到结果后手动检查几个典型时段的数据。比如深夜低负荷时段电储能不能同时充电、P2G同时耗电、碳捕集同时捕碳——这三者都在抢电系统会不会出现电不够用的情况如果调度结果告诉你这一刻购电量比正常午高峰还高你就该回头看看是不是约束方向写反了。这种“对业务直觉的复盘”比任何调试工具都靠谱。
1. 为什么P2G和碳捕集要放进同一个模型问题的起点与耦合逻辑做综合能源系统优化调度的人早晚都会遇到这么一个问题电热联供机组CHP一开热是满足了但电出力跟着被绑死风电光伏大发的时候电又多到没处去只能弃风弃光碳排放又随着煤/气消耗一路上涨。单独加一套电转气P2G设备能把多余的电变成天然气存起来但P2G制甲烷恰恰需要二氧化碳作为原料。单独加一套碳捕集设备能降低火电/CHP的碳排放但捕集本身要消耗大量电能和热能又反过来挤压系统的供电能力。两个设备单看都有价值放在一起却产生了天然的双向互补关系——碳捕集提供CO2P2G消耗CO2P2G消耗电能碳捕集也消耗电能但同时降低了系统碳排放。我最初做这个课题时想法很朴素在传统的热电联供综合能源系统里把电储能、热储能、气储能都加上再引入P2G设备和碳捕集设备然后看它们之间如何协同调度。但真正动手建模、写MATLAB代码、调结果之后才发现这个模型的核心难点根本不在“加设备”而在于三个层面一是两种新设备与CHP机组之间形成了复杂的能量-物质耦合闭环二是碳捕集设备的能耗特性会显著改变系统电热平衡关系三是这些设备的运行约束直接影响了日前调度的可行域处理不好就会出现一堆匪夷所思的劣质结果。这篇博文就把我实际搭建这个调度模型的过程、踩过的坑、MATLAB代码的组织方式以及一些值得注意的建模细节原原本本整理一遍。适合正在做综合能源系统优化调度、碳捕集与P2G耦合建模的同行参考也适合刚接触这类模型、想用MATLAB快速搭一个可运行版本的朋友做地图式导览。我不会把代码逐行粘贴因为那会变成一份没有灵魂的代码说明书我会重点讲清楚“每一部分代码背后的建模决策是什么为什么这么写换了场景要改哪里”。2. 模型骨架目标函数与约束条件的数学化拆解2.1 目标函数的费用构成别只盯着购能成本典型的日前优化调度目标函数一般是总运行成本最小化。我实际使用的目标函数包含四块购能成本购电、购气、碳交易成本、设备运维成本、弃风弃光惩罚成本。之所以单独列出碳交易成本是因为碳捕集设备的收益必须在经济性上体现出来如果不能把“捕了多少碳”换算成钱模型就会本能地“不想开”碳捕集设备。购电成本按分时电价计算这个比较简单——逐时电价向量点乘逐时购电量。购气成本类似按天然气价格计算注意CHP耗气和P2G产气都要在气平衡中体现。运维成本我习惯按出力比例计算单位是“元/MWh”而不是固定启停费用因为对于日前连续调度来说边际运维成本对调度结果的影响更平滑。弃风弃光惩罚成本是为了让模型在“多发电”与“经济性”之间取得合理折中惩罚系数通常设置为略高于上网电价否则模型会倾向于大量弃风。碳交易成本这里容易踩坑。很多人直接把碳排放量减去免费配额再乘以碳价但这个系统里存在两种碳流向CHP燃烧产生碳排放碳捕集设备从烟气中捕集一部分CO2后排放到大气中的才是净排放P2G设备如果制甲烷会消耗掉一部分捕集到的CO2。因此在目标函数里碳交易成本应当基于“净排放量”即净排放量 总碳排放量 - 碳捕集量 - 免费配额其中碳捕集量进入P2G制甲烷流程的那部分虽然物理上最终以甲烷燃烧的形式回到大气但在本模型的时间尺度内视为“已捕集暂存”。不同文献对这个问题的处理方式不太一样有的把P2G消耗的CO2作为二次排放处理我认为在日前调度这种短时间尺度下简单等效为“捕集即减排”是合理且符合直觉的。2.2 电、热、气三条平衡约束CHP是中间的枢纽节点热电联供系统的基本平衡关系并不复杂电能平衡、热能平衡、天然气平衡外加各设备的出力上下限和爬坡约束。但真正写起来CHP机组的建模方式是关键分水岭。我这里采用抽凝式CHP的简化可行域模型。所谓“可行域”就是在给定电出力Pe和热出力Hh的条件下机组能否稳定运行的区域。抽凝式机组的典型特征是可调热电比电出力提高时热出力可以在一个范围内变化而不是固定随动。这个可行域通常用几个线性不等式描述例如0 ≤ Hh ≤ H_max Pe_min(Hh) ≤ Pe ≤ Pe_max(Hh)其中电出力上下限是热出力的线性函数。写MATLAB约束时要把这些线性关系写成矩阵不等式放在同一个约束块里。如果用YALMIP直接对每个时段用循环append约束就行如果用纯intlinprog就要把所有时段的变量拼接成大向量然后用稀疏矩阵批量构建约束。前者代码可读性好后者求解效率高我后续会详细说这个取舍。电平衡约束购电 风电 光伏 CHP电出力 P2G消耗取负数 碳捕集电耗取负数 电负荷 电储能充电取负/放电取正。注意P2G和碳捕集在公式里不是“负荷”而是“可变电需求”它们的耗电量由调度决策内生决定这就让电平衡变得非线性耦合——但好消息是如果P2G效率是常数、碳捕集单位电耗是常数那么它们各自只是线性项问题仍然是线性规划。热平衡约束CHP热出力 热储能放热/充热 燃气锅炉供热如果有 热负荷。P2G在电转甲烷过程中产生的大量余热是否回收会直接影响热平衡。很多模型忽略这一项实际上P2G电解水过程的余热相当可观如果系统热负荷较大余热回收相当于免费热源。我建议建模时把这个回收系数加进去哪怕只是0.3左右也会明显改变调度结果。气平衡约束购气 P2G产气 储气罐放气 气负荷 CHP耗气 储气罐充气。这里有一个单位统一的问题——P2G的输入是MWh电输出是m³天然气而CHP耗气单位也是m³购气单位是m³但价格单位是元/m³。所有涉及能量转换的地方要么统一折算到MWh要么统一折算到m³不要混合使用。我当时在这个量纲问题上反复对了好几遍数据才意识到错在哪里。2.3 碳排放约束从简单系数到碳捕集参与平衡没有碳捕集时碳排放约束无非是“总排放 ≤ 配额 可购买额度”线性且简单。加入碳捕集后问题变得微妙捕集量不是独立决策变量它取决于CHP机组的烟气量和捕集设备的处理能力。如果要精确建模捕集量与CHP出力之间存在非线性关系——烟气量随燃料消耗量变化捕集能耗随CO2处理量变化。实际工程中更可行的路线是做分段线性化处理。我用的方案是以CHP电出力为横轴把碳捕集设备的捕获能力分三段描述每一段用不同的捕集率和不同单位能耗。在MATLAB中通过引入二进制变量和big-M约束来实现分段选择。这种做法避免了非线性求解器带来的收敛问题又能保留“低负荷时捕集率下降、单位能耗上升”的物理特性对于MILP求解器来说非常友好。3. 碳捕集与P2G的耦合细节从CO2流向看系统闭环3.1 碳捕集不是“即捕即走”它的能耗会反过来压制CHP出力碳捕集设备最反直觉的地方在于它捕集了碳却消耗了能量而这些能量恰好来自它服务的CHP机组。捕集系统的能耗分为两部分固定能耗开启溶剂循环泵、风机等和可变能耗再沸器再生溶剂的蒸汽消耗、压缩机的电耗。固定能耗的存在意味着碳捕集设备存在最小开启时间或最低处理量约束否则模型会为了减排目标频繁启停它这在物理上是不现实的。我建议把碳捕集设备建模为一个带有“最低运行负荷”的单元。假设设备设计捕集能力为C_maxt/h最低连续运行负荷为20%C_max低于这个值系统无法稳定运行。在日前调度里这意味着一旦碳捕集投入就必须保证捕集量不低于设计值的20%。加上这项约束后调度结果会呈现出更真实的模式碳捕集要么不开要么以一个可观的处理量运行而不是在低捕集率区间频繁抖动。另一个需要注意的问题是碳捕集消耗的热能来源。如果碳捕集再沸器需要用蒸汽而系统里有CHP机组那么捕集设备消耗的热能会直接挤占对外供热能力热平衡方程里要额外增加一项“碳捕集热耗”。如果忽略这一项相当于白送了碳捕集设备一份热能补贴结果会偏向于过度使用碳捕集。我做的模型中碳捕集能耗建模为电耗热耗双通道电耗进电平衡热耗进热平衡。3.2 P2G的电转甲烷路径CO2、氢气、天然气的三重耦合P2G技术路径分两步第一步电解水制氢第二步氢与CO2反应生成甲烷即甲烷化反应。从系统调度角度看第一步消耗电产出氢气和热量第二步消耗氢气与CO2产出甲烷和热量。如果只把P2G视为“电转气”中间商氢气不建模那么模型就退化为“一台电转甲烷机”这个简化对整体调度影响不大但会丢失一个重要约束——甲烷化反应需要连续的CO2供应而CO2供应量本身来自碳捕集设备二者之间存在严格的时序耦合关系。具体一点说设t时段碳捕集设备的捕集量为C(t)这部分CO2一部分直排被压缩存储另一部分输送给P2G。P2G在t时段能生产的甲烷量受到min(实际捕集到的CO2可用量氢气可用量)的约束。氢气可用量则来自电解槽产氢它有一个储氢缓冲罐所以氢气供应可以跨时段平滑。这种中间缓冲的存在使得P2G和碳捕集之间不是严格的同时刻耦合而是“当天累计量约束”写约束时要注意区分。我实际使用的是累计量限制一天内P2G消耗的总CO2不超过碳捕集设备当日送入存储的CO2总量加上初始存量。这比逐时约束更合理也更稳健因为没有必要把24小时物流都卡死在每个时段一一对应。逐时约束很容易导致模型无解或结果过于保守尤其是在碳捕集设备后半夜出力低、P2G却需要全天连续产气的场景。3.3 CO2储罐要不要建别为了好看加设备一个很容易被忽略的细节是CO2储罐。如果P2G和碳捕集在时序上不对齐大概率不对齐就必须有CO2缓存机制否则模型只能强行让两者同步运行经济性会大打折扣。我在第一版代码里没有CO2储罐结果求解器给出的碳捕集设备几乎全天满负荷只为了喂饱P2G看起来“绿色”实际上运行成本远超预期。加CO2储罐之后模型自由度增加了但也要付出代价储罐容量上限、初始CO2存量、期末存量约束都需要设置。我建议容量按“当日P2G总CO2需求量的30%~50%”来设置初始存量取容量的一半期末存量必须等于初始存量这样才能保证日滚动的连续性。这个约束在很多论文里不写但实际做连续多日调度时非常重要——如果期末CO2存量被模型悄悄清空第二天的调度就会陷入无米之炊的窘境。提示CO2储罐的建模与储气罐、储热罐几乎完全一样都是“容量上下限时段递推方程”。如果你已经有了储能设备的通用建模模板CO2储罐基本可以复用只是单位变成吨。4. MATLAB代码框架从数据结构到求解器调用的工程决策4.1 数据输入的模块化用结构体组织不要散落一堆变量写优化模型的MATLAB代码最忌讳的是把所有变量扔在工作区跑一个脚本全靠记忆。我的建议是把所有输入数据打包成一个结构体或类。例如caseData.time 1:24; caseData.load.e [...]; % 电负荷 caseData.load.h [...]; % 热负荷 caseData.load.g [...]; % 气负荷 caseData.renew.wind [...]; % 风电预测出力 caseData.renew.pv [...]; % 光伏预测出力 caseData.price.e [...]; % 分时电价 caseData.price.g [...]; % 气价 caseData.dev.chp struct(pmin,...,pmax,...,hmax,...); caseData.dev.p2g struct(eta_elec,0.75,eta_heat,0.2,co2_need,...); caseData.dev.ccs struct(capacity,...,power_cons,...,heat_cons,...); caseData.storage.co2 struct(cap,...,init,...);这样做的最大好处是换一个案例比如换典型日数据时只需要替换caseData这个结构体调度模型主体代码一个字都不用改。我后来给别人调试代码时发现很多问题都是因为数据散落在多段脚本里某一处改了另一处没改结果边际条件对不齐导致结果荒谬。4.2 YALMIP建模范式与纯intlinprog的选择可读性vs性能MATLAB里做MILP调度主流方案有两种。一种是YALMIP工具包配合Cplex/Gurobi等求解器用sdpvar定义变量用Constraints [Constraints; 表达式]不断追加约束代码跟数学公式几乎一一对应另一种是纯MATLAB的intlinprog把所有决策变量平铺成一个向量x用矩阵Aineq、bineq、Aeq、beq表达约束代码可读性差但不需要额外安装工具包。我实际推荐的做法分阶段建模初期用YALMIP快速尝试不同约束形式调试表达式、检查量纲平衡一旦模型结构和约束确定且求解顺利再考虑是否重写为纯intlinprog版本。如果你的场景规模不大比如24时段、几十台设备YALMIP版完全够用没必要为了性能牺牲可读性。如果场景扩大到8760小时或数百节点才需要考虑用矩阵形式重写并调用Cplex直接建模。YALMIP的关键代码骨架大致是这样PL sdpvar(1, 24); % 购电量 PG sdpvar(1, 24); % 购气量 Pe sdpvar(1, 24); % CHP电出力 Hh sdpvar(1, 24); % CHP热出力 Pp2g sdpvar(1, 24); % P2G电耗 Gp2g sdpvar(1, 24); % P2G产气量 Cccs sdpvar(1, 24); % 碳捕集量 Pccs sdpvar(1, 24); % 碳捕集电耗 % 储能变量... Constraints []; for t 1:24 Constraints [Constraints, Pe(t) PL(t) ... caseData.load.e(t) ...]; % 各类上下限约束 end Objective sum(...); ops sdpsettings(solver,cplex,verbose,2,showprogress,1); optimize(Constraints, Objective, ops);这一段看起来简单但真正常被忽略的是sdpsettings里的参数。我建议把verbose设为1或2这样求解器打印的日志能帮你在模型无解时快速定位是哪个约束导致不可行。还有mip tolerance即MIP相对间隙要设置到合理值默认值可能偏松导致结果在工程可接受误差边缘徘徊。我通常设为1e-4。4.3 求解之后的结果整理别只画几张图就完事求解完成后第一件事不是画图而是做“平衡校验”。把调度结果回代到三条平衡方程里逐时计算残差。残差超过1e-6的时段基本就是约束写错了或者变量搞混了。这种自查方式能节省大量调bug时间。我习惯用三个bar图呈现核心结果第一个图是电功率平衡堆叠图购电、风电、光伏、CHP电出力、P2G耗电、碳捕集耗电、电负荷第二个图是热平衡堆叠图CHP热出力、锅炉供热、热储放热、热负荷第三个图是CO2流向图排放、捕集、P2G消耗、储罐变化。三个图放在一张figure的不同subplot里一眼就能判断结果是否合理。我见过不少同行只输出目标值和几条关键曲线结果模型内部产生了明显的能量不平衡却无人察觉非常可惜。5. 调试中必须跨过的五道坎数值病态、非凸性与收敛问题5.1 不可行问题的定位方法一步一步“丢约束”找真凶MILP模型报不可行infeasible是最常见的挫败场景。我的排查方法很简单先把约束分组注释掉一部分看模型能否可行。比如先只保留电平衡和所有设备出力上下限跑一次如果能行再逐步加回热平衡、气平衡、储能约束加回来之后哪一步开始不可行问题就出在那部分。有一次我调试一个加入CO2储罐的版本模型一直不可行最后发现是期末CO2存量等于初始存量的约束与P2G全天消耗量、碳捕集最低运行负荷三者产生了矛盾——捕集设备出力虽然满足最低负荷但全天CO2产量填不满P2G需求和期末存量要求。解决方式是调整CO2储罐初始容量而不是放松约束。这类问题如果不对照“约束分块逐步回代”的流程很难快速锁定。另一个常见的不可行原因是孤立节点某时段储能设备SOC递推方程中的初始SOC和末尾SOC设置互相冲突。比如第1小时SOC初值0.5最后1小时要求SOC回到0.5但中间时段储能容量上限和充放电功率上限不足导致任何时候都满足不了“回到0.5”的要求。这时候要么降低期末要求要么增大储能额定容量。5.2 big-M参数不是越大越好M值失当会让求解器很痛苦MILP里用big-M处理逻辑约束非常普遍例如“P2G产气量只有在某种状态开启时才允许为正”这类约束。写起来很顺手但M取值过大会造成数值病态——Cplex或Gurobi在求解时会看到极大系数预求解器难以处理对偶间隙收敛困难最终要么超时要么结果偏差。我一般遵循一个原则M取该变量物理最大值的1.1~1.5倍。比如P2G电耗最大不超过5MW那M设为6~8即可没必要用1e6这种“安全值”。数值病态还会表现为同一次模型前后两次求解结果不一致或者在极小参数扰动下结果剧烈变化。如果怀疑M取值有问题可以用求解器自带的presolve日志看有没有“large coefficient”警告。5.3 单位不一致是最隐蔽的bug来源你必须统一能量单位综合能源系统里能量单位特别容易混电是MWh气是m³或kWh热是GJ碳是吨。每个变量单独看都没问题一旦放进同一个目标函数或同一个平衡方程就立刻出乱子。我的惯例是电、热、气全部折算成MWh气体按热值折算例如1 m³天然气热值约为10 kWh即0.01 MWh碳以吨为单位目标函数里CO2价格单位是元/t储能SOC用MWh表示不用百分比。在MATLAB里我建议在数据输入阶段就完成所有单位折算而不是在约束里反复乘系数。因为约束里系数越多越容易漏乘或错乘。折算完成后所有涉及单位转换的参数比如P2G产气热值只在数据定义处出现一次其他地方的公式看起来完全清爽。5.4 非线性环节的线性化分段逼近别太贪心碳捕集捕集率随负荷变化的曲线、电解槽效率随负载率变化的曲线本质都是非线性。有些同行为了“精确”建模把非线性保留了然后调用非线性求解器结果经常掉进局部最优或干脆无法收敛。我认为综合能源系统日前调度的精度不足以支撑“精确非线性”的投入产出比分段线性化是更务实的选择。分段线性化的关键是把非线性曲线分成3~5段每一段用线性关系近似。例如电解槽效率从0.65到0.85、随负载率变化我把它分成4段负载率0~25%、25%~50%、50%~75%、75%~100%每段用不同的线性效率常数。然后用sos2约束或引入二进制变量表示当前处于哪一段。YALMIP有直接的分段函数处理工具但手写约束也不难。注意分段过多会增加二进制变量拖慢求解。4段对大多数工程场景已经足够不必为了几个百分点的精度让求解时间翻好几倍。5.5 求解器超时与次优解这是一个工程决策而非纯粹数学问题即便用了线性化包含二进制变量的MILP在较大场景下也可能求解缓慢。我遇到过一个含50个二进制变量的24时段模型默认参数下Cplex要跑15分钟才收敛到可行解。这个速度对研究型项目还能忍但如果你想做场景对比分析比如跑20个典型日就会非常痛苦。应对方式有几种一是调大MIP相对间隙例如从默认1e-4放宽到5e-4或1e-3在工程上这通常只会带来1%以内的目标值偏差但求解时间可能缩短5~10倍二是给二进制变量设置合理的start point通过warm start引导求解器进入好的分支搜索路径三是减少不必要的整数变量数量比如把一些原本用二进制表示的“最小运行负荷”约束改为不启用先看对结果的影响再决定是否加回来。6. 让模型能落地的几个扩展方向6.1 从日前调度到滚动调度模型结构不需要大改日前调度模型往往是一天求解一次但实际系统中风电、光伏、负荷预测值会随时间修正因此更实用的做法是“滚动时域优化”即每1小时求解一次未来24小时的调度方案只执行当前时段的指令。好消息是如果你用的是YALMIP结构体数据组织方式从日前到滚动式改动非常小只需要在每次求解前更新caseData中的预测曲线并把储能/储气/CO2储罐的实际SOC或存量作为初始值传入。这里唯一要小心的是终止条件——最后一次滚动时段的最终SOC或存量要求要重新设定不能让模型“为了回初值而回初值”。6.2 考虑不确定性鲁棒优化与场景法的取舍如果觉得确定性模型太理想化下一个自然的方向就是风电/光伏出力的不确定性建模。我自己的经验是场景法scenario-based stochastic programming比分布鲁棒法更容易在MATLAB里落地只需要把风电场景集提前用历史数据生成比如k-means聚类出10个典型场景然后在目标函数里除以场景数量取平均。但场景法会迅速扩大变量和约束规模。处理好这个问题的核心在于场景缩减——把100个原始场景缩减为10个代表性场景。MATLAB自带的kmeans函数就能完成大部分工作。缩减后的场景如果仍然导致求解过慢还可以用多阶段决策分解Benders分解是最常见的不过那已经超出了绝大多数研究项目的需求。6.3 从单一园区到多园区互联新增的是联络线约束如果模型要扩展到多个综合能源系统区域互联物理结构只是增加了区域间联络线通常用一条功率交换线来建模。核心约束变为联络线传输功率不超过上限区域内各设备运行逻辑不变目标函数再增加跨区购售电成本项。拓展时最容易出错的其实是联络线功率方向约束的方向一致性A区输出到B区的功率在A区的电平衡方程里是“出力”在B区的电平衡方程里是“负荷”符号取反后必须确保两者同时被约束在同一个界面储能模型上。在MATLAB里把所有区域变量放在同一个求解向量中联络线功率用同一个变量表示并加正负号约束能最大程度避免这个坑。6.4 加入需求响应用户侧灵活性才是隐藏红利模型做到后面你会发现供给侧设备CHP、P2G、碳捕集、储能都调完了还是无法完全消纳全部风电特别是某些极端气象日。这时候需求响应就非常有价值一部分可转移负荷比如工厂电解铝、充电桩聚合可以平移1~2小时可削减负荷可以接受削峰指令。在调度模型里加入需求响应本质上是增加电负荷平衡方程的调整变量。可转移负荷要建模为“总量不变、时段可变”引入二进制变量和时段关联约束可削减负荷则相对简单只需在规定上下限范围内允许削减并支付补偿费用。这些扩展对原模型的影响是局部性的不会推翻CHP、P2G、碳捕集的耦合约束因此非常值得加。写在最后的一些实际体会这套模型我前后迭代了三个版本从第一版不考虑CO2储罐、直接用逐时耦合约束到第三版引入CO2储罐累计约束和分段线性化最大的体会是综合能源系统优化调度的价值恰恰藏在那些“看似多余”的细节里。碳捕集的最低运行负荷、P2G的余热回收、CO2储罐的初始状态——每一条都在约束模型但每一条背后都是真实物理系统对调度计划的无情限制。模型不是越复杂越好而是越贴近设备真实运行约束越好。如果你正在写类似的MATLAB代码我的建议是先把“无解调试”这一关想清楚。准备一个约束分块开关系统比如一个struct里放一些布尔量每个布尔量控制一组约束是否启用。这样调试时你就不需要不断注释代码而是通过修改这个struct快速定位问题。这个小小的习惯能替你省下大量“盲人摸象”式的排查时间。最后一个实操小技巧求解得到结果后手动检查几个典型时段的数据。比如深夜低负荷时段电储能不能同时充电、P2G同时耗电、碳捕集同时捕碳——这三者都在抢电系统会不会出现电不够用的情况如果调度结果告诉你这一刻购电量比正常午高峰还高你就该回头看看是不是约束方向写反了。这种“对业务直觉的复盘”比任何调试工具都靠谱。
