做综合能源系统优化电、热、气、碳四条流叠在一起设备一多Matlab里的优化变量很容易就到两三百个。这个课题我前后改了四版模型最卡壳的不是数学推导而是怎么把热电联产、电转气、碳捕集三个子系统放进同一个优化框架里还能自洽。这篇文章就按实际做项目的顺序把建模思路、设备数学化、Matlab代码写法、算例效果和排错经验完整过一遍适合正在做综合能源系统调度优化、或者准备用Matlab加Yalmip做多能互补模型的同学直接参考。先给结论P2G和CCS一旦接入系统就不再是简单的电热联供而是一个带物质循环的“能量-碳闭环”。优化时除了盯电功率和热功率平衡还要盯天然气流量和CO2流量这是电气背景的同行第一次建模时最不适应的点。文中的参数按常见论文取值的量级设定代码骨架可以直接改成你自己的设备容量和负荷数据。1. 综合能源系统的整体思路为什么要把P2G和CCS绑在一起1.1 传统热电联产的问题在哪里热电联产机组CHP最大的特点是“以热定电”当冬季热负荷高时机组必须维持较高的电出力来满足供热否则热就不够用。问题随之而来——如果此时新能源大发风电和光伏的出力加上CHP的强制电出力会远超电负荷电网又没那么大消纳能力只能弃风弃光。我见过不少项目里弃风率最高的时段不是半夜就是冬季低温天核心矛盾就是CHP的热电耦合把电出力下限抬得太高。电转气就是用来消化这部分富余电能的。多出来的电不直接并网而是送去电解水制氢再和CO2反应合成天然气氢气不好储存、不好运输但天然气可以进管网、进储气罐等用电高峰或者气价高的时候再放出来。相当于给综合能源系统装了一个大号“能量缓存的阀门”。碳捕集的角色同样清晰。CHP和燃气锅炉都要烧天然气烟气里的CO2是主要碳排放来源。CCS把这些CO2捕下来之后有两种去处一种是压缩封存一种是送给P2G做甲烷化的原料。前一种纯花钱后一种能把碳变成产品还能减少外购的CO2成本。所以把P2G和CCS放在同一个系统里不是两个技术简单叠加而是让CO2流、天然气流、电力和热力四条物质能量流形成一个闭环。1.2 系统拓扑与能量流梳理在搭建模型之前我建议先画一张拓扑图哪怕只是手绘的方框图也行。我的系统里主要包含这些单元外部电网可购电设分时电价购电上限根据变电站容量定。风力发电、光伏发电预测出力曲线作为输入超出消纳部分直接弃掉弃风弃光有惩罚成本。燃气轮机CHP机组输出电和热消耗天然气烟气送CCS。燃气锅炉作为热负荷补充和CHP共同满足供热需求。电锅炉和储热罐电锅炉在低谷电时段把电转成热储起来高峰时段放热。P2G单元消耗电能和CO2产出天然气可进入储气罐或直接供气。碳捕集系统捕集CHP烟气中的CO2捕集过程本身消耗电产出的CO2供给P2G或封存。储气罐缓冲天然气供需波动。调度周期取24小时时间粒度1小时。这个粒度是论文里最常用的因为电价、气价、负荷预测、新能源预测都能方便地对齐到小时级。想细化到15分钟也行但模型维度和求解时间会明显上升第一版不建议。从能量流角度看电平衡约束要接进CHP、风电光伏、P2G、CCS、储能和购电热平衡约束要接进CHP、燃气锅炉、电锅炉和储热气平衡约束要接进天然气管网购气、P2G产气、燃气锅炉耗气和CHP耗气CO2平衡约束则接进CHP排放、CCS捕集、P2G消耗和封存。四条平衡线一拉出来整个优化模型的大纲就出来了。1.3 优化目标和求解框架怎么选优化目标我选的是日运行成本最小包括购电费用、购气费用、碳交易费用、弃风弃光惩罚和设备启停费用。碳排放不单独作为硬约束而是用碳配额交易的方式放进目标函数给定每日免费配额实际排放超过配额就需要按碳价购买配额低于配额则可以把富余配额卖出碳价固定为80元/吨。这样碳排放和经济效益直接挂钩论文里也更好解释结果。为什么不用非线性模型因为P2G的电解效率、CHP的可行域、CCS能耗这些关系在线性化之后完全能满足工程精度而且MILP能保证收敛到全局最优。非线性模型就算能跑通也很难判断你拿到的解是不是局部最优。我第一版用的就是MILP配合Yalmip建模调用Gurobi求解稳定、快速、后处理也方便。确定性MILP先跑对之后再扩展鲁棒优化或者随机优化都不迟。2. 核心设备模型怎么数学化2.1 热电联产机组别再用固定热电比糊弄了很多入门代码喜欢把CHP写成固定热电比比如H1.2P这样省事但问题很大。真实抽凝式机组的电热可行域是一个平面多边形运行点只能落在这个多边形里面固定热电比相当于把可行域压成一条线会丢失大量可行的调度空间。做优化课题这一步最好一次到位。最通用的做法是顶点凸组合法。把CHP在“电功率-热功率”平面上的可行域近似成一个凸四边形取四个顶点设每个顶点的权重为λ1到λ4则任意运行点可以表达为四个顶点坐标按权重叠加权重之和为1且非负。这种写法在Yalmip里非常直观约束数量少求解器也喜欢。我常用的近似参数取这样一组顶点坐标单位MW顶点号电出力P热出力H140102100031009044060对应约束就是40≤P≤1000≤H≤90再加上两条上下包络线的不等式。实际项目中机组厂家会提供可行域图你直接把顶点坐标换成真实数据即可。凸组合的好处是自动保证P和H之间的耦合关系不需要额外写一堆分段线性近似。爬坡约束也要写进去。我一般用一阶差分相邻两个时段电出力之差不超过机组爬坡速率比如±20MW/h。注意热出力一般跟随电出力变化如果不想引入太复杂的联合爬坡可以先只对电出力加爬坡限制热出力通过可行域约束间接限制工程上够用。2.2 电转气效率、热值和CO2消耗量的换算电转气系统内部有两段反应电解水制氢再到甲烷化反应器里让氢气和CO2合成甲烷。从优化模型的角度不需要关心中间氢气的储能细节直接做能量层面的单输入单输出处理就行。P2G的输入是电功率输出是天然气功率中间用综合效率η_p2g联系起来取值一般0.55到0.65我代码里取0.6。这里最容易翻车的是单位换算。电功率单位是MW天然气如果按体积算就是m³/h两者必须通过热值桥接。天然气低热值LNG一般取9.7kWh/m³左右也就是说1MWh的天然气能量约等于103m³的天然气。如果P2G满发50MW输入电能50MWh输出燃气能量就是0.6×5030MWh折合气量约3093m³/h。这会直接决定后面CO2消耗量的计算。按甲烷化反应的化学计量关系每生产1m³的CH4大概需要1.96kg的CO2。换算到能量单位就是每MWh的天然气产品约消耗0.2吨CO2。我代码里写成alpha_p2g0.203 t/MWh这个系数一乘P2G的CO2需求量就和天然气产量直接挂钩了。还有一个小细节P2G不是随便什么时候启停都行的存在最小稳定出力一般取额定容量的20%左右。如果P2G容量50MW最小出力就是10MW。这可以用一个二元变量来控制或者如果算例日里P2G大多数时间都在满发或关停也可以简化成出力下限约束具体看你的场景复杂度。2.3 碳捕集能耗和碳流要线性配对碳捕集系统的建模重点是两个系数捕集率和能耗强度。捕集率指从烟气中实际捕集下来的CO2占排放量的比例我取0.85。能耗强度指每捕集1吨CO2消耗多少电胺法捕集一般每吨CO2耗电0.25到0.35MWh代码里取0.3MWh/t。在模型中CHP的碳排放量由燃料消耗量乘以排放系数得到。天然气排放系数大约0.21吨CO2/每MWh燃料热值。假设CHP每小时消耗212MWh天然气热值对应排放约44.5吨CO2CCS按0.85捕集率能捕到38吨CO2而捕集这38吨CO2需要消耗11.4MWh电。你看CCS本身就是一个耗电大户如果不把它放在电平衡约束里整个系统的电量计算就会对不上。捕集下来的CO2有三个去向送到P2G做甲烷化原料、封存到地下、或者直接外购/外售。我建模时把CO2流当成一种可分配的资源配平约束写成“捕集量 P2G消耗量 封存量”。如果P2G需要更多CO2而捕集量不够允许从外部购买CO2补充只不过价格稍高这样模型始终可行不至于因为碳源不足直接无解。2.4 储能设备与四条平衡约束储热罐模型相对标准储热量下一时段等于上一时段乘一个自损系数加上充热功率减去放热功率。充放热不能同时进行用big-M约束加二元变量控制。储气罐同理只不过充的是天然气单位统一用MWh热值方便和气网购气直接加减。四条平衡约束是整个模型的骨架。电平衡电负荷等于风电光伏出力加CHP电出力加储热放电对应的电功率折算加购电减去P2G耗电和CCS耗电。注意P2G和CCS都是耗电侧我习惯把它们放在等号右边作为负荷项这样不容易漏。热平衡热负荷等于CHP热出力加燃气锅炉热出力加储热放热减储热充热。气平衡气网购气加P2G产气加储气放气等于CHP耗气加燃气锅炉耗气加储气充气。CO2平衡则是2.3节讲的捕集量分配关系。模型内部的变量全部以MW和MWh为单位只有最终做结果展示的时候才按热值换算成m³。这个统一单位的习惯帮我避免了很多次量纲混乱强烈建议你也这么做。3. Matlab代码实现从参数初始化到求解器调用3.1 代码架构怎么组织我写这种优化项目习惯分四个文件主程序、参数文件、约束构建文件、结果后处理文件。主程序负责读参数、调用Yalmip建模、求解、调用后处理。参数文件里放所有设备参数和负荷曲线。约束构建文件专门写sdpvar变量声明和约束拼接。结果后处理文件负责计算成本明细、绘制条形图和曲线图。这样拆开的好处是改参数不用翻代码改约束不会误伤数据。整套模型升级成多场景或鲁棒优化时只需要在约束构建层加循环其他部分基本不用动。3.2 Yalmip建模的核心代码片段Yalmip建模的第一步是声明优化变量。我按设备分组声明方便检查。下面是一段核心片段T 24; p_gt sdpvar(1, T); % CHP电出力MW h_gt sdpvar(1, T); % CHP热出力MWth p_p2g sdpvar(1, T); % P2G输入电功率MW p_ccs sdpvar(1, T); % CCS捕集能耗MW c_cap sdpvar(1, T); % CCS捕集CO2量t/h y_p2g sdpvar(1, T); % 捕集CO2用于P2G的量t/h y_sto sdpvar(1, T); % 捕集CO2封存量t/h p_buy sdpvar(1, T); % 购电量MW q_buy sdpvar(1, T); % 购气量MWth h_st sdpvar(1, T 1); % 储热罐储热量MWh v_gas sdpvar(1, T 1); % 储气罐储气量MWh这里注意储热和储气要多声明一个时段因为动态约束要体现从t到t1的状态转移。Yalmip里变量维度是(1,T)写约束时可以直接用向量表达式也可以写for循环。我建议先写for循环把逻辑理顺跑通后再优化成向量写法避免一开始就维度出错。约束构建片段Constraints []; % 电功率平衡 Constraints [Constraints, p_load p_wind p_pv p_gt ... p_dis_heat p_buy - p_p2g - p_ccs - p_ch_heat]; % 热功率平衡 Constraints [Constraints, h_load h_gt h_gb h_dis - h_ch]; % P2G产出天然气的能量等于输入电功率乘效率 q_p2g eta_p2g * p_p2g; % CCS捕集量与耗电关系1吨CO2耗电0.3MWh Constraints [Constraints, p_ccs 0.3 * c_cap]; % CO2平衡 Constraints [Constraints, c_cap y_p2g y_sto]; Constraints [Constraints, y_p2g alpha_p2g * q_p2g];目标函数写成% 购电成本 购气成本 碳交易成本 弃风惩罚 obj sum(price_e .* p_buy) sum(price_g .* q_buy) ... carbon_price * (sum(e_gt e_gb) - co2_quota) ... penalty_wind * sum(p_wind_avail - p_wind_use);然后调用求解器ops sdpsettings(solver, gurobi, verbose, 2, gurobi.MIPGap, 0.005); diagnostics optimize(Constraints, obj, ops);3.3 求解后的结果核验跑通优化只是第一步最重要的其实是核验。我每次求解完都会写一段脚本把每一条平衡约束的残差算一遍看看最大残差是不是在10的负6次方量级。如果哪个平衡约束的残差明显不为零说明建模时某个变量的符号或者索引写错了这类错误靠肉眼看代码很难发现但残差计算会一抓一个准。举个核验电平衡的例子residual_e value(p_load) - (value(p_wind_use) value(p_pv) ... value(p_gt) value(p_dis_heat) value(p_buy) ... - value(p_p2g) - value(p_ccs) - value(p_ch_heat)); max(abs(residual_e))如果结果里有负值或者数量级异常就要回去检查对应边界约束。我以前遇到过P2G出力算出来是负数排查发现是p_p2g漏写了非负下界Yalmip默认变量自由上下界这种错误特别隐蔽。解决办法是每个有物理意义的变量都显式加上下界约束不要偷懒。4. 算例分析三种方案对比看清协同价值4.1 典型日场景与参数设置为了把问题说清楚我构造一个冬季典型日。冬季热负荷高CHP被迫高发风电在凌晨大发昼夜负荷差明显是最能体现P2G和CCS耦合价值的场景。风电装机200MW光伏装机100MWCHP额定电功率100MW、热功率90MWP2G容量50MWCCS捕集率0.85储气罐容量500MWh储热罐容量200MWh。电价采用分时电价峰时段10:00-15:00和18:00-22:00为0.8元/kWh平段0.5元/kWh低谷时段23:00-次日7:00为0.2元/kWh。气价取2.6元/m³折合每MWh天然气热值约268元。碳配额150吨/日碳价80元/吨弃风惩罚200元/MWh。这些参数不一定和你手头案例一样但量级是合理的。4.2 方案对比表格我跑三个方案Case1只有CHP和燃气锅炉Case2在Case1基础上加CCSCase3把CCS和P2G同时接入。结果如下方案日运行成本万元弃风率CO2排放吨/日P2G日产气m³/日Case1CHP基线46.221.3%1850Case2CHPCCS49.816.5%1230Case3CHPCCSP2G47.38.7%118186004.3 结果说明什么Case2最直观的教训是CO2捕下来了排放也确实降了但成本反而升高。原因很简单CCS本身要耗电捕集一吨CO2要烧掉0.3MWh电日捕集60吨就要多耗18MWh这部分负荷增加了购电压力尤其在峰时电价时段把成本拉上去了。所以CCS单独上马经济性并不好看。Case3把P2G接进来之后凌晨低谷时段原本弃掉的风电被P2G消化掉P2G发电需要CO2而CCS捕下来的CO2刚好直接喂给它省掉了外购CO2的费用也省掉了部分封存成本。生产出的18600m³天然气存入储气罐白天热负荷高的时候放出来给燃气锅炉用又抵消了一部分购气成本。这就实现了“低谷弃电→合成天然气→高峰供能”的跨时段转移成本比Case2低2.5万元/日弃风率从16.5%降到8.7%。更关键的是CO2排放。Case3的排放比Case2还低因为部分天然气产品被“回收碳”替代等效降低了系统净碳排放。单看P2G或者单看CCS都会觉得经济性一般但两者耦合起来成本和排放是双赢的这才是这个课题最有说服力的结论。5. 常见问题与排错经验5.1 求解器报Infeasible怎么办我踩过最多的坑就是模型不可行。Yalmip会返回diagnostic信息但通常不会告诉你具体哪条约束导致不可行。我的排查顺序是固定的第一步先把储热和储气的末状态约束去掉很多模型在周期闭环约束上翻车。第二步检查功率平衡约束的符号尤其是购电、购气变量是不是漏了非负约束。第三步用二分法思想先跑一个只有CHP和燃气锅炉的简化模型确认可行后再逐步加入CCS和P2G加到哪一步挂了问题就在哪一步。另外一个常见来源是P2G的CO2需求约束。如果CCS捕集量设得太小而P2G额定容量又大模型会为了满足CO2平衡被迫让P2G不能满发如果P2G的出力下界又限制了它不能低于10MW就会直接不可行。解决办法是把CO2外购变量加上或者允许P2G关停加二元变量。5.2 量纲不统一导致的“结果离谱”结果跑出来后最怕的不是报错而是结果不合理但没报错。比如P2G产气量算出来有几十万方明显超出常识十有八九是量纲问题。我在模型内部全部用MWh做能量单位气网购气价格也用元/MWh只有最后展示结果时才除以9.7转成m³。这样所有等式约束里的“气”本质上都是能量不会出现“电功率MW减气功率m³/h”这种荒唐组合。如果你拿到一套别人的代码第一件事就是检查每个变量的单位注释把所有常数和公式用维数分析过一遍。我见过有同学把0.203吨CO2/MWh和1.96kgCO2/m³混在一起用最后算出来的CO2需求量差了近一个数量级这种错误在审稿阶段是致命的。5.3 求解时间失控时怎么降维24时段的MILP模型通常几秒钟就能解完真正头疼的是全年8760小时或者双层的多场景模型。如果整数变量太多求解时间会暴涨。我常用的手段有三个一是把储热和储气同时充放的互补约束用线性不等式替代减少二元变量充放功率上限之和的约束有时就能覆盖需求二是给Gurobi设置MIPGap比如0.5%让它提前接受次优解三是先做典型日聚类全年问题聚成4到6个典型场景别直接上8760。另一个容易被忽视的是big-M的取值。储热充放约束里M取得太大会导致数值病态Gurobi求解稳定性下降取太小又可能把可行域切掉。我的习惯是先看储热罐最大充放功率是多少M取它的1.5倍左右然后加上非负约束模型基本不会因为big-M出问题。5.4 结果里出现负出力或资源不平衡有时候跑出来的解里CCS捕集量很高但P2G消耗量很低封存成本也高居不下这是正常的但不正常的是封存量算出来是负的。出现负值要么是变量缺下界要么是CO2平衡约束里被减项的位置搞反了。把每一个设备的变量下界都显式加上包括购电量、购气量、P2G耗电、CCS耗电、储热充放功率、储气充放功率基本上能防住这类问题。还有一种情况是储气罐的初始和末端储量不衔接。如果初始储量设40%末端要求不低于40%P2G产气多的时段储气就必须能存得下储气罐容量设小了会直接把模型锁死。先放开末端约束看解的变化再决定是改容量还是改约束别一上来就怀疑求解器。6. 个人实操中的几个深层体会6.1 建模阶段最值得花时间的三个点第一个是CHP可行域。固定热电比模型虽然能跑通但算出来的调度方案在实际机组上根本执行不了审稿人也一眼能看出来。多边形顶点凸组合法值得花半小时改进去。第二个是单位体系。把所有物理量统一成MW和MWh之后模型的通篇一致性会好很多后处理再转体积单位。第三个是结果核验脚本。我强烈建议在写主优化代码的同时写一个残差检查脚本每次跑完先看残差再看图这样能省下大量排查时间。我踩过的最深的一个坑是CCS耗电放在电平衡的左侧而CHP的净出力又没有扣除这部分导致晚上峰时段“凭空”多出一块电模型为了让平衡成立会让CHP多发电碳排放异常升高。这种错误从目标函数数值上很难发现但残差检查一跑就露馅了。6.2 后续可以继续扩展的方向这套模型跑通之后扩展方向其实很自然。你可以在确定性MILP的基础上加风电出力的多场景随机优化把不确定集和决策变量的关系处理成两阶段也可以做日前-日内两时间尺度的滚动调度日前定机组启停和储能计划日内修正出力还可以把热网动态和建筑热惯性引入进来让热负荷变成柔性负载此时P2G和CCS的调节空间会更大。最后再分享一个小技巧如果只是写论文算例数据不要追求极端好看。把CCS容量设得刚刚好P2G容量也选在能和弃风量匹配的量级结果才会有说服力。现实中我见过不少模型P2G容量给到200MW远超系统富余电能结果整天满发看着效果很好但一和实际工程对比就站不住脚。建模时候多想一步后面对你的工作评价会完全不同。