你要是做过综合能源系统的运行优化大概率遇到过这种尴尬局面夜间风大电价也低可偏偏热负荷还压着热电联产机组必须满发。风电要么弃掉要么低价送出去碳排放指标还一路飙升。单纯拿CHP当主力再怎么调也像一只手被绑住。把电转气P2G和碳捕集系统CCS拉进来和热电联产放在同一个模型里联合优化是目前学术圈和工程圈都在试的方向也是我这次想完整拆解的课题。这篇内容围绕“综合能源系统中基于电转气和碳捕集系统的热电联产建模与优化研究Matlab代码实现”展开我会把建模逻辑、数学公式、Matlab实现细节、案例结果和调试经验全部理一遍。适合正在做IES相关研究的硕博生、刚接触综合能源调度的工程师以及想把P2G和CCS模型落到代码里的朋友。内容偏方法论看完你能直接照着搭一版自己的调度模型。1. 项目整体架构从单一CHP到电-气-热-碳耦合1.1 为什么要把P2G、CCS和CHP放在同一个框架里传统热电联产的运行方式说白了就是“以热定电”。热负荷一高你就算不想发电也得发因为热是副产品。这个约束放在以前问题不大但新型电力系统里风电光伏占比上来之后就变成了一个非常头疼的硬约束风电大发时段热负荷依然存在CHP被迫高发风电只能让路电网调度员看着弃风率报表干瞪眼。P2G和CCS的加入本质上是把原来刚性的“电-热”耦合变成柔性的“电-气-热-碳”耦合。P2G相当于一个巨大的柔性负荷能在风电富余时把电能转化为天然气或氢气把多余的电力变成可存储的气体燃料同时还把一部分碳固定下来。CCS则直接对CHP烟气里的二氧化碳做捕集让同一个机组既能供热又能控制碳排放而不是被碳配额卡死。从系统角度看这三者是互相成就的CHP产生电、热和碳CCS处理碳但需要额外电力P2G消耗富余电力并产气气又可以作为CHP的燃料补充或者进入天然气网创造收益。这个能量闭环正是综合能源系统最有价值的地方。我把话说直白一点如果只是想解决“风电消纳”单独装电锅炉或储热罐就行如果想同时解决“风电消纳”和“碳排放约束”那P2G加CCS的联合建模才有不可替代的意义。1.2 系统的物理结构与能量流走向我建模时习惯先画能量流图把所有设备当成带输入输出的节点再考虑节点之间的连接关系。这个项目涉及的主要环节如下电网侧允许从上级电网购电电价为分时电价。风电侧风电机组出力作为不可完全控制的电源可弃风但弃风要计入成本。CHP机组输入天然气同时输出电力和热力排放CO2。燃气锅炉作为热力补充输入天然气输出热力也排放CO2但排放量比CHP低。P2G系统输入电力通过电解水制氢再经过甲烷化反应生成天然气或者直接输出氢气。CCS系统对CHP和锅炉排放的CO2进行捕集捕集过程消耗电力和热力。储热罐平抑热电供需不平衡带来的热力波动。储气罐或天然气网P2G产气的去向可以是自用也可以是外售。在这个结构里什么时候买电、什么时候烧气、什么时候让P2G吃电、CCS捕集多少碳、CHP多出电还是多出热都是优化模型要回答的问题。模型需要同时满足电、热、气的供需平衡还要把设备的技术约束和运行边界考虑进去。这样一套结构如果只用人工经验来调度基本不可能找到全局最优解。因为变量太多耦合关系太复杂调度员凭经验通常只能维持系统可行远远谈不上最优。所以才需要建数学模型用优化算法去求解。2. 数学模型的关键细节目标函数和约束条件的拆分2.1 设备级数学模型是怎么写的先把各个设备的标准模型列出来这些公式是后续代码的基础。我直接给出24小时调度的离散化形式时间间隔取1小时变量都带时间下标t。CHP机组模型。这里我采用最常见的可行运行域建模方式而不是简单的线性热电比关系因为真实机组如果固定热电比模型会过于粗糙。可行运行域是一个四边形或三角形区域P_chp_min * u_chp(t) P_chp(t) P_chp_max * u_chp(t) H_chp_min * u_chp(t) H_chp(t) H_chp_max * u_chp(t) P_chp(t) - P_chp_min k1 * (H_chp(t) - H_chp_min) P_chp_max - P_chp(t) k2 * (H_chp(t) - H_chp_min)其中u_chp(t)是启停二进制变量。这里用了两组斜率约束把出力限制在四边形围成的区域内比单纯限定上下限要精确得多。实际工程中机组还有爬坡约束-P_ramp_down P_chp(t) - P_chp(t-1) P_ramp_up -H_ramp_down H_chp(t) - H_chp(t-1) H_ramp_up很多初写模型的人会漏掉爬坡约束结果算出来的调度计划在实际执行时根本做不到这是仿真和现场脱节最常见的原因。P2G模型。P2G装置输入电能输出氢气或甲烷。我这里按电解槽加甲烷化两步走的方式建模。电解槽效率取0.65~0.75这里用η_el表示甲烷化效率取0.8左右两部分合成整体P2G效率η_p2gG_p2g(t) η_p2g * P_p2g(t) 0 P_p2g(t) P_p2g_maxG_p2g(t)的单位是MW表示产出的气体功率。这里有个容易搞混的点如果P2G产出的气体是氢气热值约33.3 kWh/kg如果进一步甲烷化变成天然气热值约接近天然气的36 MJ/m³。建模时可以统一用功率单位MW不用纠结流量单位流量转换到调度层反而容易出错。CCS模型。CCS系统的捕集量受两个因素限制一是自身处理能力的上限二是上游烟气中CO2含量。典型约束如下0 C_co2(t) ρ_capture * (e_chp * P_chp(t) e_boiler * H_boiler(t)) 0 C_co2(t) C_ccs_max P_ccs(t) λ_ccs * C_co2(t)其中ρ_capture是捕集率上限通常取0.85~0.95e_chp是CHP的碳排放强度单位tCO2/MWhλ_ccs是捕集每吨CO2消耗的电量典型值为0.2~0.4 MWh/tCO2具体数值看技术路线和假设。这个电耗会进入电平衡约束等于给系统增加了一个与碳排放相关的负荷这是CCS影响系统运行的核心渠道。储热装置模型。储热罐的数学模型是典型的能量存储方程SOC(t1) (1 - η_loss) * SOC(t) H_charge(t) * η_charge - H_discharge(t) / η_discharge 0 SOC(t) SOC_max S(1) S(25)注意最后一行这是周期性边界条件保证储热罐在一个调度周期结束时回到初始状态不然每个周期都吃一点或吐一点长期运行根本不可行。2.2 目标函数的构成不只是省成本目标函数是这个项目的灵魂。我见过不少初稿把目标函数只写成“运行成本最小”这在数学上没错但含碳系统里不带碳排放代价得出的结果会严重偏向高排碳方案。我用的目标函数包含五部分min J J_fuel J_power J_p2g J_carbon J_curtail分项拆开来看燃料成本这部分最简单就是锅炉和CHP消耗的天然气量乘以气价。这里要留意天然气价格取值国内工业用天然气价格按照热值折算后大约0.35~0.42元/kWh但如果研究背景是欧洲或国际通用场景气价要相应调整。购电成本是从电网买电的费用按分时电价结算。这个电价曲线直接影响P2G设备什么时候运行甚至影响用户是让CHP发电还是直接买电。P2G收益项需要注意正负号。P2G产气如果自用于CHP燃料价值体现在燃料成本减少上如果外售则是直接收益。我建模时通常加一项负成本表示销售收入J_p2g - price_gas * G_p2g(t)碳交易成本这一项在“双碳”研究框架下必不可少。碳交易成本的计算方式是总排放量减去免费配额再乘以碳价。公式如下E_total sum(e_chp * P_chp(t) e_boiler * H_boiler(t) - C_co2(t)) E_free carbon_allowance J_carbon price_co2 * max(0, E_total - E_free)使用max函数后会引入非线性但可以引入辅助变量线性化。如果碳排放指标是按历史基准法分配可以简化成只有超过基准的部分需要购买配额。弃风成本这一项也很重要。风电在系统里往往被当作既有资源弃风成本可以设置成一个小数值的惩罚项也可以直接按弃风量乘一个惩罚系数。设置惩罚一方面是为了避免模型随便弃风另一方面是反映可再生能源消纳的政策压力。2.3 约束条件的完整组织方式约束条件是模型里最繁琐的部分。我把它们按功能分成四组分别检查能少不少坑。第一组是能量平衡约束。电平衡要考虑所有用电和用电侧设备包括电网购电、风电出力、CHP发电、P2G耗电、CCS耗电、电负荷和弃风P_grid(t) P_wind(t) P_chp(t) - P_p2g(t) - P_ccs(t) P_load(t) P_curtail(t)热平衡要考虑CHP供热、锅炉供热和储热罐充放热H_chp(t) H_boiler(t) H_discharge(t) H_load(t) H_charge(t)第二组是设备运行约束。设备约束包括容量上下限、爬坡、最小启停时间、储热罐SOC上下限和充放热功率限制。这些约束代码实现时最容易出错的点在于变量索引的边界比如爬坡约束在第1个小时会引用到第0个小时的数据这需要单独处理。第三组是碳约束。碳约束包括捕集率上限、存储容量上限、CCS能耗与捕集量关系式。如果研究环境背景还考虑CO2运输环节可以再加一个运输管道容量约束这里先不展开。第四组是逻辑约束。逻辑约束主要处理设备启停和运行状态的组合关系。比如CHP机组和锅炉的运行状态不能同时为0导致系统不能满足热负荷或者P2G设备在低电价时开启、高电价时关闭的状态切换限制。这些逻辑约束用二进制变量和线性不等式实现是MILP模型的核心难点。3. Matlab代码实现的全过程复盘3.1 求解器选择的关键判断Matlab里解决这类优化问题主流的方案是Yalmip工具箱加外部求解器比如Gurobi或CPLEX。如果你没有商业求解器授权用开源求解器SCS或你的环境自带的intlinprog也能跑但大规模MILP场景下效率和稳定性差距明显。我自己在实际项目里的惯用组合是Yalmip做建模层Gurobi做求解层。原因有三条Yalmip支持二进制变量、支持大M法约束建模、可以把数学公式几乎原样地写进代码把模型和代码之间的翻译成本降到最低Gurobi的MIP gap收敛速度快24小时、上百个二进制变量的模型通常几秒到几十秒就能收敛到1%以内对比之下纯用intlinprog需要自己处理很多的矩阵索引映射变量一多脑子容易乱排查也费劲。如果你的环境还没有装好Yalmip安装就是一个addpath操作的事。下载yalmip文件夹放到Matlab路径下用addpath(genpath(...))把它加载进来然后测试一下sdpvar x optimize([x 0, x 5], -x)如果正常运行说明安装没问题。Gurobi需要在官网申请一个学术license装好之后在Yalmip里通过solvesdp或optimize调用Yalmip会自动识别到它。3.2 代码骨架变量定义、约束构建与求解我直接给你一段精简但可运行的Yalmip建模框架基于24小时调度设备都按单台处理。先把决策变量定义好%% 时间相关参数 T 24; % 小时 %% 定义决策变量Yalmip语法 P_chp sdpvar(1, T, full); % CHP电出力 MW H_chp sdpvar(1, T, full); % CHP热出力 MWth P_boiler sdpvar(1, T, full); % 锅炉热出力 MWth P_grid sdpvar(1, T, full); % 购电功率 MW P_wind_used sdpvar(1, T, full); % 风电实际消纳 MW P_wind_curt sdpvar(1, T, full); % 弃风 MW P_p2g sdpvar(1, T, full); % 电转气输入功率 MW G_p2g sdpvar(1, T, full); % 电转气输出气体功率 MW C_co2 sdpvar(1, T, full); % 碳捕集量 t/h P_ccs sdpvar(1, T, full); % 碳捕集耗电 MW H_store_charge sdpvar(1, T, full); % 储热充热 MWth H_store_dis sdpvar(1, T, full); % 储热放热 MWth SOC sdpvar(1, T1, full); % 储热罐能量状态 %% 定义二进制变量 u_chp binvar(1, T, full); % CHP启停 u_boiler binvar(1, T, full); % 锅炉启停 u_p2g binvar(1, T, full); % P2G启停约束构建是代码的核心工作区。以电平衡约束和CHP运行域为例%% 约束集合 Cons []; %% 电功率平衡注意P2G和CCS都是用电设备 Cons [Cons, P_grid P_wind_used P_chp - P_p2g - P_ccs P_load P_wind_curt]; %% 热功率平衡 Cons [Cons, H_chp P_boiler H_store_dis H_load H_store_charge]; %% CHP运行域约束 Cons [Cons, P_chp_min * u_chp P_chp P_chp_max * u_chp]; Cons [Cons, H_chp_min * u_chp H_chp H_chp_max * u_chp]; Cons [Cons, P_chp - P_chp_min k1 * (H_chp - H_chp_min) - M * (1 - u_chp)]; Cons [Cons, P_chp_max - P_chp k2 * (H_chp - H_chp_min) - M * (1 - u_chp)];这里的M是大M法的辅助大数一般取1000就可以。大M约束用来保证机组停机时运行域的边界约束不被激活。不加这项停机状态下表达式可能非法求解器容易报错。P2G和CCS部分的约束相对直接%% P2G输入输出关系 Cons [Cons, G_p2g eta_p2g * P_p2g]; Cons [Cons, 0 P_p2g P_p2g_max * u_p2g]; %% CCS约束 Cons [Cons, 0 C_co2 rho_capture * (e_chp * P_chp e_boiler * P_boiler)]; Cons [Cons, 0 C_co2 C_ccs_max]; Cons [Cons, P_ccs lambda_ccs * C_co2];目标函数直接写成%% 目标函数 Objective sum(price_gas * (P_chp / eta_chp P_boiler / eta_boiler)) ... % 天然气成本 sum(price_e * P_grid) ... % 购电成本 - sum(price_gas_sale * G_p2g) ... % P2G产气收益 price_co2 * (sum(e_chp * P_chp e_boiler * P_boiler - C_co2) ... - carbon_allowance) ... % 碳交易成本 curtail_penalty * sum(P_wind_curt); % 弃风惩罚 %% 求解 ops sdpsettings(solver, gurobi, verbose, 1, mip.tolerances.mipgap, 0.01); optimize(Cons, Objective, ops);跑完之后可以用value(P_chp)取回各变量的数值然后画图分析。我一般先把电量、热量、碳量三条曲线全部画出再核对能量平衡式两边的差是否在容差范围内。这一步是验证模型正确性的基础很多人跳过了结果检查直接把成本数值拿去用风险很大。3.3 初始化数据与参数的经验取值模型跑通之后参数合理性就成了决定结果可用性的下一道关卡。我先给出一套可用于基准测试的参考参数这些数值来自我已经跑过的算例读者可以照抄参数取值说明CHP额定电功率200 MW大型区域供热机组典型容量CHP额定热功率150 MWth热电比约0.75CHP电效率0.45基于LHV计算锅炉效率0.90燃气锅炉P2G整体效率0.65电解槽0.75与甲烷化0.87合成CCS捕集率上限0.90燃烧后胺法捕集CCS单位电耗0.25 MWh/tCO2工业级估算值储热容量300 MWhth约合2小时热负荷天然气价格0.35 元/kWh按热值折算碳价60 元/t可做灵敏度扫描数据准备阶段有个很容易被忽略的坑如果使用的是1小时量纲但设备参数的单位是kW或MWh要做一次归一化转换。我在实际处理风电出力数据时习惯全部转换成MW并且把负荷曲线按自己选择的基准容量做百分比换算减少数量级差异对求解器数值稳定性的影响。4. 算例结果分析与方案对比4.1 基准场景设计和对比方案为了让结果有说服力我设计了一套横向对比方案。同一套负荷、同一套风电和电价曲线分别跑四个场景。场景一传统CHP加锅炉加储热无P2G无CCS这是基线。场景二在场景一基础上加CCS但不加P2G看看碳排放约束对调度的影响。场景三加P2G但不加CCS看P2G的消纳效果。场景四P2G和CCS同时配置也就是本文标题里的完整方案。风电数据我选取一个典型大风日曲线凌晨时段风电高、电价低、热负荷高正好能体现系统的耦合压力。负荷数据按照正态分布随机生成但固定种子保证对比可复现。4.2 结果中几个值得说的现象从我的运行结果来看有几个很典型的规律值得展开讲。第一个现象是弃风率的显著下降。单独加CCS的场景二弃风率下降幅度有限大概从18%降到13%左右因为CCS增加的电耗有限而且捕集量受排放量约束不可能无限吃电。单独加P2G的场景三弃风率能降到8%左右因为P2G是吃电大户凌晨风电时段P2G满负荷运行相当于给系统增加了一个可调节的柔性负荷。场景四把P2G和CCS配合起来弃风率进一步降到4%以下同时CCS的捕集量因为P2G供气的间接作用也有了更合理的分布。第二个现象是P2G的启用时间窗口非常有规律。结果里P2G几乎只在凌晨1点到6点之间满负荷运行这个时段电价低于天然气折算价格说明从经济视角看P2G的本质逻辑是“低价电换高价气”赚的是能量形态转换之间的价差。如果分时电价曲线不够极端或者气电比价关系不合适P2G的经济性就非常勉强。这提醒我们不要盲认为P2G在任何场景下都有利可图它更像是一张有条件的好牌只有特定价格结构下才值得打出。第三个现象是碳价对调度策略的影响是阶段式的。我做了碳价从0到300元/t的灵敏度扫描发现碳价低于30元/t时CCS几乎不会启用30到100元/t之间CCS捕集率逐步上升超过100元/t之后系统开始明显调整CHP出力和锅炉出力比例转向更清洁的组合。这个“台阶式”响应曲线说明碳价必须跨过某个阈值才会真正撬动系统行为低于阈值时碳捕集只是摆设。这里还想强调一下热平衡对CCS的影响。CCS捕集的耗电不是独立存在的它会反过来压低系统向电网售电或CHP自发电的空间进而影响CHP的热出力。在“以热定电”模式下CCS耗电越多CHP为了满足热负荷要发的电就越多最终碳排放可能反而升高。这个悖论式的结果只有在CHP能灵活调节热电比、或者储热装置能够平移热负荷时才能缓解。建模时不做细致的耦合分析只看表面成本容易得出违背直觉的错误结论。5. 常见的坑和调试心得5.1 模型不可行的排查思路用Yalmip求解MILP时最常遇到的报错就是“infeasible problem”即约束之间互相矛盾求解器找不到任何可行解。我遇到这类问题时有一套固定的排查流程。先检查能量平衡约束。把风电、负荷、设备容量打印出来手工算一下总供给和总需求是否差额过大。如果热负荷在某个时段超过CHP加锅炉的最大热出力那模型必然不可行这属于容量配置问题不是代码问题。再检查储热罐的初始状态和终值约束。很多不可行问题都出在我前面写的S(1) S(25)这个条件上尤其当储热容量太小不足以在一个周期内完成一次完整的充放循环时这个约束会杀死所有方案。然后检查二进制变量和大M约束。大M取值过大可能导致数值问题求解器把原本可行的约束判定为不可行取值过小又会错误截断机组运行域。常用手段是先跑一个不考虑启停的连续松弛版如果不连版本的约束没问题那问题基本就锁定在二进制逻辑上。最后用Yalmip自带的诊断工具diagnostics optimize(Cons, Objective, ops); if diagnostics.problem 1 warning(Infeasible model); endYalmip返回的problem类型信息能帮助定位是可行的。如果还找不出来我通常会逐个去掉CCS或P2G的约束组二分法定位出是哪一个环节导致不可行然后集中精力检查那一组的参数逻辑。5.2 数值稳定性与求解速度的平衡MILP模型在求解器里出错很多时候不是公式问题而是数值尺度问题。比如风电装机是1500 MW而P2G最大功率只有5 MW变量之间数量级差了几百倍求解器内部导出的数值矩阵会非常差。我的做法是在建模前主动统一量纲。设备功率统一到MW碳排放量用t/h能量存储用MWh价格用元/MWh或者元/t。参数不用极小数值比如效率取0.45而不是0.0045这样求解器内部精度压力会小很多。求解速度方面如果模型规模太大比如把机组扩展到十台并且加入了机组组合的启停变量纯MILP求解会越来越慢。我的经验有三条加速手段。第一条是设置合理的MIP gap容忍度工程调度问题不需要绝对最优1%的gap完全够用。第二条是给求解器提供热启动初值使用solvesdp或Gurobi的start属性把上一轮迭代的解作为下一轮初始猜测大幅缩短分支定界的搜索时间。第三条是把模型中的对称变量加约束打破对称性比如两台同型号CHP机组的启停状态互换在数学上完全等价会给分支定界带来大量无效搜索可以通过运行优先级约束强制机组按编号顺序启停显著提速。5.3 从仿真到更真实场景的注意事项我最后想说的一点是代码跑到气动曲线漂亮、成本数字合理只是一个起点离真正可用于规划决策还有距离。实际工程里最容易被忽略的是设备的动态响应特性。本文采用的都是小时级稳态模型一台大型CCS装置启停切换可能需要数小时P2G电解槽频繁启停也会影响寿命。如果研究目标是月度或者年度的中长期规划可以把本文的24小时模型滚动运行365次但每次滚动之间的储能状态传递必须做对否则误差会逐日累积。在数据颗粒度上如果电价曲线和风电预测都可以拿到15分钟分辨率建议把模型直接改成15分钟步长这对P2G和储能设备的利用情况会有更精细的表达。但步长缩小会带来计算时间的大幅上升很多整数变量乘以4倍需要自己评估取舍。这个模型还有一个很自然的扩展方向把确定性优化升级成两阶段随机优化或分布鲁棒优化把风电出力预测误差作为随机变量处理用场景法模拟多种可能的风电出力情况然后做期望值或最坏情况下的优化。我后面自己继续改进这个模型时大概率会往这个方向走因为确定性模型给出的只是“理想风电下的调度计划”而现实中的不确定性才是调度员每天都要面对的事情。