做综合能源系统调度优化的人对这套组合应该不陌生热电联产机组CHP、电转气设备P2G、碳捕集系统CCS。单个拎出来都是老话题但要把三者放进同一个优化模型里在Matlab中实现一套低碳经济调度事情就完全不一样了。CCS的捕集能耗会重新分配CHP的电热出力P2G既要吃电又要吃CO2还得跟储碳罐的库存联动热电联产那套经典的以热定电约束在这种结构下基本要重写。这篇文章记录了我复现含P2G与碳捕集系统的热电联产建模与优化的完整过程从系统架构、数学建模到Matlab代码和算例结果最后把我在YALMIPCPLEX求解中踩过的坑也一并倒出来。适合正在做综合能源系统方向研究的学生、准备给园区能源系统做低碳改造的工程师以及想拿一个真实优化案例练手Matlab建模的朋友。1. 从排碳-耗能到柔性化生产P2G与碳捕集耦合的系统逻辑1.1 这套系统到底解决什么问题传统热电联产机组的核心矛盾在于热和电被强耦合在一起。冬天供暖需求大时机组必须多供热连带发出来的电远超实际用电需求多余的电在风光大发时段直接变成弃电。更要命的是烧煤烧气必然排碳在碳约束越来越紧的背景下老师傅们都在琢磨怎么让CHP既能灵活调节出力又能把碳排放压下去。P2G和CCS恰好是从两个方向补这个短板。CCS把烟气里的CO2抓下来不让它直接排到大气里P2G把多余电能转成天然气相当于给电网装了一个可以随时调节大小的电胃口。如果把两者连起来CCS捕集到的CO2恰好是P2G甲烷化反应需要的原料这就形成了一条CHP排碳到CCS捕碳再到P2G用碳的内部循环链。听起来很完美但工程建模的复杂度就在这里——它们不是三个独立设备而是一个互相牵制的整体。我做复现时第一个顿悟就是这个题目表面是加设备实际是改变整个系统的调度自由度。P2G充当柔性负荷去吸收波动性可再生能源CCS充当碳的调节阀两者配合后CHP的热电耦合约束被大大放松了。理解到这一层后面建模才有方向。1.2 能量流与碳流怎么走全局耦合关系拆解要建模先把系统中的物质流和能量流捋清楚。整个系统站在电网、气网、热网三条母线的交汇点上看电网上CHP和风电场发电负荷用电P2G电解槽用电CCS的吸收剂再生也用电。气网上P2G产出的甲烷注入气网或供给燃气负荷。热网上CHP供热给热负荷同时分出一部分蒸汽给CCS的再生塔做热源。碳流是另一条线CHP燃烧产生的烟气进入CCS吸收塔富液送到再生塔加热解吸出高浓度CO2一部分CO2去封存或外售另一部分送到P2G的甲烷化反应器与电解水产生的H2合成CH4。我把这些关系画成表格式的对应关系建模时就不会漏约束设备输入输出中间产物CHP机组燃料煤/气电、热、烟气CO2-CCS系统烟气、电、热浓CO2、净烟气吸收富液P2G电解槽电、水H2-P2G甲烷化H2、CO2CH4、H2O-注意一个细节P2G两步反应里甲烷化这一步才消耗CO2电解水制氢本身不需要碳。所以建模时如果只写P2G耗电产气就丢掉了CO2这个中间纽带。我建议要么把P2G拆成电解和甲烷化两段要么用化学计量比把CO2消耗量直接折算成产气量的线性函数下文会展开。1.3 耦合之后传统CHP调度模型哪里不够用了我最早学CHP调度时用的是最朴素的背压式模型热出力等于热电比乘电出力H c_m·P。一个等式就把问题定死了机组只能沿着一条线运行调度起来非常省心。但加了CCS之后这条线必须松绑。原因很简单当CCS投入运行CHP的净上网电功率不再是P本身而是P减去CCS消耗的电功率。换句话说同一台机组在同一个热出力下因为碳捕集负荷不同对外表现的净电出力也不同。这样一来热和电的解耦不再靠抽凝式机组的物理结构而是靠碳捕集装置这个可调负荷来腾挪。P2G的影响则体现在电功率平衡上。以前风电多了只能弃掉因为常规机组压不下去、CHP又受供热限制压不下去现在有了P2G它可以像一个大功率电锅炉一样把多余电能吃掉转化天然气。这个吃电的能力是可调的调度模型里它就成为一个决策变量和CHP、风电一起做联合优化。所以传统模型不够用的点有三个一是CHP出力可行域需要重写二是碳捕集能耗必须作为与捕集量耦合的变量进入平衡方程三是电平衡里多了P2G和CCS两个大功率用电项它们不是常数是由优化决定的。这三条正是后文数学模型的核心。2. 数学建模的关键环节CHP、CCS、P2G各自的约束长什么样2.1 热电联产机组可行域要比以热定电复杂一点实际研究里抽凝式CHP比背压式更常用因为可调范围大。抽凝式机组的电热可行域可以近似描述为一个凸多边形用一组线性不等式表示电出力上下限P_min ≤ P_chp ≤ P_max热出力范围0 ≤ H_chp ≤ H_max热电耦合约束P_chp ≥ P_min β₁·H_chp保证供热足够时电出力不低于下限最大出力受限P_chp ≤ P_max - β₂·H_chp供热抽汽多凝汽发电少背压式则直接简化为H_chp c_m·P_chp适合做机理分析但做调度优化时我首选抽凝式因为系统自由度更大能看出P2G和CCS带来的调节价值。燃料消耗与电热出力之间我习惯用线性函数逼近F_chp a₀ a₁·P_chp a₂·H_chp。单位是kW或MW的燃料功率。线性化肯定有误差但对24小时日前调度来说足够而且可以避免二次规划带来的求解负担。如果机组数据里有明确的燃耗曲线系数直接用最小二乘拟合出a₀、a₁、a₂即可。2.2 碳捕集系统捕集量、净排放与辅助能耗的折算CCS建模最核心的是三个量总排放量、捕集量、净排放量再加上捕集能耗。CHP产生的总CO2量与燃料消耗成正比简化写成E_total e_int·P_chpe_int是单位电出力的碳排放强度。捕集系统投入运行时设t_c为捕集率则捕集量为E_cap t_c·E_total净排放为E_net E_total - E_cap。捕集过程本身要耗电耗热。以目前最成熟的燃烧后化学吸收法为例典型的再生热耗约3~4 GJ/t CO2电耗约100~200 kWh/t。把这些折算成与捕集量成比例的辅助负荷P_ccs λ_el·E_cap电耗H_ccs λ_heat·E_cap热耗来自CHP抽汽这两项会分别进电平衡和热平衡。我计算时通常把E_cap的单位统一为t/hλ_el取0.15 MWh/tλ_heat取1.0 MWh/t数值上对应150 kWh/t电耗和3.6 GJ/t热耗和公开文献数据基本吻合。还有一个容易被忽略的约束捕集系统不要面面俱到烟气可以直接旁路。也就是说捕集率t_c不必固定而是一个可调节的变量可以理解为运行在0到上限之间的连续值。这给了调度系统一个新的控制自由度——碳价高时多捕碳价低时少捕和机组出力一起优化。这个可调捕集率是整个模型灵活性的关键。2.3 电转气系统从电解水到甲烷化的化学计量约束P2G分两步电解水制氢再加氢甲烷化。总反应可以写成2H₂O → 2H₂ O₂电解 CO₂ 4H₂ → CH₄ 2H₂O甲烷化理论上1 mol CH4需要1 mol CO2和4 mol H2。建模时我把两段合并成一步法模型已知输入电功率P_p2g总效率η_p2g电解效率乘甲烷化效率范围大约45%~60%则产甲烷功率G_ch4 η_p2g·P_p2g这里G_ch4按热值计单位是MW。如果要同时看出气体积需要除以LHV约0.00994 MWh/m³即P2G产气体积约等于G_ch4/0.00994。按化学计量折算CO2消耗量产生1 m³CH4大约需要1.96~2.0 kg CO2。如果产甲烷功率为G_ch4MW则耗碳量可以写成C_co2 β_co2·G_ch4β_co2按上述折算大约0.2 t/MWh因为1 MW功率一小时产气约100 Nm³耗碳约200 kg。这个线性系数非常实用把一个化学反应约束变成了一个简单的比例关系YALMIP里直接写一行约束。如果追求更精细可以把P2G拆成两段建模电解段P_e2h到H2的转换甲烷化段由H2和CO2生成CH4后者受限于CO2供应量。但在24小时调度模型里一步法精度已经足够关键是别漏了P2G耗碳这一项。我见过不少实现里把P2G建模成单纯的电力负荷完全不写CO2消耗约束那等于把系统的碳流逻辑丢了。2.4 系统级功率平衡与储碳罐动态设备级约束凑齐后靠系统平衡把大家绑在一起。我用了四个平衡电功率平衡 P_chp P_wind P_load P_p2g P_ccs热功率平衡 H_chp H_load H_ccs气网平衡P2G产气全部外送气网作为一个次级模块处理不承担系统内燃气负荷时只记录产气量收益。碳平衡捕集量、P2G消耗量、储碳罐充放必须闭合。储碳罐是容易被忽略的元件。为了应对P2G和CCS运行节奏不一致系统里会有一个CO2缓冲罐动态约束为S_co2(t1) S_co2(t) E_cap(t) - C_co2(t)0 ≤ S_co2 ≤ S_max且S_co2(1) S_co2(T1)保证一天内碳库存回归初始值。这个约束在Matlab里用sdpvar定义一个T1维变量就能实现YALMIP天然支持时间索引写起来很顺手。储碳罐让CCS和P2G不必实时匹配CCS多捕的碳可以先存着等P2G有电可吃时再消耗大大增加系统灵活性。3. Matlab实现路线从目标函数到YALMIPCPLEX求解3.1 目标函数选择运行成本、碳交易与弃风惩罚的权重优化目标我选择低碳经济调度把运行成本和碳排放在同一个目标里权衡。具体包含四部分CHP燃料成本Σ c_coal·F_chp(t)c_coal按煤价折算碳交易成本Σ c_co2·E_net(t)碳价设为100元/t。净排放越多成本越高弃风惩罚Σ c_wind·(P_wind_avail(t) - P_wind(t))让优化器尽量消化风电不消纳就罚钱P2G售气收益-Σ π_gas·G_ch4(t)/LHV产出的甲烷按气价出售这是负成本项整体目标函数就是cost 燃料成本 碳成本 - 售气收益 弃风惩罚为什么弃风要放目标函数而不是硬约束因为如果硬性要求风电全消纳在某些时段CHP受热负荷限制压不下去时模型会无解。用惩罚项代替硬约束既保证可行又有经济解释这是调度建模里很实用的处理技巧。3.2 核心代码骨架变量定义与约束组装的规范写法我的Matlab环境是MATLAB 2023b加YALMIP加CPLEX。先定义决策变量T 24; P_chp sdpvar(1, T, full); % CHP电出力 H_chp sdpvar(1, T, full); % CHP热出力 P_wind sdpvar(1, T, full); % 风电实际并网功率 P_p2g sdpvar(1, T, full); % P2G耗电 G_ch4 sdpvar(1, T, full); % 产甲烷功率 P_ccs sdpvar(1, T, full); % CCS电耗 E_cap sdpvar(1, T, full); % 捕集CO2量 C_co2 sdpvar(1, T, full); % P2G消耗CO2量 S_co2 sdpvar(1, T1, full); % 储碳量约束组装的核心思路是逐条追加到Constraints变量里最后一次性丢给求解器Constraints []; % 电平衡 Constraints [Constraints, P_chp P_wind P_load P_p2g P_ccs]; % 热平衡CCS热耗从CHP热出力中扣除 Constraints [Constraints, H_chp H_load H_ccs]; % CHP抽凝式可行域 Constraints [Constraints, P_min P_chp P_max]; Constraints [Constraints, 0 H_chp H_max]; Constraints [Constraints, P_chp P_min beta1 * H_chp]; Constraints [Constraints, P_chp P_max - beta2 * H_chp]; % CCS捕集与能耗 Constraints [Constraints, E_cap t_c_max * e_int * P_chp]; Constraints [Constraints, P_ccs lambda_el * E_cap]; Constraints [Constraints, H_ccs lambda_heat * E_cap]; % P2G和耗碳 Constraints [Constraints, G_ch4 eta_p2g * P_p2g]; Constraints [Constraints, C_co2 beta_co2 * G_ch4]; % 储碳罐动态 for t 1:T Constraints [Constraints, S_co2(t1) S_co2(t) E_cap(t) - C_co2(t)]; end Constraints [Constraints, 0 S_co2 S_max]; Constraints [Constraints, S_co2(1) S_co2(T1)];这套骨架我每次搭新模型都从它开始改。YALMIP这层封装比较友好约束维度就是单纯的1×T向量比较不容易写错。唯一的提醒是H_ccs需要先用系数算出来再放进约束不然等式两边都有决策变量展开时会乱。我习惯先定义所有辅助变量再组装约束最后写目标函数条理最清晰。3.3 求解器配置、可行性检查与结果导出求解器配置看似小事其实踩坑不少。我的标准模板是ops sdpsettings(solver, cplex, verbose, 2); sol optimize(Constraints, cost, ops); if sol.problem 0 fprintf(求解成功\n); P_chp_opt value(P_chp); H_chp_opt value(H_chp); else disp(sol.info); pause; end这里有两个重点。第一多目标别直接加警示和权重后就让求解器闷头跑最好先单独跑一次不考虑弃风惩罚的模型看目标函数量级再设定惩罚系数至少比正常量级高一个数量级。否则会出现一种很尴尬的情况弃风惩罚设定太低优化器觉得弃点风比开P2G更划算结果明明能消纳却故意弃风。第二value()函数在YALMIP里用来提取解勘察解时我习惯一次性把所有变量的value结果塞进一个小结构体后续画图和分析都方便。可行性检查也是我必须做的一步。如果sol.problem不为0先不要怀疑求解器坏了先检查约束是否写成了不一致的等式比如电平衡里忘了加负荷或者储碳罐初值和末值约束跟实际数据冲突。我通常用sdpsettings(debug,1)打开YALMIP的调试模式它能定位到具体哪条约束不可行省去大量盲试时间。4. 算例验证三套场景跑完碳循环到底省在哪4.1 算例数据与场景设计复现研究论文时算例数据是最大的坎。我采用了一套典型日数据风电数据取自某风电场冬季典型日功率曲线电负荷和热负荷也采用典型冬季日曲线。CHP参数参考了一台100 MW级抽凝机组的公开数据P2G总效率取50%CCS捕集效率上限90%。我设计了三个场景场景A传统CHP调度不含CCS和P2G风电可弃场景B加CCS但P2G不投运捕集到的CO2外售或封存场景C加CCS加P2G捕集的CO2部分供P2G合成甲烷这样设计的好处是结果能一层层解释B相对A体现CCS的减排代价C相对B体现P2G对碳循环和弃风消纳的贡献。4.2 运行成本、碳排放与弃风率的横向对比算例结果以A场景各项指标为基准1.00做归一化参数为典型文献规模数据指标场景A传统CHP场景BCHPCCS场景CCHPCCSP2G总运行成本1.001.100.96碳排放量1.000.430.31弃风率12%18%2%这个结果我最初看到时还愣了一下加了CCS后弃风率反而从12%涨到18%仔细一想才明白CCS本身吃电挤占了风电的上网空间在没有柔性负荷帮忙消纳的情况下弃风自然更严重。而场景C引入P2G后P2G变成一个巨大的灵活电负荷弃风率被压到2%同时捕集的CO2被拿去产甲烷也算创造收入总的运行成本甚至比传统基准场景还低了一点。碳排放方面场景B直接把排放压到基准的43%场景C更进一步到31%。注意这31%是净排放即总排放减去捕集并封存的碳被P2G利用的那部分碳虽然最终燃烧后又回到大气但已经算过一次燃料的能量利用系统层面的净碳排账就是这样算的。4.3 结果背后的机理CCS调度的双刃剑效应这个算例最值得琢磨的是CCS的双重身份。一方面它是减排功臣另一方面它是电耗大户。场景B里CCS的运行完全受制于机组出力风电大发时CHP压低出力CCS跟着没烟气可捕捕集量也上不去风电小发时段CHP顶上来CCS又消耗大量电能抬高用电高峰。结果就是CCS和风电错位弃风不降反升。场景C能解决这个错位靠的是储碳罐和P2G的配合。风电大发时P2G加大吃电产气需要CO2这时候储碳罐里之前存的CO2就派上用场了等到风电乏力、CHP顶上去时CCS大量捕碳一边补充储碳罐一边继续供P2G或其他用途。换句话说P2G不再需要和CCS同频运行中间加了一个缓冲罐整个系统的时序耦合就被打通了。这就是我把储碳罐动态约束单独列为一个小节的原因。很多简化模型把P2G和CCS直连默认捕多少用多少完全丢失了这个时间解耦能力得到的结论自然比完整模型差一大截。5. 复现过程中最容易被卡住的地方我的排查思路5.1 非线性约束的线性化处理Big-M与分段线性这个模型里最危险的非线性源有三个CHP燃料成本曲线里的二次项、P2G效率随负荷变化的分段特性、还有设备启停状态和出力的乘积项。我处理的原则是能线性化就不保留非线性。燃料成本曲线如果给的是二次函数就按工作点做分段线性化把二次项拆成几个线性段用YALMIP的线性分段功能或者自己写Big-M约束。P2G效率我直接取常数因为目前文献报道的总效率差异在45%~60%之间对日前调度的结论影响远小于碳价和负荷预测误差。启停状态最容易把人坑进去。如果CHP可以启停就需要二进制变量u_chp约束里会出现u_chp·P_chp这类乘积。你说它是线性的吧变量的变量相乘它不是凸的CPLEX会直接报错或者退化成极端慢的混合整数二次规划。正确做法是用大M法把机组停机时出力为0写成P_chp ≤ M·u_chp用M表示一个足够大的数。这个M不要取太大否则数值稳定性变差一般取机组最大出力的1.2倍左右就够。5.2 YALMIL/CPLEX版本与环境不一致的坑Matlab做优化最让人崩溃的问题往往不是模型错而是环境错。我用MATLAB 2023b时遇到过YALMIP旧版本和CPLEX新版本不兼容的情况具体表现是求解器报Unknown solver或Index exceeds array bounds这种完全摸不着头脑的错。排查了两次才明白是YALMIP的求解器清单solver list里没认到CPLEX。我的解法是先把CPLEX装好然后在YALMIP命令行里跑yalmiptest看那些诊断信息确认solver后面出现cplex几个字再去model里调。如果还是认不到手动加一句solverPath设置CPLEX可执行文件路径。另外要注意2023b以后的matlab把一些老函数换了名直接影响第三方工具箱我建议直接装最新版YALMIPGitHub上那个R2023xxxx版本别用老版本硬抗。5.3 模型规模膨胀从24时段到多机组怎么控制求解时间第一个版本我写完直接跑24时段CPLEX大概几秒钟出结果当时觉得顺畅得很。后来有人问我加多台CHP和多台P2G怎么办我把变量维数翻了3倍求解时间暴涨到十几分钟。问题出在储碳罐约束里每个时段的S_co2和C_co2相互依赖加上二进制变量的组合爆炸一时很难收敛。我做的优化有两点。一是尽量用向量化约束代替for循环组装YALMIP里对于逐时段相同的约束如果约束矩阵是带Toeplitz结构可以整个块写进去速度提升明显。二是把非必要的二进制变量全部换成连续变量。比如P2G如果不需要考虑启停成本就完全用连续变量建模不做0-1决策整个问题立刻从MILP简化成LP秒收敛。还有一个通用技巧先用粗粒度跑一遍比如把24小时聚合成6个典型时段确认模型逻辑没有bug再切成24时段跑精修的。我犯过的错是直接用24时段跑结果储碳罐初始值约束写反跑了十分钟才发现结果全错。5.4 数据来源与参数标定别随便从论文里抄数字做这类复现最花时间的其实是数据。我总结了三类来源。第一类是机组公开技术手册CHP的可行域参数、效率曲线可以从厂家的典型运行数据里找。第二类是电网和气象公开数据风电出力曲线可以从一些开源数据集里拿。第三类是文献里的模型参数比如CCS的再生热耗、P2G总效率可以从综述论文里获取。但这里要特别提醒文献里给的参数往往是在特定假设下的直接抄会跟自己的系统对不上。我的做法是把关键参数都设成一个参数结构体放在文件最前面用注释标清楚来源后续调参只需要改这一个文件。比如碳价现在按100元/t算如果政策变了或交易所价格波动改一个数字整个系统重新跑一遍就行不需要动模型。另外数据的量纲一致性绝对不能出问题。我最狼狈的一次是把CO2捕集量的单位写成了kg而P2G耗碳量单位是t结果目标函数里碳成本凭空大了1000倍求解器倒没报错只是跑出来的调度策略傻得不忍直视。后来我所有的单位换算都集中在一个unit_convert.m文件里再也不用到处填系数。复现这种综合能源系统的优化模型说到底不是把设备和约束堆起来那么简单而是要在模型里保留系统的真实灵活性。我最大的体会是P2G和CCS这个组合真正的价值不是多装了两台设备而是给调度系统增加了两个自由度——一个可调的碳捕集率一个可调的柔性电负荷再加上储碳罐的时间解耦CHP的以热定电约束才算真正被打破。建模时如果把这三个自由度丢了那代码写得再漂亮算出来的也只是一个新瓶子装旧酒。