前阵子给一个园区级虚拟电厂做优化调度白天光伏一上来晚间风电又满发偏偏深夜负荷往下掉燃气轮机只能压到最低技术出力甚至停机。起初只看电功率平衡问题勉强靠弃风解决可一旦把碳成本、P2G-CCS耦合和燃气掺氢这些环节加进去调度决策就完全变了——深夜的多余风电不再只是“多余的电”它可以变成氢气、变成天然气甚至变成碳减排额度。这就是我写下这篇分享的原因。这篇内容围绕“基于阶梯碳交易的含P2G-CCS耦合和燃气掺氢的虚拟电厂优化调度”展开把我实际用Matlab跑通的一套实现思路完整讲清楚从阶梯碳交易怎么建模到P2G-CCS能量流和物质流怎么耦合再到燃气掺氢的热值修正、优化模型的目标函数与约束最后是YalmipGurobi环境下的代码细节和排坑经验。适合正在做虚拟电厂、综合能源系统优化调度尤其是想把碳交易机制和电转气/碳捕集纳入模型的研究生和工程师参考。1. 把调度对象先拆清楚P2G-CCS、掺氢燃气轮机在虚拟电厂里的位置1.1 一次具体的调度困境低谷时段为什么愁假设一个典型日园区负荷凌晨时段在40-50MW左右风电出力却有70MW。没有储能、没有P2G时风电只能压出力或弃掉。如果燃气轮机还要带最低负荷那电功率平衡就更紧要么燃气轮机停要么风电弃。这时候如果把“碳”作为约束加进来单纯弃风反而是“安全”的选择因为燃气轮机一启动就会产生碳排放碳配额超了要付出额外成本。可问题是弃风本身就是一种经济损失。风电边际成本接近零弃掉等于白白浪费了绿色电力。而P2G设备的出现让“多余电能”多了一个出口电解水制氢氢气可以储起来也可以直接和捕集到的CO2甲烷化生成天然气供燃气轮机使用。这样一来低谷时段的弃风可以通过P2G变成“燃料”而燃气轮机掺氢燃烧又能降低单位出力的碳排放。整个虚拟电厂的调度就从单纯的“电功率平衡”变成了“电、气、氢、碳”四个网络耦合的优化问题。1.2 虚拟电厂各单元的功能与能量流这类系统中各单元的任务可以简单归类风电、光伏可再生电源出力具有不确定性优先消纳。燃气轮机可控电源承担基荷或调峰但会产生碳排放。P2G单元电解槽制氢氢气一部分进入储氢罐一部分进入甲烷化反应器。CCS单元从燃气轮机或外部烟气中捕集CO2捕集后的CO2送入甲烷化反应器作为合成天然气的碳原料。甲烷化反应器将H2和CO2合成为CH4生成合成天然气可注入天然气网络存储。储氢、储气、蓄电池提供时间维度的解耦能力。电负荷、气负荷终端用户需求。这里面最关键的耦合是P2G需要电CCS也需要电甲烷化需要CO2和H2燃气轮机可以使用天然气和氢的混合物碳排放总量影响碳交易成本。每一条能量流背后都有一个成本项或约束项建模时最怕的是只考虑功率平衡而忽略物质流平衡。1.3 耦合关系引发的建模难点耦合关系带来的直接难点是变量维度膨胀。传统机组组合问题里有连续变量、0-1变量加了P2G-CCS后还需要同时跟踪氢气流量、CO2流量、天然气流量。如果每个时段都定义一套变量模型的规模会翻倍。另一个难点是时间尺度。电解槽可以在分钟级调整出力甲烷化反应器则有一定惰性燃气轮机掺氢比例受限于设备改造上限。这些限制都需要用约束表达否则优化结果会“天马行空”比如让掺氢比例在1%和40%之间瞬间跳变。所以我在建模时特意加入了调节速率约束和比例上下限约束宁可保守一些也不能让结果在工程上不可执行。2. 阶梯碳交易建模从统一碳价到分段惩罚差别在“边际压力”2.1 统一碳价为什么不够用很多入门模型假设碳价是固定常数例如每吨CO2 60元这样碳成本就是排放量与单价的简单乘积。这种模型的好处是线性、好求解但缺陷是无论排放量多高边际碳价始终不变系统只要排放量小于某个阈值就没有进一步减排的动力。实际碳市场中配额免费发放超出部分按阶梯价格购买超得越多碳价越贵减排压力逐级增加。阶梯碳交易的核心思想就是对“配额内排放”和“超配额排放”区别对待。配额内排放不产生购买成本甚至可以通过出售富余配额获得收益超配额部分划分若干区间每一区间的价格高于上一区间。这样优化的边际成本是分段上升的调度模型会自动想办法把排放量压到更低的阶梯内。我在项目中采用的是“购买阶梯惩罚”模型不涉及出售配额这样偏保守也更容易收敛。如果要把出售配额也考虑进去可以再加一个对称的收益分段函数但要注意买卖价格不等时模型会套利需要增加约束防止同时买卖。2.2 阶梯碳交易的分段函数表达设调度周期内有T个时段免费碳配额总量为E0可以是全周期总量的配额也可以分摊到每个时段。为了简化我按整个调度周期核算总碳排放量[ C_{CO2} \begin{cases} 0 E \le E_0 \ p_1 (E-E_0) E_0 E \le E_0\Delta_1 \ p_1 \Delta_1 p_2 (E-E_0-\Delta_1) E_0\Delta_1 E \le E_0\Delta_1\Delta_2 \ ... \end{cases} ]其中 (E) 为总碳排放量(\Delta_1, \Delta_2) 为阶梯区间长度(p_1, p_2,...) 为对应阶梯碳价。实际文献中阶梯碳价常取 (p_21.2p_1)、(p_31.4p_1) 等递增系数。关键点在于这个函数是分段线性的直接写进目标函数会引入非线性需要线性化。我在Matlab里采用的是“区间变量拆分”方式将超出配额部分拆成若干区间变量 (E_{seg,k})满足[ 0 \le E_{seg,1} \le \Delta_1,\quad 0 \le E_{seg,2} \le \Delta_2,\quad ... ] [ E - E_0 \sum_k E_{seg,k} ]目标函数中加入 (\sum_k p_k E_{seg,k})。为了保证低阶梯先被使用除了约束区间上限外还可以依赖p递增的特性优化器自然会优先填满低价区间。这个线性化方式比big-M引入二进制变量更简单求解速度也更快。2.3 配额分配与参数设置的细节点配额计算有两种常见口径按历史排放强度分配祖父法和按设定基准线分配基准法。我在仿真中采用的是“按出力-排放基准线”方式即根据燃气轮机的有功出力和设计排放因子估算一个免费额度再乘以一个系数 (\alpha)比如0.9表示碳配额逐渐收紧。这样做还有个好处配额跟出力水平挂钩不会因为燃气轮机开得少而配额过低导致模型过度惩罚机组启停。细节点是阶梯区间长度的单位是吨调度周期若取24小时配额也应是当日配额不要混用年度配额口径。2.4 分段函数线性化big-M还是sos2除了区间拆分还可以用文献中常见的0-1变量方法引入 (S_k) 指示是否进入第k阶梯然后用big-M约束。但这种方法在每个时段都要加多个二元变量对大型机组组合问题负担很大。我个人建议优先试试区间拆分法因为阶梯价格天然递增模型会自动按顺序用。如果你的碳价并非严格递增例如某些情况下阶梯价格先高后低那就只能老老实实用二进制变量。另一个替代方案是Yalmip自带的piecewise函数。虽然用起来方便但它在背后会引入辅助变量有时候会带来额外的数值问题。我更习惯自己搭结构至少出了问题容易排查。3. P2G-CCS耦合单元碳不再是终点而是下一个反应物的原料3.1 P2G的两步反应与能量效率P2G通常分两步。第一步是电解水制氢[ 2H_2O \rightarrow 2H_2 O_2 ]第二步是甲烷化Sabatier反应[ CO_2 4H_2 \rightarrow CH_4 2H_2O ]从反应式可以看出来制备1个单位的CH4按物质的量需要4个单位的H2和1个单位的CO2。但在能量层面电解槽消耗电能效率一般取65%-75%按氢气低热值折算甲烷化反应本身放热但工程上更关心的是H2到CH4的转化效率通常按80%-90%考虑。整体P2G能量效率在50%-60%左右。建模时我会定义三个关键变量电解槽输入电功率 (P_{P2G,t})、氢气产量 (H_{H2,t})、甲烷产量 (G_{CH4,t})。用转换系数关联[ H_{H2,t} \eta_{ele} P_{P2G,t} / \lambda_H ]其中 (\lambda_H) 是单位氢量对应的热值。如果做纯“功率”建模直接把电、氢、气都换算成功率的p.u.值数学上更整齐但会在物质流上出问题。我的做法是内部统一用kW和kWh氢量和气量用热值kWh/m³或kWh/kg换算成能量单位这样碳流计算时再用整数化学反应计量比折算CO2消耗量。3.2 CCS捕集量与甲烷化需要的CO2的平衡关系CCS环节捕集CO2需要额外电能。设捕集每吨CO2的耗电为 (e_{ccs}) kWh那么CCS运行时除了捕集设备的固定电耗外还要在电功率平衡中扣除一部分“自身消耗”。耦合约束的核心是甲烷化反应消耗的CO2不超过CCS捕集量与外部供给CO2量之和。由于碳捕集成本通常高于外购CO2优化器会倾向捕集烟气中的CO2相当于燃气轮机排放的碳又被回收了一部分形成碳循环。但需要注意的是捕集量不能超过燃气轮机实际排放量即[ E_{ccs,t} \le E_{gt,t} ]否则模型会无中生有地“产生”减排量。这个约束我一开始漏掉了结果碳捕集量超过了机组排放碳交易成本变成负的出现明显的不符合物理的套利。这也是很多初稿里常见的错误。3.3 储氢储气环节为什么必须加电解槽产氢和甲烷化用氢不一定同步所以必须有储氢罐缓冲。存储环节建模与蓄电池类似[ S_{H2,t1} S_{H2,t} H_{pro,t} - H_{use,t} - H_{sale,t} ]同时要设置储氢容量上下限和充放速率。储气也类似用于存储甲烷化生成的合成天然气以及从天然气网购得的天然气。加了存储后P2G可以在夜间风电富余时低电价电解水把氢气或合成天然气存储起来白天燃气轮机掺氢燃烧或直接供气实现“削峰填谷”。存储模型最容易被忽略的是自损耗。氢分子小储存会有泄漏但24小时尺度内泄漏比例很小可以忽略如果做长时间尺度建议加一个0.1%-0.5%/h的自损耗系数。3.4 实际建模中的能量自耗和响应速度P2G和CCS本身也要消耗能量。电解槽除了消耗直流电还有辅助系统电耗捕获系统在再生塔需要热能。如果系统里有热电联产CHP捕集系统的热耗可以由余热提供如果忽略这部分热平衡会高估捕集系统的净减排收益。响应速度方面电解槽功率调节速率一般能达到每分钟10%-20%额定功率甲烷化装置调节速率稍慢。在日尺度调度中如果不考虑机组组合可以只加最大爬坡约束如果涉及实时调度建议用一阶惯性环节描述。我在参数设置时P2G最大输入功率设为10MWCCS捕集速率上限设为对应燃气轮机额定排放的80%储氢容量设置为3MW按热值折算这样一天内的灵活性基本够用。4. 燃气掺氢燃烧建模别只用体积比例热值会骗人4.1 混合燃料热值计算与标幺问题燃气轮机掺氢是当前降低碳排放的重要技术路线。掺氢比例 (\alpha) 通常指体积分数。氢气低热值LHV约10.8 MJ/Nm³天然气LHV约35.8 MJ/Nm³二者差异巨大。混合燃料的低热值按下式计算[ LHV_{mix} (1-\alpha)LHV_{CH4} \alpha LHV_{H2} ]注意这里的体积是同一标准状态下的体积。如果按质量比例算热值关系又会不同。所以写代码时一定要统一单位。我见过不少文献直接用“掺氢比例20%”然后假设燃料成本下降20%这是错的。20%体积掺氢后混合燃料热值变为[ LHV_{mix} 0.8 \times 35.8 0.2 \times 10.8 30.8 \text{ MJ/Nm³} ]相比纯天然气低热值下降约14%这意味着产生同样的热量需要的混合燃料体积要增加约16%。4.2 掺氢后的机组燃料费用和出力效率设定燃气轮机在纯天然气模式下的燃料耗量函数为[ F_{CH4,t} a P_{gt,t} b u_{gt,t} ]其中 (a)、(b) 是耗量系数(u) 是启停状态0-1变量。掺氢后由于燃料热值变化燃料消耗体积会变化但机组所需输入的热量基本不变。因此可以先将耗量函数换算成输入热量[ Q_{in,t} F_{CH4,t} \cdot LHV_{CH4} ]然后混合燃料体积消耗[ F_{mix,t} Q_{in,t} / LHV_{mix,t} ]燃气轮机效率本身也会随掺氢比例变化这与燃烧室设计有关。保守做法假设电效率不变只修正燃料热值更精细的做法加一个效率修正系数 (\eta_{gt}(\alpha))例如文献中掺氢30%时效率下降约0.5-1个百分点。我采用的是线性修正[ \eta_{gt,t} \eta_{gt,0} \cdot (1 - \beta \alpha_t) ]其中 (\beta) 取0.02左右表示掺氢带来的效率惩罚。这个参数直接影响结果如果忽略效率惩罚模型会无限堆高掺氢比例加上后结果才会出现一个合理的拐点。4.3 掺氢后的碳排放核算碳排放核算的基础是实际燃用的碳摩尔数。混合燃料中只有天然气含碳。单位体积混合燃料燃烧产生的CO2质量可以表示为[ EF_{mix} (1-\alpha) \cdot EF_{CH4} ]其中 (EF_{CH4}) 是每标准立方米纯天然气燃烧产生的CO2量。这里要特别提醒不能用“单位热值排放因子”乘以混合燃料的热值来算因为氢气的热值不为零但没有碳排放这样会高估碳排放。正确思路是先算出混合燃料的体积消耗再乘以 ((1-\alpha))从而得到碳排量。4.4 掺氢比例上限安全与设备约束掺氢比例不能无限提高。现有燃气轮机的燃烧室、燃料阀、密封和控制系统通常只能承受一定的氢份额某些机型改造后可以到50%甚至更高但日常运行中为了安全一般限制在20%-30%体积比例内。我的模型中加入[ 0 \le \alpha_t \le \alpha_{\max} ]并且加了掺氢比例变化速率约束防止每个时段跳变。还要注意如果虚拟电厂里有多个燃气轮机每一台的掺氢上限可能不同不能统一用一个比例否则问题约束不够准确。5. 完整优化模型目标函数、约束与线性化处理5.1 目标函数四类成本怎么放进一个式子里整个优化调度的目标函数可以写为最小化调度周期内的总成本。我分四块购气和购电成本从天然气网购气从主网购电。设备运行维护成本风电、光伏、燃气轮机、P2G、CCS等按出力线性计算维护费用。碳交易成本阶梯碳价产生的费用。弃风弃光惩罚成本单位弃电量设一个惩罚因子让优化器尽量消纳可再生能源。目标函数写出来就是[ \min \sum_{t1}^{T} [C_{gas,t} C_{buy,t} C_{om,t} C_{CO2,t} C_{curt,t}] ]这里面最难处理的是碳交易成本 (C_{CO2,t})因为它是全周期碳排放量的函数不是每个时段独立。所以我把碳排放总量写成所有时段的排放之和然后用2.2节的区间拆分方式建模。5.2 功率平衡、气量平衡、碳量平衡电功率平衡[ P_{wt,t} P_{pv,t} P_{gt,t} P_{buy,t} P_{dis,t} P_{load,t} P_{p2g,t} P_{ccs,t} P_{ch,t} ]其中 (P_{ch,t}) 是蓄电池充电功率(P_{dis,t}) 是放电功率。P2G和CCS的用电都作为负荷参与平衡。气量平衡[ G_{buy,t} G_{ch4,t}^{p2g} G_{out,t}^{storage} G_{load,t} G_{gt,t}^{fuel} G_{in,t}^{storage} ]碳平衡并不是系统级的总量守恒而是追踪燃气轮机燃烧产生的CO2、CCS捕集量、以及最终排放到大气中的CO2[ E_{gt,t} E_{ccs,t} E_{emit,t} ]总排放配额核算使用 (E_{emit,t}) 的累加。5.3 机组运行约束与储能SOC约束燃气轮机约束包括出力上下限、爬坡速率、最小启停时间、启停变量逻辑。由于是MILP最小启停时间会增加计算负担如果只是分析性算例可以先忽略最小启停时间只保留爬坡和出力范围。蓄电池SOC约束[ SOC_{t1} SOC_t \eta_{ch}P_{ch,t} - P_{dis,t}/\eta_{dis} ][ SOC_{\min} \le SOC_t \le SOC_{\max} ]还要防止同时充放电加一个互斥约束或者通过补充0-1变量。如果模型规模大可以用简化每时段要么充电要么放电但约束不能漏。P2G、CCS、储氢、储气都需要对应的容量、上下限、爬坡约束。这里不再一一列公式但代码里必须完整。5.4 整体模型梳理与求解器选择最终模型是一个混合整数线性规划MILP问题。因为掺氢比例、碳交易分段、充放电互斥这些都需要0-1变量线性关系可以通过分段线性函数保留。求解器方面Matlab下最方便的组合是Yalmip Gurobi或Cplex。我实际用的是Gurobi 9.5和Matlab R2021bYalmip版本为2023版。为什么不用fmincon因为MILP需要分支定界fmincon只能解连续非线性。如果强行用罚函数处理0-1容易陷入局部最优而且分段碳交易函数的非光滑性会让梯度法失效。6. Matlab代码实现与排坑YalmipGurobi下的真实操作6.1 数据结构与参数表写代码前强烈建议先把参数整理成一个结构体数组。我习惯用 struct字段名可读性强比散落的变量清晰很多。% 设备参数示例 para.gt.Pmax 30; % 燃气轮机最大出力 MW para.gt.Pmin 6; % 最小出力 MW para.gt.a 2.5; % 耗量系数 MMBTU/MWh para.gt.b 1.2; % 空载耗量 MMBTU para.gt.ramp 6; % 爬坡速率 MW/h para.gt.alpha_max 0.3; % 最大掺氢体积比例 para.p2g.eta_ele 0.7; % 电解槽效率 para.p2g.eta_meth 0.85; % 甲烷化效率 para.p2g.Pmax 10; % P2G最大电功率 MW ...时间序列数据例如风电、光伏、负荷曲线用T*1列向量。每个时段为1小时T24。6.2 关键建模代码变量、约束和分段碳成本变量定义T 24; P_gt sdpvar(1,T); % 燃气轮机出力 u_gt binvar(1,T); % 启停状态 alpha sdpvar(1,T); % 掺氢比例 P_p2g sdpvar(1,T); % P2G输入电功率 E_emit sdpvar(1,T); % 实际碳排放量 E_ccs sdpvar(1,T); % CCS捕集量 S_h2 sdpvar(1,T1); % 储氢状态 ... % 碳交易阶梯变量 dE1 sdpvar(1,1); dE2 sdpvar(1,1); dE3 sdpvar(1,1);约束添加示例——P2G产气量% 电解槽产氢量按热值折算为MW H_pro para.p2g.eta_ele * P_p2g / 3.6 * 34; % 每MWh电对应氢气热值实际要按单位换算 % 甲烷化耗氢与产气 G_ch4 para.p2g.eta_meth * H_pro / 4.2; % 简化比例 % 实际使用中建议直接把所有能量单位统一为MW反应计量比用摩尔比换算这里不展开完整物理量转化但提醒一句氢气和天然气的热值单位、燃料耗量单位如果不统一很容易在甲烷化环节算出差数倍的结果。我自己就犯过错一度以为模型可以“凭空产气”。分段碳成本约束E_total sum(E_emit); dE1 dE2 dE3 E_total - E0; 0 dE1 D1; 0 dE2 D2; 0 dE3 D3; C_co2 p1*dE1 p2*dE2 p3*dE3;Yalmip中这样直接添加即可。注意如果E_total E0差值应为0所以在约束左边需要加 max分段区间变量本来就是非负的当E_total E0时等式右边负数会导致无解。处理方法是在等式前加一个“负偏差变量”表示配额富余或者直接用max表达式。更稳妥的方式是引入 (dE_neg)dE1 dE2 dE3 - dE_neg E_total - E0; dE_neg 0;配额富余时 dE_neg 0碳交易成本为0不奖励。6.3 让模型可解的三个细节第一给所有连续变量设定合理上下限。Yalmip里不设置boundsGurobi虽然能处理但预处理效果差求解速度明显下降。尤其是SOC变量一定要bind到容量上下限。第二储能互斥约束不要用乘积。正确的线性化方式是P_ch M * u_ch; P_dis M * (1 - u_ch); u_ch binvar(1,T);M取一个足够大的数比如500但别太大否则数值不稳定。第三模型里所有的0-1约束尽量用二进制变量配合big-M而不是用ceil等其它函数。Gurobi对二元变量有很好的处理但要避免变量名中的NaN和重复定义。6.4 调试中的常见错误与解决我实际踩过几个坑坑一等式约束两边量纲不一致。比如功率是MW气体量是m³直接相加导致无解。解决办法是全部折算成MW或者标幺值。坑二碳配额E0用全年数据而调度周期是24小时导致E_total - E0永远是负数模型认为没有任何碳成本。解决方法是把配额按天数折算或按比例折算。坑三Yalmip中sdpvar和binvar的维度没对齐约束矩阵维度不一致报“Inconsistent dimensions”。解决方法是先初始化空约束再循环添加避免用矩阵批量乘法时弄错维度。坑四求解器gap设置太紧MILP长时间不收敛。实际算例中我设置相对gap为0.5%或1%完全够用。7. 算例结果给我的三个直接结论7.1 阶梯碳价倒逼燃料结构变化我用一个典型日负荷和风电数据进行测试算例规模为1台燃气轮机、10MW P2G、8MW CCS、蓄电池5MWh、储氢3MWh。在统一碳价60元/吨时累计碳排放为128吨改用电价35/55/80元/吨的三阶梯碳价后相同负荷条件下累计碳排放降到了96吨下降约25%。原因是模型把一部分天然气消耗转移到了“夜间风电电解水制氢甲烷化”环节白天燃气轮机的出力被压低高碳价的边际压力真正起了作用。这个现象说明阶梯碳价不仅增加成本更重要的是改变了不同时段机组出力的相对经济性。夜间风电便宜P2G制气成本低合成的天然气白天用整体碳排放自然降低。7.2 P2G-CCS对弃风消纳的贡献比想象中大同样场景下纯电调度弃风率为18.6%加入P2G后夜间多出来的2MW风电直接被电解槽吸收弃风率降到了10.2%再加上CCS捕集和储气环节弃风率进一步降到6.1%。P2G-CCS耦合的价值不在于“多发了几度电”而在于给系统增加了一个跨时段、跨能源品种的调节通道。要注意的是弃风率下降不等于总成本一定下降。P2G和CCS的设备折旧与维护费用不低在我的参数设定下含P2G-CCS的总运行成本比纯电调度高约4%但如果把碳交易成本算进去总成本反而略低。所以评估时要看综合成本而不是只看弃风率。7.3 掺氢比例不是越大越好我把最大掺氢比例从0逐步调到30%观察变化。掺氢从0到20%时碳排放下降明显因为氢气替代了部分天然气但超过20%后由于混合燃料热值下降燃料体积流量增加加上效率小幅下降总燃料成本上升碳排放下降速度趋缓。在阶梯碳价下最优掺氢比例基本停在20%-25%之间而不是一直冲到上限。这说明掺氢是一个“边际收益递减”的技术。实际工程中不必追求最大掺氢比例而是要根据氢气来源成本、设备效率和碳价水平做联合优化。虚拟电厂调度模型的优势就在这里它能在每个时段给出最优掺氢比例而不是人为拍一个固定值。如果你也想复现这套程序我的建议是先搭一个最简版——只有燃气轮机、风电、负荷和阶梯碳交易跑通后再逐步加入P2G、CCS、储氢和掺氢。每加一个环节就对比一下前后的功率平衡和碳排放变化这样能快速定位模型问题。我最初就是跳过这一步直接把所有单元一起建模结果模型不收敛时根本不知道是哪个环节出了问题。这个项目还有很多可以扩展的地方比如加入碳捕集率优化、考虑风电预测误差的鲁棒调度、把掺氢比例做成机组组合的一部分而非固定值。先把手里的模型跑出可信结果再往更复杂的方向走不会走弯路。