最近在整理综合能源系统优化调度的实验代码又把电转气、碳捕集和热电联产这三个模块完整搭了一遍。说实话这个题目看起来长核心其实就一句话把富余风电通过电解水制氢再与CO2甲烷化生成天然气同时把热电联产机组排放的CO2捕集回来作为甲烷化的原料形成一个“电-热-气-碳”相互耦合的小闭环然后用Matlab把这套模型写成可求解的优化代码。这篇文章我会从建模动机、数学建模、Matlab实现、算例设置到调试坑位按我自己做项目的顺序完整梳理一遍。如果你正在做综合能源系统方向的开题、写论文或者复现代码应该能直接省掉不少试错时间。这套代码我已经在不同参数下反复跑过也踩了不少坑下面写的都是实际可以落地的内容。1. 为什么要把电转气、碳捕集和热电联产放进同一个优化模型先说一个很现实的背景在我国北方很多地区冬天供暖季的热负荷非常高而热电联产机组因为“以热定电”的运行方式夜间电出力降不下来。这时候偏偏是风电大发时段电网消纳空间被CHP挤占弃风就成了家常便饭。过去我们习惯加电锅炉让多余的电变成热来缓解弃风但这种方式本质上是“把电能降级成热能”只解决了部分问题并没有给系统提供真正的灵活性。电转气Power to GasP2G的思路不一样它是把多余的电能转化为氢气或天然气把电能“升级”成可存储、可运输、可回燃的气体燃料。热电联产机组之所以在夜间压缩风电空间是因为它的电出力由热出力牵制着P2G接入后相当于给CHP的刚性电出力开了一个侧路风电不用再和机组硬挤低谷电量可以被消化成气体储存起来到高峰时段再释放。碳捕集Carbon Capture SystemCCS则从另一个维度切入它把CHP烟气中的CO2分离出来大大降低了系统碳排放。单独看CCS是一个减排设备但放在P2G系统里它捕集的CO2刚好是甲烷化反应的碳源。就这样三条看似独立的技术路线在同一个模型里形成了互补闭环。1.1 “以热定电”这个老问题为什么P2G能撬动“以热定电”背后的物理逻辑是热电联产机组必须保证供热而抽汽供热后汽轮机还要维持一定的发电量。于是热负荷越高电出力越压不下去。夜间风电大发时系统里没有足够的负荷去消纳这些电只能弃风。P2G接入后低谷时段的电能可以被电解槽吃掉变成氢气和后续的天然气。关键在于P2G的耗电是灵活可调的它像一个“弹性负荷”专门在低谷时段开启在高峰时段关小。这样CHP不必为了跟随风电而大幅波动风电出力也找到了去处。如果再配置储气罐P2G产出的气体可以存到第二天高峰再燃烧供热等于实现了电能在时间维度上的转移。为什么模型里选择“电转甲烷”而不是“电转氢”我的理由是如果只产氢氢气很难直接并入现有燃气管网和燃气轮机因为设备对氢气比例有严格限制。甲烷化的产物是合成天然气可以和常规天然气一样被CHP直接利用系统耦合最顺畅。这也意味着CO2成了一个不可或缺的原料于是自然引出了碳捕集系统的接入需求。1.2 碳捕集在这里不只是减排设备还是甲烷化的原料供应商单独看碳捕集它的作用可以用一句话概括从烟气中把CO2分离出来不让它排到大气里。但在P2GCCS耦合模型里CCS的意义远不止减排。甲烷化反应的化学式是 CO2 4H2 → CH4 2H2O没有CO2这个反应根本走不了。如果从空气中直接捕集CO2成本高得离谱而从自己系统里的CHP烟气中捕集等于就地取材。于是整个系统形成了一个有趣的碳循环CHP烧天然气排放CO2 → CCS从烟气中捕集CO2 → P2G用氢和CO2合成甲烷 → 甲烷又被CHP烧掉。外排到大气的CO2只剩下捕集率没有覆盖到的那一小部分。建模时如果忽略这条碳流CCS就只是个昂贵的减排末端P2G也缺一个便宜的碳源二者叠加的优势完全体现不出来。这里要注意CCS本身是耗电大户。胺法捕集一吨CO2大约需要0.2到0.4 MWh的电耗捕集量越大厂用电率越高系统电平衡里的“负荷侧”就越重。所以建模时不能把CCS当成一个单纯的减排约束必须把它用电量显式写进电平衡方程。我在代码里会把捕集电耗写成捕集量的线性函数单独占一个变量这样结果才真实。1.3 热电联产机组是电热碳三条流交汇的核心节点在电、热、气、碳四条流里CHP是交叉点它同时生产电和热消耗天然气排放CO2。很多初学者的模型把CHP简化成“电出力效率×燃料”再乘一个热电比得到热出力这其实只适用于背压式机组。更常见的抽汽凝汽式机组电出力和热出力之间是一个二维可行运行域热出力增加时电出力上限会下降且电出力下限还会随热出力改变。所以CHP的建模质量决定了整个优化结果能不能落地。如果可行运行域写得不准确优化器可能给出一个在实际机组上根本跑不出来的运行点后面无论怎么算都是白搭。后文我会专门说明CHP可行运行域在Matlab里的写法。这三个小节合起来已经回答了一个问题为什么这个题目要把三者放在同一个模型里因为风电消纳、碳减排、热电供应这三个矛盾单独靠任何一块都解决不干净只有让P2G消纳电力、CCS回收碳、CHP提供热电基础同时把三条流耦合起来才能得到一个既降碳又降弃风的系统方案。2. 建模框架从物理过程到可求解的数学表达2.1 我习惯用“母线式能量平衡”组织约束我最早做这类模型时喜欢一台设备写一组约束结果十几个设备下来约束又长又乱查错特别痛苦。后来换成了母线式写法把电网、气网、热网、碳流分别当成虚拟母线每台设备只是从某条母线取功率或者向某条母线注入功率。代码里只需要写四条平衡方程所有设备约束再单独列。电平衡风电出力 CHP电出力 购电 电负荷 P2G耗电 CCS捕集耗电 其他电耗。热平衡CHP热出力 锅炉热出力 热负荷。气平衡外购气 P2G产甲烷 储气罐放气 CHP燃气消耗 储气罐充气。碳平衡CHP燃烧产生的CO2 直接排放部分 CCS捕集部分CCS捕集到的CO2 P2G甲烷化消耗 封存或外供部分。这样组织的好处很明显每一行都有明确的物理意义新增设备时只需要在对应母线上增减一个变量不会牵一发动全身。后面我给的Matlab核心代码也是按这个思路写的。调试时如果结果不平衡直接逐个母线检查残差很容易定位是哪条约束漏了项。2.2 CHP机组的可行运行域怎么表达CHP的可行运行域是我觉得全模型最需要谨慎处理的地方。一个相对完整的抽汽凝汽式机组模型电出力和热出力关系不是一条直线而是一个由多个线性不等式围成的凸多边形。我算例里用的是简化三约束0 P_chp 300; % 电出力上下限MW 0 H_chp 350; % 热出力上下限MW P_chp 0.15 * H_chp 300; % 电热耦合上限 P_chp - 0.45 * H_chp 100; % 最小电出力随热出力变化这里第一、二条是设备本身的物理边界第三条表示抽汽量越大凝汽发电能力越小第四条表示热出力较低时机组仍有一定的最低技术出力。实际工程中你可以把厂商工况图离散成若干拐点再用凸包方法生成一组线性不等式原理完全一致。再加上爬坡约束-50 P_chp(t) - P_chp(t-1) 50; -80 H_chp(t) - H_chp(t-1) 80;CHP部分就算基本成型了。有一点要提醒如果原机组的可行运行域不是凸的直接做MILP会非常难求实际项目里一般会先做凸近似损失一点精度但换来可靠的可解性。2.3 CCS和P2G的“碳-电”双向耦合关系P2G我按“电解水甲烷化”全流程建模。电解水环节的关键参数是单位电耗工程上大约4.5到5.5 kWh可以产1 Nm3氢气效率大致在60%到80%。甲烷化环节需要氢气和CO2H2和CO2的摩尔比接近4:1也就是按体积算大约4份氢气配1份CO2。为了方便优化求解我会把化学反应细节抽象成两组线性关系氢气产量由P2G电功率和电解效率决定甲烷产量由氢气输入量和甲烷化效率决定同时甲烷产量也对应一个确定的CO2消耗量。CCS侧的模型则是捕集到的CO2量等于CHP烟气排放量乘捕集率捕集能耗正比于捕集量。F_co2_cap eta_capture * E_fuel_emission; P_ccs_energy alpha_ccs * F_co2_cap;这里E_fuel_emission取决于CHP的天然气消耗量和排放因子。捕集到的CO2是去甲烷化还是去封存由优化器根据碳价和气价自行选择。这样做的好处是碳流从CHP烟气走到CCS再从CCS走到甲烷化甲烷再回到CHP气源整个闭环在数学上完整闭合。如果少了这条碳平衡优化器很容易凭空产生CO2产甲烷量就会虚高。2.4 储气罐让时间耦合真正落地如果P2G产出的甲烷必须当场用完耦合效果会大打折扣因为低谷时段产的气没法留到高峰。实际系统中应该有一个储气罐它承担了时间解耦的作用。储气罐的状态转移约束很简单S_gas(t1) S_gas(t) G_in(t) - G_out(t) - loss; 0 S_gas(t) S_max; S_gas(1) S_gas(25);最容易被忽略的是初始时刻和最终时刻储气量一致。如果不加这个日周期边界条件优化器会把储气罐当成免费的垃圾场低谷大量充气、高峰全部放空最后一天的末时刻储量为零第二天还得从零开始加气这不符合连续运行逻辑。我在所有算例里都把S(1)S(T)设成硬约束这样结果才是真正可循环的调度策略。3. 优化目标、约束陷阱与线性化处理3.1 目标函数并不等于“总成本最小”那么简单我常用的优化目标都是系统日运行成本最小但每一项怎么计价需要想清楚。常见的成本项包括CHP燃料成本天然气消耗量乘以天然气价购电成本峰谷电价不同时用分时电价向量逐个时段乘P2G运行维护成本按输入电功率计也可以按产气量计CCS运行维护成本按捕集CO2量计弃风惩罚按弃风电量乘以惩罚单价碳交易成本实际排放量与免费配额之差乘以碳价。有一个很隐蔽的坑碳捕集虽然降低了系统外排CO2但捕集消耗的电如果来自主网主网侧可能存在间接排放。如果你的系统边界只到园区或微网内部就按系统内CHP排放计算如果要把主网纳入考核就必须给购电加一个碳排放因子。我在算例里为了让模型聚焦在局部耦合上只计算系统内部排放但在分析结论时会把边界条件写清楚避免读者误读。3.2 耦合约束是模型的灵魂也是最容易出现冗余的地方构建约束时我会反复检查四个不变量电平衡闭合、热平衡闭合、气平衡闭合、碳平衡闭合。比这更隐蔽的是设备连接关系约束。比如P2G的耗电必须来自同一个电母线不能单独指定一个“虚拟电源”CHP的燃气可以是外购气也可以是本地P2G产甲烷但不能同时让储气罐充气和放气CCS捕集的CO2参与甲烷化时必须先证明这部分CO2确实来自CHP烟气。我检查约束是否正确的方法很简单把目标函数的所有成本系数都设成0只求一个可行解。如果连可行解都找不到再逐个放宽约束看看是哪条母线的平衡被打破。在YALMIP里用optimize(Constraints,0)就能快速拿到诊断信息。这个方法比盯着公式硬想要高效得多。3.3 非线性项的线性化三个实用技巧很多优化问题之所以从非线性变成MILP是因为求解器对MILP的处理已经非常成熟。我有三个常用技巧第一P2G效率曲线分段线性化。电解槽在低功率段效率很低高功率段效率接近饱和可以用三到四段线性折线近似。YALMIP里可以用implies或iff描述折线段的激活逻辑但求解性能不如手动标准建模所以我一般自己引入0-1变量和辅助连续变量。第二Big-M处理启停逻辑。比如判断CHP是否开机z_chp binvar(1,24); P_chp 300 * z_chp; P_chp 20 * z_chp;M值不要取1e6紧凑的界就是300最多取到360。后面会专门讲为什么M值太大会导致数值灾难。第三惩罚函数线性化。如果目标里出现弃风惩罚的平方项模型会变成QP加入0-1变量后变成MIQP求解速度明显变慢。工程上我更推荐用分段线性罚函数替代二次惩罚即在弃风超过某个阈值后再增加单位惩罚效果和二次罚函数接近但求解效率高很多。4. Matlab代码实现从数据结构到求解器调用4.1 我常用的代码目录结构一套能反复修改的代码目录结构最好这样组织CHP_P2G_CCS/ 00_data/ base_profile.m unit_parameters.m 01_model/ build_balance.m build_chp.m build_p2g_ccs.m build_storage.m 02_solve/ run_case.m 03_result/ plot_results.m我坚持模块化而不是一个几百行脚本跑到底原因有两个一是换数据时只需要改00_data目录二是做多个对照组时不需要复制整个主程序只需要切换模型构建函数。参数文件里我通常用结构体保存所有设备参数unit.P_CHP_max 300; unit.P_CHP_min 20; unit.H_CHP_max 350; unit.coef_p_h 0.15; unit.eta_p2g_elect 0.65; unit.alpha_ccs 0.32; % MWh/tCO2结构体字段名清楚约束里写unit.P_CHP_max比裸数字可读性好得多。时间长了你会发现后续做灵敏度分析时只需要在循环里改结构体字段再重跑求解函数非常方便。4.2 核心建模代码从sdpvar到optimize下面这段是简化版的核心建模逻辑我用的是YALMIPGurobi。变量按母线分类约束逐个模块追加%% 变量定义 P_wind_use sdpvar(1, 24); % 实际消纳风电 P_chp sdpvar(1, 24); % CHP电出力 H_chp sdpvar(1, 24); % CHP热出力 P_grid sdpvar(1, 24); % 网购电 P_p2g sdpvar(1, 24); % P2G电耗 F_co2_cap sdpvar(1, 24); % CCS捕集CO2量 F_co2_use sdpvar(1, 24); % 用于甲烷化的CO2量 F_ch4 sdpvar(1, 24); % P2G产甲烷量 G_sto_in sdpvar(1, 24); G_sto_out sdpvar(1, 24); S_gas sdpvar(1, 25); % 储气量含初始时刻 C []; %% 电母线平衡 C [C, P_wind_use P_chp P_grid ... data.P_load P_p2g unit.alpha_ccs * F_co2_cap]; %% 热母线平衡 C [C, H_chp data.H_load]; %% CHP可行运行域 C [C, unit.P_CHP_min P_chp unit.P_CHP_max]; C [C, 0 H_chp unit.H_CHP_max]; C [C, P_chp unit.coef_p_h * H_chp unit.P_CHP_max]; C [C, P_chp - unit.coef2_p_h * H_chp unit.P_CHP_min]; %% P2G与CCS耦合 C [C, 0 P_p2g unit.P2G_max]; C [C, F_ch4 unit.eta_meth * P_p2g / unit.hv_ch4]; C [C, F_co2_use unit.CO2_per_CH4 * F_ch4]; C [C, F_co2_use F_co2_cap]; C [C, 0 F_co2_cap unit.CCS_max]; %% 储气罐 C [C, S_gas(1) unit.S_init]; C [C, S_gas(25) unit.S_init]; C [C, 0 S_gas unit.S_max]; C [C, S_gas(2:25) S_gas(1:24) G_sto_in - G_sto_out ... - unit.loss * S_gas(1:24)]; %% 目标函数 Cost sum(unit.price_gas * F_gas) sum(unit.price_grid .* P_grid) ... sum(unit.om_p2g * P_p2g) sum(unit.om_ccs * F_co2_cap) ... sum(unit.penalty_wind * (data.wind_avail - P_wind_use)) ... unit.price_CO2 * (sum(E_emit) - data.CO2_quota); %% 求解 ops sdpsettings(solver,gurobi,verbose,2,... gurobi.MIPGap,1e-4,... gurobi.TimeLimit,300); optimize(C, Cost, ops);这段代码为了演示已经做了大量简化实际跑算例时还需要补F_gas和E_emit的具体表达式但骨架就是变量、约束、目标、求解四步。初学者最容易卡住的是变量类型定义CHP的出力是连续变量只有启停状态才用binvar储气量也是连续变量不需要离散化。另外data.wind_avail - P_wind_use就是弃风电量因为这里风电出力上限是已知曲线实际消纳不能超过它。4.3 数据准备中容易被忽略的单位和边界我踩过的最大的坑是单位换算。论文里碳捕集电耗常用kWh/kg CO2产气量常用Nm3/h热负荷用GJ/h如果不统一模型结果差一个数量级还完全看不出。我的做法是全部统一到下面这套单位制电功率用MW热功率用MW气体量按热值折成MWhCO2量用t。比如1 Nm3甲烷的低位热值约35.8 MJ相当于0.00994 MWh我会在参数文件里写一个换算函数function q_mwh ncm3_to_mwh(V_nm3) q_mwh V_nm3 * 35.8 / 3600; end外部数据进来后全部先过一遍换算函数能有效杜绝热值单位MJ和MWh混用的经典错误。另外一天24个小时功率乘以1小时就是能量所以代码里所有sum(P)自动就是日电量MWh如果后续改成15分钟一个时段T96所有与时间有关的系数都要除以4。这个细节常常决定结果能不能对上工程实际。5. 算例设计与结果讨论耦合方案到底带来了什么5.1 三组对照方案的设置为了把P2G和CCS各自的作用拆开我设计了三个对照方案。所有方案都用同一天的风电、电负荷和热负荷曲线机组参数也保持一致唯一区别是碳捕集和电转气是否接入。这样后面对比结果时指标差异就能明确归因于对应子系统。方案配置说明ACHP 风电不加任何低碳子系统作为基准BCHP 风电 P2G只用电转气消纳风电碳捕集不参与CCHP 风电 P2G CCS完整的碳-气耦合闭环典型日的设置上我让夜间风电高出、电负荷峰谷明显、热负荷全天保持高位这样才能暴露“以热定电”和弃风矛盾。5.2 结果数据与基本结论先说明一点不同论文的参数和基准选择会造成结果差异下面给的是用合理参数跑出来的趋势性结果具体数值只作为量级参考。我的一轮典型结果如下方案日运行成本碳排放弃风率A36.2万元214吨18.7%B34.8万元198吨6.2%C33.5万元151吨3.4%趋势很清楚P2G单独接入后主要收益是降低弃风但对碳排放的直接影响有限因为多余风电变成甲烷后可能又被CHP烧掉整体碳循环没有闭合。CCS接入后CHP烟气中的CO2被重新送回燃料流碳排放大幅下降同时甲烷化有了廉价碳源P2G的经济性进一步提升弃风也因此压得更低。这就是“耦合”比“单项叠加”更有意义的原因两个子系统的收益不是简单的加法而是相互促进的。5.3 碳价和风电渗透率变化的灵敏度规律我随后把碳价从50元/t抬到300元/t发现方案C的优势越来越大。高碳价下CCS捕集CO2后参与P2G甲烷化等于用本地碳源替代外购天然气外购气成本下降减排还能获得碳收益。方案B在高碳价下也能增加产气量但由于没有捕集它只能消纳风电无法降低CHP自身排放碳成本还是会走高。所以碳价越高CCS和P2G的耦合价值越明显。风电渗透率提高时弃风风险也会变大。方案B和C都会增加P2G电功率但受到P2G容量和甲烷化CO2需求的限制。如果CCS捕集能力不够P2G即便有电也不能无限产气这时碳约束变成了系统瓶颈。实际项目中如果风资源特别丰沛可能需要配置更大CCS容量或者引入外部CO2源来补足碳供应。5.4 结果里两个容易误读的指标“弃风率下降等于碳排放下降”这句话并不一定成立。如果P2G产甲烷只是替代了外购天然气而CHP总产热量不变整体碳排放变化不大只有CCS接入、外排CO2被重新利用时碳排才会明显下降。所以解读结果时最好把弃风率和碳捕集量分开看。另外“总成本下降”可能来自很高的弃风惩罚系数。如果惩罚单价取得比上网电价还高优化器会不计成本消纳风电得出的成本结论会有偏差。我建议在结果输出里单独列出惩罚项金额然后做一次惩罚系数敏感性验证看看结论是否稳定。6. 复现这套模型最容易踩的坑6.1 时间尺度不一致模型和物理脱节电平衡本质上是实时的但储气罐和燃气管网是慢过程。如果模型只有24个点且没有储气罐P2G产出的甲烷只能当下被消耗等于是把“气体储能”排除了优化结果会严重低估P2G的灵活性。反过来如果加了储气罐但忘记初始储量和日末约束优化器会把储气罐当成免费垃圾场大量在低谷充气、高峰放气结果看似漂亮实际跑不出来。我的处理方式是一定保留储气罐状态变量并把S(1)S(25)写成硬约束让调度成为一个完整日循环。如果研究对象是几天或一周就要把日末约束改成周期末约束同时加入储气罐初始水平的历史数据。6.2 Big-M取值不合理求解器性能天差地别很多复现者喜欢写M1e6觉得这样一定不会出错。但在Gurobi里大M会引入大量数值误差有些约束明明满足却因为容差被判违反求解速度也会急剧下降。我的原则是M取到物理边界的1.2倍左右。比如CHP电出力上限300 MW就用300或360储气量上限100 MWh就用120。YALMIP里如果写了implies或iff会自动引入Big-M建议查一下翻译后的约束必要时手动建模替换。6.3 数值缩放和求解容差导致的“随机差异”同一套代码在不同电脑上结果不同最常见原因不是随机化而是单位量级差太远。比如捕集CO2量只有几十吨电功率却是几百MW目标函数里几个数量级差太多的项会让求解器数值病态不同MIP gap下给出的取舍也完全不同。解决办法是在数据准备阶段做统一缩放把数量级控制在1e-2到1e2之间。例如碳价50元/t日成本几十万元不妨把成本单位改成万元排放量用百吨为单位数值问题会少很多。求解器参数方面我会设置合理的MIPGap和时间限制ops sdpsettings(solver,gurobi,verbose,2,... gurobi.MIPGap,1e-3,... gurobi.TimeLimit,300);Matlab里我还会在求解前做一次模型诊断[detail, infeas] analyze(Constraints, Cost);新手拿到无解报错后最忌讳直接删约束“试到能跑为止”那样会得到物理不完整的模型。正确做法是用analyze定位是哪条约束导致无解再针对性修正参数或线性化方式。7. 从这套代码还能往哪些方向扩展7.1 把确定性优化改成两阶段鲁棒或随机优化现在的模型里风电是给定曲线属于确定性优化。真实调度中风电预测误差很大我建议下一步引入场景生成和削减。比如用历史风电数据生成100个场景再同步回代削减成10个典型场景对应预测场景和极端场景。目标函数可以改成期望成本约束里加上第一阶段的日前决策和第二阶段的实时调整变量。YALMIP处理这种规模仍然可行但要控制场景数否则求解时间会暴涨。两阶段鲁棒优化的更成熟做法是结合对偶和CCG算法工程量会大不少但结果更有说服力。7.2 把模型往实际园区项目搬如果这套模型要用于真实园区还需要补充网络潮流约束尤其是气网压力和天然气组分以及设备的部分负荷效率曲线。P2G装置的启动时间、最低负荷率、响应速度都要加进来。我的建议是先从“能量平衡机组容量约束”的MILP版本跑通拿到结果后再用MatlabSimulink做多时段动态仿真验证。优化证明趋势仿真验证动态两条腿走路最稳。7.3 一些个人实操体会我自己调试这类模型时最管用的习惯是每加一个模块就先跑一次可行解不要让多个新约束一起进入模型。因为无解时你根本不知道是碳平衡写错还是储气罐约束写错。另外结果文件里一定要输出每条母线的残差比如电不平衡量、热不平衡量是多少。残差接近0模型基本正确残差不是0但求解器说最优那一定是约束漏项了。很多论文里画的能量流图如果你自己从优化结果里按母线累加会发现根本对不上这往往就是约束缺失。如果以后你要在这个基础上做实时调度我建议再把P2G和CCS的启停成本写进目标函数。电转气设备频繁启停成本很高确定性日调度里看不出来一旦放到滚动优化里就原形毕露。到那时可以把P2G的爬坡约束改成小时级最小连续运行时间约束效果会好不少。先跑通再一点一点加约束永远比一口气建复杂模型舒服得多。