做这方向的人应该都有同感看到“基于阶梯碳交易的含P2G-CCS耦合和燃气掺氢的虚拟电厂优化调度”这个标题第一反应不是兴奋而是“又要开始处理一堆非线性和整数变量了”。但说句实话这种模型恰恰是现在虚拟电厂VPP研究里最值得啃的一块骨头——风电、光伏出力一波动弃风弃光靠什么消化燃气轮机的高碳排怎么压下来碳配额成本怎么传导到每个时段的运行决策这几个问题串起来就是一台完整的“电-气-氢-碳”多能流优化调度机器。这篇文章我打算把整个项目从建模到 Matlab 落地的思路完整走一遍重点放在阶梯碳交易机制如何线性化、P2G-CCS 耦合和燃气掺氢怎么放进约束、以及用 Yalmip 建模时哪些地方最容易翻车。适合正在做综合能源系统、虚拟电厂调度方向的研究生也适合刚接触碳交易建模、需要快速出一版可运行代码的工程师参考。1. 虚拟电厂为什么要把 P2G、CCS 和掺氢捆在一起先说个最直白的问题虚拟电厂本身就是一个聚合体风电、光伏、燃气轮机、储能、负荷捏在一起按整体最优来调度。既然已经有储能了为什么还要引入 P2G、CCS 和掺氢我的理解是这三张牌解决的是储能解决不了的三件事。1.1 三张技术牌的各自角色风电光伏大发时电多用不完最朴素的办法是弃掉或者给储能充电。但储能的容量和放电时长有限多出来的大量富余电力放不下这时候 P2G 的价值就出来了——电解水制氢把电能转化成氢能既可以存进氢罐也可以直接送到燃气轮机掺烧还能拿去和 CO2 反应制甲烷。P2G 实际上是把电力系统的“富余电量”转化成“气网/燃料侧的替代能源”相当于给电力系统开了一个额外的出口。CCS 的角色则完全相反它盯的是排放侧。燃气轮机烧天然气排出来的 CO2 是碳交易账单的主要来源。CCS 把这些 CO2 捕集下来净排放量降下来碳交易成本跟着降。如果只装 CCS捕下来的 CO2 还要考虑封存成本但如果在同一套系统里配上 P2G这些 CO2 可以直接送到甲烷化反应器和氢气合成甲烷重新变成燃料。这就是 P2G-CCS 耦合最核心的价值碳捕集不再是纯成本项而是给 P2G 提供碳源。燃气掺氢则是第三条降碳路径。天然气里掺入一定比例的氢气因为氢气燃烧不产生 CO2掺氢后的混合燃料碳排放强度会明显降低而且不需要对燃气轮机做大的改造——现有机组大多数能支持 10%~20% 的体积掺氢比例。从调度角度看掺氢比例本身就是一个连续决策变量相当于给系统加了一个“燃料侧降碳”的旋钮。1.2 耦合之后系统多了哪些可调度量把技术说清楚之后咱们再看调度模型。没有 P2G、CCS、掺氢时虚拟电厂的决策变量大概是燃气轮机出力、储能充放电、向上级电网买电卖电、弃风弃光量。加入这三套设备之后决策变量一下子增加了很多。P2G 侧有电解槽功率、氢气分流比例多少进储氢罐、多少直接送燃气轮机、多少去甲烷化。CCS 侧有捕集功率、捕集到的 CO2 量。掺氢侧有混合燃料中氢气的热值占比。燃气轮机侧除了出力大小还要决定燃料里天然气和氢气的配比。再加上储氢罐的动态变化、甲烷化反应器的运行功率整个问题从普通的机组组合变成多能流耦合优化。这些变量之间互相牵制牵制的纽带就是能量平衡和物料平衡。例如P2G 多产氢储氢罐的压力就上去了但燃气轮机的掺氢比例也有上限多出来的氢不能无限燃烧要么进罐要么去甲烷化。CCS 捕集的 CO2 多了甲烷化反应器如果开不了那么大功率CO2 就会富余富余的 CO2 要么储存要么排放。这些关系在数学上全是约束也是建模时最花时间的地方。1.3 碳政策是怎么变成调度指令的最后一个问题为什么偏偏要引入阶梯碳交易单一固定碳价也能反映碳成本但效果不够。固定碳价时只要边际减排成本低于固定碳价系统就会一直减一旦边际减排成本高于碳价调度就直接放弃减排该排就排。阶梯碳价的逻辑是递进式的——实际排放量超过免费配额越多单位碳排放的价格越贵直到减排带来的收益超过减排设备的运行成本。在这个机制下碳排放不是简单乘以一个常数加入目标函数而是一个随排放总量变化的分段线性函数。它会对调度产生非常直接的影响高排放时段燃气轮机出力会被压一压P2G 和 CCS 的出力会提一提掺氢比例也会尽量往上拉因为燃料侧降碳直接降低总排放量从而避开了阶梯碳价的高价区间。所以说这三套设备和阶梯碳交易不是各管各的它们共同构成了一个闭环新能源富余电力给 P2G 供能P2G 产氢给燃气轮机掺烧燃气轮机装配 CCS 捕集 CO2CO2 又回供给甲烷化反应器最终结果是碳排放总量下降、碳交易成本下降、弃风弃光率也降下来。而这个闭环的运行策略恰恰由阶梯碳价来定调。2. 阶梯碳交易机制怎么进优化模型阶梯碳交易在文献里很常见但真正动手写 Yalmip 代码时不少人会在这里卡住。原因很简单碳交易成本不是一个简单的线性函数而是一个分段函数而且分段点的斜率是递增的。处理不好要么模型非线性要么结果明显不符合预期。2.1 为什么是“阶梯”而不是固定碳价全国碳排放权交易市场的基本逻辑是政府给控排企业发放免费配额企业实际排放低于配额可以把富余配额拿出来卖实际排放高于配额需要购买缺口。阶梯的意思就是购买缺口这件事不是统一价而是分档计费——缺口在某个范围内的部分按低价算超过这个范围的部分按更高的价格算再往上还有第三档、第四档。举例来说假设免费配额是 A 吨第一档缺口是 0~a1 吨价格 c1第二档缺口是 a1~a2 吨价格 c2第三档是超过 a2 吨的部分价格 c3并且 c1 c2 c3。这样排放越多超出配额越多单位碳成本越贵呈现“多排多付、累进惩罚”的效果。这种结构的好处是调度目标里的碳成本会对高排放行为产生更强的抑制作用比固定碳价更贴近实际碳市场的运行特点。2.2 免费配额怎么给免费配额的计算方法直接影响模型结果。目前学术文献里常用的基准线法公式大致是E_a δ · P_VPP其中 δ 是单位电力碳排放基准强度tCO2/MWhP_VPP 是虚拟电厂的对外供电量。有的文献也会把 P2G 合成的天然气折算成对外供气量按供气和供电两条线发配额。我自己的习惯是取一个基准强度系数比如 δ 0.6~0.8 t/MWh再乘以虚拟电厂的总发电量。实际项目里依据当地碳配额核算规则进行替换即可模型结构不用改动只需要把 δ 换成实际标定值。这里想提醒一句免费配额如果给得太宽松阶梯碳价根本触发不了P2G-CCS 和掺氢的减排价值就体现不出来给得太紧系统为了减碳会过度投资设备运行成本总成本反而上升。做算例时建议把配额强度设成中间值再做一个敏感性分析观察碳价区间变化时各设备出力的变化趋势。2.3 目标函数里的碳成本表达式在优化模型里碳交易成本 C_carb 可以写成实际碳排放 E_total 和免费配额 E_a 之差的函数。设 d E_total - E_a则当 d ≤ 0 时配额富余系统可以通过卖配额获得收益碳交易成本为负值可简化按固定回购价 c0 计算当 d 0 时需要购买配额购买成本按阶梯价格计算。如果不做处理这是一个分段线性函数。但因为每一段斜率递增整体是一个凸函数所以存在一个很经典的线性化技巧碳交易成本 C_carb 会大于等于每一个分段线性函数在该点的取值。用数学表达就是C_carb ≥ c1 · d C_carb ≥ c1 · a1 c2 · (d - a1) C_carb ≥ c1 · a1 c2 · (a2 - a1) c3 · (d - a2)这里要求 d ≥ 0。如果 d 可能为负那就需要把配额富余的收益单独计算再和阶梯惩罚部分叠加。凸函数的这个性质允许我们用连续变量线性约束来逼近不需要引入二进制变量求解速度会快很多。2.4 阶梯碳价的线性化求解技巧实际写 Yalmip 代码时我会分两种情况处理。情况一排放配额缺口 d 始终大于 0即系统是净买配额方。这种情况下直接用上面三行不等式约束即可% 阶梯碳价参数 a1 30; a2 60; c1 80; c2 120; c3 160; % 设 E_total 和 E_a 均为 sdpvar 或常数d 为净排放缺口 d E_total - E_a; C_carb sdpvar(1, 1); constraints [constraints, C_carb c1 * d, C_carb c1 * a1 c2 * (d - a1), C_carb c1 * a1 c2 * (a2 - a1) c3 * (d - a2)];目标函数里加 C_carb最小化问题会自动让 C_carb 取到满足这些约束的最小值。因为 d 在第二段时第一条和第二条约约束同时有效在第三段时三条都有效C_carb 会自动落在正确的线性段上。情况二d 可正可负需要同时考虑配额收益。此时可以拆成 n_d max(d,0) 和 surplus max(-d,0)再用 big-M 和二进制变量表示分段归属。推荐用连续配额收益项加上二进制精确建模的阶梯购买项保证模型在非凸情况下也能准确描述碳成本。论文里很多直接用凸松弛算例里一般没问题但如果你要做极端场景测试建议用二进制版本兜底。M 1e4; z1 binvar(1,1); z2 binvar(1,1); z3 binvar(1,1); n1 sdpvar(1,1); n2 sdpvar(1,1); n3 sdpvar(1,1); % 净缺口非负 d_pos max(d, 0); constraints [constraints, n1 0, n1 a1 * z1, n2 0, n2 a2 * z2, n3 0, n3 M * z3, d_pos n1 n2 n3, sum(z1 z2 z3) 1]; C_carb c1 * n1 c2 * n2 c3 * n3 ... - c0 * max(-d, 0);注意这段代码里的 z1、z2、z3 要和 d_pos 落到的具体区间对应实际使用时要对边界条件做更严格的定义。我这里展示的是建模思路跑正式算例前建议先用小规模测试数据验证分段边界是否正确。3. 系统结构与优化调度模型模型的核心是把 VPP 内部各条能量流的关系理清楚。我设计的结构是风电、光伏作为不可控但可削减的电源燃气轮机作为可控主力机组支持掺氢P2G 电解槽作为电力消纳和氢气生产设备配储氢罐CCS 捕集燃气轮机烟气中的 CO2捕集到的 CO2 供给甲烷化反应器甲烷化反应器把氢气和 CO2 合成为甲烷回补气网或直接供燃气轮机使用。3.1 VPP 内部能量流关系先看电功率平衡这是最基本的约束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_cha(t) P_sell(t) P_curt(t)其中 P_wt、P_pv 是风电光伏实际出力和削减后并入系统的功率P_gt 是燃气轮机净电出力P_buy 和 P_sell 是向上级电网购电和售电P_p2g 是电解槽耗电功率P_ccs 是碳捕集系统耗电功率P_curt 是弃风弃光功率P_cha/P_dis 是储能充放电功率。氢能平衡稍微复杂一点。P2G 电解槽产氢速率记为 F_h2(t)按能量单位 MW 计产出的氢气有三个去处送入储氢罐记为 F_h2_st(t)直接供燃气轮机掺氢燃烧记为 F_h2_gt(t)送入甲烷化反应器记为 F_h2_meth(t)。所以F_h2(t) F_h2_st(t) F_h2_gt(t) F_h2_meth(t)储氢罐的动态方程S_h2(t1) S_h2(t) η_st · F_h2_st(t) - F_out(t)注意不是所有的制氢都进罐也不是罐里的氢随时都能无限流出要加充放氢速率上限和容量上下限约束。我习惯在模型里加一个调度周期首末储氢量相等约束否则求解器会把储氢罐在最后时段“搬空”结果不物理。3.2 目标函数综合成本最小化目标函数是所有成本的加总。我的定义如下C_total C_elec C_fuel C_om C_carb C_curt其中C_elec Σ[price_buy(t)·P_buy(t) - price_sell(t)·P_sell(t)]购电成本和售电收益C_fuel Σ[price_NG(t)·V_ng(t)]天然气购气成本掺氢后耗气量按热值换算C_om Σ[k_gt·P_gt(t) k_p2g·P_p2g(t) k_ccs·P_ccs(t)]设备运行维护成本C_carb 是阶梯碳交易成本C_curt Σ[k_curt·(P_wt_pre(t)-P_wt(t)) k_curt·(P_pv_pre(t)-P_pv(t))]弃风弃光惩罚。碳交易成本这一项是标题的核心前面已经单独分析过。其它几项虽然简单但参数给得合理与否直接决定结果好坏——尤其是天然气价格和电网购电价格如果峰谷价差不够大P2G 在谷时制氢、在峰时发电/售气的套利空间就出不来。3.3 约束条件里的四个关键耦合模型里最需要仔细写的是下面这四个耦合约束。第一个是燃气轮机掺氢约束。我采用热值掺氢比 h_ratio(t) 来描述掺氢程度上限设为 15%Q_gt(t) Q_ng(t) Q_h2(t) Q_h2(t) h_ratio(t) · Q_gt(t) Q_ng(t) (1 - h_ratio(t)) · Q_gt(t) 0 ≤ h_ratio(t) ≤ H_MAX其中 Q_gt(t) 是燃气轮机消耗的总燃料热量MWQ_ng 是天然气贡献的热量Q_h2 是氢气贡献的热量。氢气的碳排放因子为零所以燃气轮机实际碳排放量按 Q_ng 计算E_emit(t) EF_ng · Q_ng(t)第二个是 CCS 捕集约束。CCS 不是想捕多少就捕多少捕集量受设备容量和运行功率限制E_capture(t) η_ccs · P_ccs_norm(t)·E_emit(t) 0 ≤ P_ccs_norm(t) ≤ 1 E_capture(t) ≤ P_ccs_max · Δt这里的 P_ccs_norm(t) 是 CCS 归一化运行水平可以理解成捕集装置在 t 时段的负荷率。捕集等量 CO2 需要消耗电能P_ccs(t) β_ccs · E_capture(t)β_ccs 的单位是 MWh/tCO2我算例里常用 0.2~0.4 MWh/t。别小看这个耗电项它会进入电功率平衡约束直接影响系统对风电光伏的消纳能力。第三个是 P2G-CCS 耦合约束。CCS 捕集到的 CO2 并不一定全部被利用甲烷化反应器用多少 CO2 取决于氢气的供应量和反应器容量E_co2_used(t) ≤ E_capture(t) E_co2_st(t) E_co2_used(t) ≤ CO2_reactor_max · Δt进入甲烷化反应器的 CO2 和氢气满足反应比例关系。我按能量转化效率处理甲烷化产出的甲烷功率 η_meth · F_h2_meth(t)同时消耗的 CO2 要满足化学计量关系换算成简化式为E_co2_used(t) γ_co2_per_mwh · Q_meth(t)γ_co2_per_mwh 是每 MW·h 甲烷对应的 CO2 耗量t。这样既保持线性又能量化“捕集下来的 CO2 是否够用”这个关键约束。第四个是阶梯碳交易约束。实际碳排放总量 E_total Σ E_emit(t) - Σ E_capture(t)用这个总量去和免费配额 E_a 比较再代入前文的分段函数。4. 用 Matlab 落地Yalmip 建模与求解器选择很多同学在这个方向卡住不是因为模型不会列而是 Matlab 环境、Yalmip 语法、求解器之间的配合不熟练。下面按我的实操顺序来讲。4.1 环境准备Matlab Yalmip Gurobi做优化调度我的首选组合是 Matlab R2020 以上版本 Yalmip 最新版 Gurobi 求解器。Yalmip 是一个免费的建模工具箱提供 sdpvar、binvar、optimize 等接口特点是非常适合论文里常见的优化问题建模约束可以直接用表达式累加不用手动拼大矩阵。Gurobi 对学术用户免费性能和稳定性比 Matlab 内置的 intlinprog 好不少尤其是混合整数线性规划问题规模上去之后差距非常明显。安装时记得把 Gurobi 的 Matlab 接口路径添加到 Matlab 路径里并在代码开头用 yalmiptest 验证一下求解器是否被识别。Yalmip 的默认输出里能看到当前可用的求解器列表确认 Gurobi 出现在列表中再往下写。如果你暂时没有 Gurobi也可以用 Matlab 内置的 intlinprog 跑小规模算例。Yalmip 里只要设置 sdpsettings(solver,intlinprog) 即可切换。只是模型一旦超过几百个变量和几百个约束intlinprog 的速度和稳定性就很难让人满意。4.2 基础数据与变量定义我习惯把 24 小时作为调度周期时间间隔 Δt1h。首先定义时间序列数据风电预测出力、光伏预测出力、负荷曲线、分时电价、气价、风光削减惩罚系数等。T 24; dt 1; % 预测数据示例按实际算例替换 P_wt_pre [50, 48, 45, 42, 40, ...]; P_pv_pre [0, 0, 0, 0, 0, ...]; P_load [70, 68, 65, 62, 60, ...]; % 电价峰段、平段、谷段 price_buy ones(1, T) * 0.5; price_buy(9:12) 0.9; price_buy(18:21) 0.9; price_sell price_buy * 0.8; price_NG ones(1, T) * 0.35; % 元/kWh热值口径然后定义决策变量。Yalmip 里连续变量用 sdpvar整数变量用 binvar/intvarP_wt sdpvar(1, T, full); P_pv sdpvar(1, T, full); P_gt sdpvar(1, T, full); P_buy sdpvar(1, T, full); P_sell sdpvar(1, T, full); P_p2g sdpvar(1, T, full); P_ccs sdpvar(1, T, full); F_h2_gt sdpvar(1, T, full); F_h2_st sdpvar(1, T, full); F_h2_meth sdpvar(1, T, full); S_h2 sdpvar(1, T1, full); h_ratio sdpvar(1, T, full); E_emit sdpvar(1, T, full); E_capture sdpvar(1, T, full); P_curt sdpvar(1, T, full);sdpvar(1,T,full)里的 full 表示这是一个 1×T 的普通矩阵而不是对称矩阵。不写 full 时 Yalmip 默认创建对称方阵非常容易踩坑。4.3 各模块约束的 Yalmip 写法变量定义完之后建议用 cell 数组或矩阵拼接来累积约束。我常用的是一个 constraints 数组从空开始逐步 appendconstraints []; % 电功率平衡 constraints [constraints, P_wt P_pv P_gt P_buy P_load P_p2g P_ccs P_sell P_curt];风电光伏的实际出力不能超过预测值削减量非负这是弃风弃光的来源constraints [constraints, 0 P_wt P_wt_pre, 0 P_pv P_pv_pre, P_curt (P_wt_pre - P_wt) (P_pv_pre - P_pv)];燃气轮机的出力和燃料热量之间建立关系。假设机组效率 η_gt燃料总热量 Q_gt P_gt / η_gt。然后掺氢比约束eta_gt 0.45; H_MAX 0.15; Q_gt P_gt / eta_gt; Q_ng (1 - h_ratio) .* Q_gt; Q_h2 h_ratio .* Q_gt; constraints [constraints, 0 P_gt P_gt_max, 0 h_ratio H_MAX];P2G 产氢eta_p2g 0.65; F_h2 eta_p2g * P_p2g; constraints [constraints, 0 P_p2g P_p2g_max, F_h2 F_h2_gt F_h2_st F_h2_meth];CCS 捕集与耗电beta_ccs 0.25; % MWh/tCO2 eta_ccs 0.90; % 捕集率 constraints [constraints, E_capture eta_ccs * E_emit, P_ccs beta_ccs * E_capture, 0 E_capture E_capture_max];碳排放量按天然气消耗量计算EF_ng 0.2; % tCO2/MWh天然气低位热值口径 constraints [constraints, E_emit EF_ng * Q_ng];储氢罐S_h2_min 0; S_h2_max 80; constraints [constraints, S_h2(1) S_h2_init, S_h2(T1) S_h2_init, 0 S_h2 S_h2_max, S_h2(2:T1) S_h2(1:T) eta_st * F_h2_st - F_h2_gt];注意最后一行储氢动态可以把 F_h2_gt 看成从储氢罐流向燃气轮机的流量但如果你希望“P2G 直接产氢即产即用不需要经过储氢罐”就要把平衡式拆开。我的做法是F_h2_gt 一部分直接来自电解槽一部分来自储氢罐二者通过另一个分流变量控制。简化版里可以直接按“罐进罐出”建模做机理验证够用。4.4 求解与结果提取约束搭好、目标函数写好后调用求解器的代码如下ops sdpsettings(solver, gurobi, verbose, 2, showprogress, 1); sol optimize(constraints, Objective, ops); if sol.problem 0 disp(求解成功); else disp([求解失败: , sol.info]); end P_wt_opt value(P_wt); P_gt_opt value(P_gt); h_ratio_opt value(h_ratio); E_capture_opt value(E_capture); C_carb_opt value(C_carb);5. 典型日算例设计与结果分析算例设计直接决定论文或项目能不能讲出自洽的故事。我建议至少设置 4 个场景来对比这样 P2G-CCS 耦合和掺氢的价值能被单独拆出来看。5.1 场景设置场景一是基准场景不带 P2G、CCS、掺氢燃气轮机纯烧天然气弃风弃光只能靠储能和向上级电网售电来缓解。场景二在基准上加入 P2G 和 CCS但不掺氢用来考察“电-气-碳耦合”单独作用的效果。场景三加入燃气掺氢但碳价是固定值而不是阶梯值用来对比阶梯碳价和固定碳价的差别。场景四用完整模型阶梯碳价 P2G CCS 掺氢全上。这种对比矩阵可以把每个模块的边际贡献都拆出来而且在论文里也很容易画成表格或者柱状图。5.2 算例参数建议下面是我常用的一组参数供参考设备/参数数值风电场容量120 MW光伏电站容量80 MW燃气轮机容量60 MW燃气轮机效率45%P2G 电解槽容量50 MWP2G 制氢效率65%储氢罐容量80 MWh能量口径掺氢比上限15%热值比CCS 捕集率90%CCS 耗电系数0.25 MWh/tCO2甲烷化效率70%免费配额强度0.65 t/MWh阶梯区间30 t / 60 t阶梯碳价80 / 120 / 160 元/t注意这里的所有能量单位都统一到 MWh碳单位统一到吨。如果不统一后面的约束很容易出现数量级错误。我自己在调试时吃过亏储氢罐和甲烷化用不同的能量口径导致氢平衡怎么调都不满足最后发现是单位没换算。5.3 从模型机理看的典型结论我按照这套模型跑下来的经验有几个比较稳定的趋势你可以拿自己的算例验证。第一阶梯碳价下燃气轮机的高出力时段会被明显压缩。尤其是傍晚负荷高峰和电价高峰时段燃气轮机本来应该满发但因为碳排放缺口已经逼近第二档甚至第三档继续高发会触发高碳价调度会把一部分负荷转向上级电网购电。这个转变在场景一里完全看不到因为固定碳价不会产生这种累积惩罚效应。第二掺氢比例在负荷高峰时段倾向于走到上限。原因很直接高峰时段燃气轮机出力大天然气消耗量大排放量也大提高掺氢比例可以压低天然气消耗从而压低碳排放。这个收益在阶梯碳价下比固定碳价更明显因为减少的排放处在高价区间的边际价值更高。第三P2G 的运行策略受“电-氢-气”价差驱动。凌晨风电大发、电价低谷时P2G 出力高制出来的氢一部分进储氢罐一部分用于甲烷化白天电价高峰时段储氢罐放氢和天然气掺混发电甲烷化也优先在碳价压力大的时段运行。如果气价和电价没有足够大的峰谷差P2G 的收益率会明显下降这在参数敏感性分析里能看得很清楚。第四CCS 捕集到 CO2 能否被完全利用取决于甲烷化反应器的容量和氢气供应量。算例里经常出现的情况是CCS 在燃气轮机满发时捕集了大量 CO2但因为 P2G 产氢量不够或者反应器容量限制CO2 用不完只能储存在 CO2 储罐或封存。这时如果模型里没有 CO2 储存约束就等效于“免费排放”结果会过于乐观。6. 调试经验与常见坑最后聊聊我在调试此类模型时踩过的一些坑都是真实花费过时间的地方希望能帮你少走弯路。6.1 掺氢口径不一致导致热值错乱这是最隐蔽的坑。工程上燃气掺氢常用体积比比如掺氢 15%意思是混合燃气中氢气体积占 15%。但相同体积下天然气的热值比氢气高不少直接用体积比乘以燃料热值会低估天然气贡献。我的建议是在模型里用热值掺氢比换算公式是β_energy β_vol · LHV_H2 / (β_vol · LHV_H2 (1 - β_vol) · LHV_NG)其中 LHV_NG 约 36 MJ/m³LHV_H2 约 10.8 MJ/m³。做参数标定时先用这个公式把体积比换算成热值比再代入模型。否则碳排放强度、天然气耗量全都会偏离真实值。6.2 CCS 耗电项被漏掉很多初版模型容易忘记把 P_ccs 放进电功率平衡方程。漏掉之后的结果是CCS 仿佛在不耗电地捕碳捕集量再大也不影响电力调度系统会过度使用 CCS碳捕集成本被严重低估。我习惯在每个设备运行约束旁边同步检查它在哪些平衡方程里出现——P2G 出现在电平衡和氢平衡CCS 出现在电平衡和碳平衡燃气轮机出现在电平衡、气平衡和碳平衡。6.3 CO2 物料守恒约束过于松散如果模型里没有约束“甲烷化消耗的 CO2 不能超过 CCS 捕集量和 CO2 储罐放出量之和”求解器会私自“凭空变出 CO2”让甲烷化白白产气结果特别好看但不物理。这个问题的排查方法很简单把结果里的 CO2 利用量拉出来看看它是否等于捕集量如果不等于说明耦合约束没生效。6.4 储氢罐首末状态约束为了避免调度策略把储氢罐当成免费垃圾桶首末状态约束必须加S_h2(1) S_h2(T1)。不加的话模型会在最后一个时段把所有储氢清空把氢全部烧掉或者卖掉人为拉低总成本。同样的问题也会出现在 CO2 储罐和储能电池上建议所有有存储环节的统一加首末平衡约束。6.5 Yalmip 求解失败时的调试顺序遇到 infeasible 时不要慌按下面的顺序排查通常能定位问题。先用 sdpsettings(verbose,1) 看求解器输出的是哪类错误。如果是“infeasible”把目标函数临时设成常数 0先求可行解如果还是不可行说明约束本身就矛盾。然后逐条“注释”约束同时固定变量取值看哪条约束和其余约束冲突。最容易出问题的通常是电平衡和储氢动态方程尤其是下标维度对不上时Yalmip 不会立刻报错但会生成一个形状不对的约束矩阵。另外Yalmip 的 sdpvar 定义中一个常见错误是忘记 full导致变量变成对称矩阵。写代码时建议用 assert 检查每个变量的 sizeassert(isequal(size(P_wt_pre), [1, T]), P_wt_pre dimension mismatch);检查一遍所有预测数据的维度做数据驱动建模时能省很多事。我在实际项目里还有一个习惯把整个调度模型封装成一个 function输入是风光预测、电价、气价、碳价、设备参数输出是各设备的出力计划、总成本和碳排放量。这样跑敏感性分析和场景对比时不需要反复复制代码只需要在循环里改参数调用函数即可效率和可维护性都会高很多。function [result] vpp_schedule(P_wt_pre, P_pv_pre, P_load, price_buy, price_NG, cparam) % 该函数返回调度结果结构体 result包含各出力序列和成本 ... end参数敏感性分析是这个方向绕不开的一步也是论文里最出彩的部分。我通常会对碳价曲线、掺氢上限、CCS 耗电系数三个参数分别做多组扫描观察总成本、碳排放强度、P2G 利用小时数三个指标的变化这些结果比单独的算例更有说服力。跑完模型再回头看这套系统的核心其实不是某个单一设备而是碳价格信号如何同时影响电力调度、燃料配比和碳捕集水平这三层决策。把阶梯碳价这个关键信号建模清楚了P2G、CCS 和掺氢的调度策略就都跟着顺理成章了。 基于阶梯碳交易的含P2G-CCS耦合与燃气掺氢虚拟电厂优化调度
做这方向的人应该都有同感看到“基于阶梯碳交易的含P2G-CCS耦合和燃气掺氢的虚拟电厂优化调度”这个标题第一反应不是兴奋而是“又要开始处理一堆非线性和整数变量了”。但说句实话这种模型恰恰是现在虚拟电厂VPP研究里最值得啃的一块骨头——风电、光伏出力一波动弃风弃光靠什么消化燃气轮机的高碳排怎么压下来碳配额成本怎么传导到每个时段的运行决策这几个问题串起来就是一台完整的“电-气-氢-碳”多能流优化调度机器。这篇文章我打算把整个项目从建模到 Matlab 落地的思路完整走一遍重点放在阶梯碳交易机制如何线性化、P2G-CCS 耦合和燃气掺氢怎么放进约束、以及用 Yalmip 建模时哪些地方最容易翻车。适合正在做综合能源系统、虚拟电厂调度方向的研究生也适合刚接触碳交易建模、需要快速出一版可运行代码的工程师参考。1. 虚拟电厂为什么要把 P2G、CCS 和掺氢捆在一起先说个最直白的问题虚拟电厂本身就是一个聚合体风电、光伏、燃气轮机、储能、负荷捏在一起按整体最优来调度。既然已经有储能了为什么还要引入 P2G、CCS 和掺氢我的理解是这三张牌解决的是储能解决不了的三件事。1.1 三张技术牌的各自角色风电光伏大发时电多用不完最朴素的办法是弃掉或者给储能充电。但储能的容量和放电时长有限多出来的大量富余电力放不下这时候 P2G 的价值就出来了——电解水制氢把电能转化成氢能既可以存进氢罐也可以直接送到燃气轮机掺烧还能拿去和 CO2 反应制甲烷。P2G 实际上是把电力系统的“富余电量”转化成“气网/燃料侧的替代能源”相当于给电力系统开了一个额外的出口。CCS 的角色则完全相反它盯的是排放侧。燃气轮机烧天然气排出来的 CO2 是碳交易账单的主要来源。CCS 把这些 CO2 捕集下来净排放量降下来碳交易成本跟着降。如果只装 CCS捕下来的 CO2 还要考虑封存成本但如果在同一套系统里配上 P2G这些 CO2 可以直接送到甲烷化反应器和氢气合成甲烷重新变成燃料。这就是 P2G-CCS 耦合最核心的价值碳捕集不再是纯成本项而是给 P2G 提供碳源。燃气掺氢则是第三条降碳路径。天然气里掺入一定比例的氢气因为氢气燃烧不产生 CO2掺氢后的混合燃料碳排放强度会明显降低而且不需要对燃气轮机做大的改造——现有机组大多数能支持 10%~20% 的体积掺氢比例。从调度角度看掺氢比例本身就是一个连续决策变量相当于给系统加了一个“燃料侧降碳”的旋钮。1.2 耦合之后系统多了哪些可调度量把技术说清楚之后咱们再看调度模型。没有 P2G、CCS、掺氢时虚拟电厂的决策变量大概是燃气轮机出力、储能充放电、向上级电网买电卖电、弃风弃光量。加入这三套设备之后决策变量一下子增加了很多。P2G 侧有电解槽功率、氢气分流比例多少进储氢罐、多少直接送燃气轮机、多少去甲烷化。CCS 侧有捕集功率、捕集到的 CO2 量。掺氢侧有混合燃料中氢气的热值占比。燃气轮机侧除了出力大小还要决定燃料里天然气和氢气的配比。再加上储氢罐的动态变化、甲烷化反应器的运行功率整个问题从普通的机组组合变成多能流耦合优化。这些变量之间互相牵制牵制的纽带就是能量平衡和物料平衡。例如P2G 多产氢储氢罐的压力就上去了但燃气轮机的掺氢比例也有上限多出来的氢不能无限燃烧要么进罐要么去甲烷化。CCS 捕集的 CO2 多了甲烷化反应器如果开不了那么大功率CO2 就会富余富余的 CO2 要么储存要么排放。这些关系在数学上全是约束也是建模时最花时间的地方。1.3 碳政策是怎么变成调度指令的最后一个问题为什么偏偏要引入阶梯碳交易单一固定碳价也能反映碳成本但效果不够。固定碳价时只要边际减排成本低于固定碳价系统就会一直减一旦边际减排成本高于碳价调度就直接放弃减排该排就排。阶梯碳价的逻辑是递进式的——实际排放量超过免费配额越多单位碳排放的价格越贵直到减排带来的收益超过减排设备的运行成本。在这个机制下碳排放不是简单乘以一个常数加入目标函数而是一个随排放总量变化的分段线性函数。它会对调度产生非常直接的影响高排放时段燃气轮机出力会被压一压P2G 和 CCS 的出力会提一提掺氢比例也会尽量往上拉因为燃料侧降碳直接降低总排放量从而避开了阶梯碳价的高价区间。所以说这三套设备和阶梯碳交易不是各管各的它们共同构成了一个闭环新能源富余电力给 P2G 供能P2G 产氢给燃气轮机掺烧燃气轮机装配 CCS 捕集 CO2CO2 又回供给甲烷化反应器最终结果是碳排放总量下降、碳交易成本下降、弃风弃光率也降下来。而这个闭环的运行策略恰恰由阶梯碳价来定调。2. 阶梯碳交易机制怎么进优化模型阶梯碳交易在文献里很常见但真正动手写 Yalmip 代码时不少人会在这里卡住。原因很简单碳交易成本不是一个简单的线性函数而是一个分段函数而且分段点的斜率是递增的。处理不好要么模型非线性要么结果明显不符合预期。2.1 为什么是“阶梯”而不是固定碳价全国碳排放权交易市场的基本逻辑是政府给控排企业发放免费配额企业实际排放低于配额可以把富余配额拿出来卖实际排放高于配额需要购买缺口。阶梯的意思就是购买缺口这件事不是统一价而是分档计费——缺口在某个范围内的部分按低价算超过这个范围的部分按更高的价格算再往上还有第三档、第四档。举例来说假设免费配额是 A 吨第一档缺口是 0~a1 吨价格 c1第二档缺口是 a1~a2 吨价格 c2第三档是超过 a2 吨的部分价格 c3并且 c1 c2 c3。这样排放越多超出配额越多单位碳成本越贵呈现“多排多付、累进惩罚”的效果。这种结构的好处是调度目标里的碳成本会对高排放行为产生更强的抑制作用比固定碳价更贴近实际碳市场的运行特点。2.2 免费配额怎么给免费配额的计算方法直接影响模型结果。目前学术文献里常用的基准线法公式大致是E_a δ · P_VPP其中 δ 是单位电力碳排放基准强度tCO2/MWhP_VPP 是虚拟电厂的对外供电量。有的文献也会把 P2G 合成的天然气折算成对外供气量按供气和供电两条线发配额。我自己的习惯是取一个基准强度系数比如 δ 0.6~0.8 t/MWh再乘以虚拟电厂的总发电量。实际项目里依据当地碳配额核算规则进行替换即可模型结构不用改动只需要把 δ 换成实际标定值。这里想提醒一句免费配额如果给得太宽松阶梯碳价根本触发不了P2G-CCS 和掺氢的减排价值就体现不出来给得太紧系统为了减碳会过度投资设备运行成本总成本反而上升。做算例时建议把配额强度设成中间值再做一个敏感性分析观察碳价区间变化时各设备出力的变化趋势。2.3 目标函数里的碳成本表达式在优化模型里碳交易成本 C_carb 可以写成实际碳排放 E_total 和免费配额 E_a 之差的函数。设 d E_total - E_a则当 d ≤ 0 时配额富余系统可以通过卖配额获得收益碳交易成本为负值可简化按固定回购价 c0 计算当 d 0 时需要购买配额购买成本按阶梯价格计算。如果不做处理这是一个分段线性函数。但因为每一段斜率递增整体是一个凸函数所以存在一个很经典的线性化技巧碳交易成本 C_carb 会大于等于每一个分段线性函数在该点的取值。用数学表达就是C_carb ≥ c1 · d C_carb ≥ c1 · a1 c2 · (d - a1) C_carb ≥ c1 · a1 c2 · (a2 - a1) c3 · (d - a2)这里要求 d ≥ 0。如果 d 可能为负那就需要把配额富余的收益单独计算再和阶梯惩罚部分叠加。凸函数的这个性质允许我们用连续变量线性约束来逼近不需要引入二进制变量求解速度会快很多。2.4 阶梯碳价的线性化求解技巧实际写 Yalmip 代码时我会分两种情况处理。情况一排放配额缺口 d 始终大于 0即系统是净买配额方。这种情况下直接用上面三行不等式约束即可% 阶梯碳价参数 a1 30; a2 60; c1 80; c2 120; c3 160; % 设 E_total 和 E_a 均为 sdpvar 或常数d 为净排放缺口 d E_total - E_a; C_carb sdpvar(1, 1); constraints [constraints, C_carb c1 * d, C_carb c1 * a1 c2 * (d - a1), C_carb c1 * a1 c2 * (a2 - a1) c3 * (d - a2)];目标函数里加 C_carb最小化问题会自动让 C_carb 取到满足这些约束的最小值。因为 d 在第二段时第一条和第二条约约束同时有效在第三段时三条都有效C_carb 会自动落在正确的线性段上。情况二d 可正可负需要同时考虑配额收益。此时可以拆成 n_d max(d,0) 和 surplus max(-d,0)再用 big-M 和二进制变量表示分段归属。推荐用连续配额收益项加上二进制精确建模的阶梯购买项保证模型在非凸情况下也能准确描述碳成本。论文里很多直接用凸松弛算例里一般没问题但如果你要做极端场景测试建议用二进制版本兜底。M 1e4; z1 binvar(1,1); z2 binvar(1,1); z3 binvar(1,1); n1 sdpvar(1,1); n2 sdpvar(1,1); n3 sdpvar(1,1); % 净缺口非负 d_pos max(d, 0); constraints [constraints, n1 0, n1 a1 * z1, n2 0, n2 a2 * z2, n3 0, n3 M * z3, d_pos n1 n2 n3, sum(z1 z2 z3) 1]; C_carb c1 * n1 c2 * n2 c3 * n3 ... - c0 * max(-d, 0);注意这段代码里的 z1、z2、z3 要和 d_pos 落到的具体区间对应实际使用时要对边界条件做更严格的定义。我这里展示的是建模思路跑正式算例前建议先用小规模测试数据验证分段边界是否正确。3. 系统结构与优化调度模型模型的核心是把 VPP 内部各条能量流的关系理清楚。我设计的结构是风电、光伏作为不可控但可削减的电源燃气轮机作为可控主力机组支持掺氢P2G 电解槽作为电力消纳和氢气生产设备配储氢罐CCS 捕集燃气轮机烟气中的 CO2捕集到的 CO2 供给甲烷化反应器甲烷化反应器把氢气和 CO2 合成为甲烷回补气网或直接供燃气轮机使用。3.1 VPP 内部能量流关系先看电功率平衡这是最基本的约束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_cha(t) P_sell(t) P_curt(t)其中 P_wt、P_pv 是风电光伏实际出力和削减后并入系统的功率P_gt 是燃气轮机净电出力P_buy 和 P_sell 是向上级电网购电和售电P_p2g 是电解槽耗电功率P_ccs 是碳捕集系统耗电功率P_curt 是弃风弃光功率P_cha/P_dis 是储能充放电功率。氢能平衡稍微复杂一点。P2G 电解槽产氢速率记为 F_h2(t)按能量单位 MW 计产出的氢气有三个去处送入储氢罐记为 F_h2_st(t)直接供燃气轮机掺氢燃烧记为 F_h2_gt(t)送入甲烷化反应器记为 F_h2_meth(t)。所以F_h2(t) F_h2_st(t) F_h2_gt(t) F_h2_meth(t)储氢罐的动态方程S_h2(t1) S_h2(t) η_st · F_h2_st(t) - F_out(t)注意不是所有的制氢都进罐也不是罐里的氢随时都能无限流出要加充放氢速率上限和容量上下限约束。我习惯在模型里加一个调度周期首末储氢量相等约束否则求解器会把储氢罐在最后时段“搬空”结果不物理。3.2 目标函数综合成本最小化目标函数是所有成本的加总。我的定义如下C_total C_elec C_fuel C_om C_carb C_curt其中C_elec Σ[price_buy(t)·P_buy(t) - price_sell(t)·P_sell(t)]购电成本和售电收益C_fuel Σ[price_NG(t)·V_ng(t)]天然气购气成本掺氢后耗气量按热值换算C_om Σ[k_gt·P_gt(t) k_p2g·P_p2g(t) k_ccs·P_ccs(t)]设备运行维护成本C_carb 是阶梯碳交易成本C_curt Σ[k_curt·(P_wt_pre(t)-P_wt(t)) k_curt·(P_pv_pre(t)-P_pv(t))]弃风弃光惩罚。碳交易成本这一项是标题的核心前面已经单独分析过。其它几项虽然简单但参数给得合理与否直接决定结果好坏——尤其是天然气价格和电网购电价格如果峰谷价差不够大P2G 在谷时制氢、在峰时发电/售气的套利空间就出不来。3.3 约束条件里的四个关键耦合模型里最需要仔细写的是下面这四个耦合约束。第一个是燃气轮机掺氢约束。我采用热值掺氢比 h_ratio(t) 来描述掺氢程度上限设为 15%Q_gt(t) Q_ng(t) Q_h2(t) Q_h2(t) h_ratio(t) · Q_gt(t) Q_ng(t) (1 - h_ratio(t)) · Q_gt(t) 0 ≤ h_ratio(t) ≤ H_MAX其中 Q_gt(t) 是燃气轮机消耗的总燃料热量MWQ_ng 是天然气贡献的热量Q_h2 是氢气贡献的热量。氢气的碳排放因子为零所以燃气轮机实际碳排放量按 Q_ng 计算E_emit(t) EF_ng · Q_ng(t)第二个是 CCS 捕集约束。CCS 不是想捕多少就捕多少捕集量受设备容量和运行功率限制E_capture(t) η_ccs · P_ccs_norm(t)·E_emit(t) 0 ≤ P_ccs_norm(t) ≤ 1 E_capture(t) ≤ P_ccs_max · Δt这里的 P_ccs_norm(t) 是 CCS 归一化运行水平可以理解成捕集装置在 t 时段的负荷率。捕集等量 CO2 需要消耗电能P_ccs(t) β_ccs · E_capture(t)β_ccs 的单位是 MWh/tCO2我算例里常用 0.2~0.4 MWh/t。别小看这个耗电项它会进入电功率平衡约束直接影响系统对风电光伏的消纳能力。第三个是 P2G-CCS 耦合约束。CCS 捕集到的 CO2 并不一定全部被利用甲烷化反应器用多少 CO2 取决于氢气的供应量和反应器容量E_co2_used(t) ≤ E_capture(t) E_co2_st(t) E_co2_used(t) ≤ CO2_reactor_max · Δt进入甲烷化反应器的 CO2 和氢气满足反应比例关系。我按能量转化效率处理甲烷化产出的甲烷功率 η_meth · F_h2_meth(t)同时消耗的 CO2 要满足化学计量关系换算成简化式为E_co2_used(t) γ_co2_per_mwh · Q_meth(t)γ_co2_per_mwh 是每 MW·h 甲烷对应的 CO2 耗量t。这样既保持线性又能量化“捕集下来的 CO2 是否够用”这个关键约束。第四个是阶梯碳交易约束。实际碳排放总量 E_total Σ E_emit(t) - Σ E_capture(t)用这个总量去和免费配额 E_a 比较再代入前文的分段函数。4. 用 Matlab 落地Yalmip 建模与求解器选择很多同学在这个方向卡住不是因为模型不会列而是 Matlab 环境、Yalmip 语法、求解器之间的配合不熟练。下面按我的实操顺序来讲。4.1 环境准备Matlab Yalmip Gurobi做优化调度我的首选组合是 Matlab R2020 以上版本 Yalmip 最新版 Gurobi 求解器。Yalmip 是一个免费的建模工具箱提供 sdpvar、binvar、optimize 等接口特点是非常适合论文里常见的优化问题建模约束可以直接用表达式累加不用手动拼大矩阵。Gurobi 对学术用户免费性能和稳定性比 Matlab 内置的 intlinprog 好不少尤其是混合整数线性规划问题规模上去之后差距非常明显。安装时记得把 Gurobi 的 Matlab 接口路径添加到 Matlab 路径里并在代码开头用 yalmiptest 验证一下求解器是否被识别。Yalmip 的默认输出里能看到当前可用的求解器列表确认 Gurobi 出现在列表中再往下写。如果你暂时没有 Gurobi也可以用 Matlab 内置的 intlinprog 跑小规模算例。Yalmip 里只要设置 sdpsettings(solver,intlinprog) 即可切换。只是模型一旦超过几百个变量和几百个约束intlinprog 的速度和稳定性就很难让人满意。4.2 基础数据与变量定义我习惯把 24 小时作为调度周期时间间隔 Δt1h。首先定义时间序列数据风电预测出力、光伏预测出力、负荷曲线、分时电价、气价、风光削减惩罚系数等。T 24; dt 1; % 预测数据示例按实际算例替换 P_wt_pre [50, 48, 45, 42, 40, ...]; P_pv_pre [0, 0, 0, 0, 0, ...]; P_load [70, 68, 65, 62, 60, ...]; % 电价峰段、平段、谷段 price_buy ones(1, T) * 0.5; price_buy(9:12) 0.9; price_buy(18:21) 0.9; price_sell price_buy * 0.8; price_NG ones(1, T) * 0.35; % 元/kWh热值口径然后定义决策变量。Yalmip 里连续变量用 sdpvar整数变量用 binvar/intvarP_wt sdpvar(1, T, full); P_pv sdpvar(1, T, full); P_gt sdpvar(1, T, full); P_buy sdpvar(1, T, full); P_sell sdpvar(1, T, full); P_p2g sdpvar(1, T, full); P_ccs sdpvar(1, T, full); F_h2_gt sdpvar(1, T, full); F_h2_st sdpvar(1, T, full); F_h2_meth sdpvar(1, T, full); S_h2 sdpvar(1, T1, full); h_ratio sdpvar(1, T, full); E_emit sdpvar(1, T, full); E_capture sdpvar(1, T, full); P_curt sdpvar(1, T, full);sdpvar(1,T,full)里的 full 表示这是一个 1×T 的普通矩阵而不是对称矩阵。不写 full 时 Yalmip 默认创建对称方阵非常容易踩坑。4.3 各模块约束的 Yalmip 写法变量定义完之后建议用 cell 数组或矩阵拼接来累积约束。我常用的是一个 constraints 数组从空开始逐步 appendconstraints []; % 电功率平衡 constraints [constraints, P_wt P_pv P_gt P_buy P_load P_p2g P_ccs P_sell P_curt];风电光伏的实际出力不能超过预测值削减量非负这是弃风弃光的来源constraints [constraints, 0 P_wt P_wt_pre, 0 P_pv P_pv_pre, P_curt (P_wt_pre - P_wt) (P_pv_pre - P_pv)];燃气轮机的出力和燃料热量之间建立关系。假设机组效率 η_gt燃料总热量 Q_gt P_gt / η_gt。然后掺氢比约束eta_gt 0.45; H_MAX 0.15; Q_gt P_gt / eta_gt; Q_ng (1 - h_ratio) .* Q_gt; Q_h2 h_ratio .* Q_gt; constraints [constraints, 0 P_gt P_gt_max, 0 h_ratio H_MAX];P2G 产氢eta_p2g 0.65; F_h2 eta_p2g * P_p2g; constraints [constraints, 0 P_p2g P_p2g_max, F_h2 F_h2_gt F_h2_st F_h2_meth];CCS 捕集与耗电beta_ccs 0.25; % MWh/tCO2 eta_ccs 0.90; % 捕集率 constraints [constraints, E_capture eta_ccs * E_emit, P_ccs beta_ccs * E_capture, 0 E_capture E_capture_max];碳排放量按天然气消耗量计算EF_ng 0.2; % tCO2/MWh天然气低位热值口径 constraints [constraints, E_emit EF_ng * Q_ng];储氢罐S_h2_min 0; S_h2_max 80; constraints [constraints, S_h2(1) S_h2_init, S_h2(T1) S_h2_init, 0 S_h2 S_h2_max, S_h2(2:T1) S_h2(1:T) eta_st * F_h2_st - F_h2_gt];注意最后一行储氢动态可以把 F_h2_gt 看成从储氢罐流向燃气轮机的流量但如果你希望“P2G 直接产氢即产即用不需要经过储氢罐”就要把平衡式拆开。我的做法是F_h2_gt 一部分直接来自电解槽一部分来自储氢罐二者通过另一个分流变量控制。简化版里可以直接按“罐进罐出”建模做机理验证够用。4.4 求解与结果提取约束搭好、目标函数写好后调用求解器的代码如下ops sdpsettings(solver, gurobi, verbose, 2, showprogress, 1); sol optimize(constraints, Objective, ops); if sol.problem 0 disp(求解成功); else disp([求解失败: , sol.info]); end P_wt_opt value(P_wt); P_gt_opt value(P_gt); h_ratio_opt value(h_ratio); E_capture_opt value(E_capture); C_carb_opt value(C_carb);5. 典型日算例设计与结果分析算例设计直接决定论文或项目能不能讲出自洽的故事。我建议至少设置 4 个场景来对比这样 P2G-CCS 耦合和掺氢的价值能被单独拆出来看。5.1 场景设置场景一是基准场景不带 P2G、CCS、掺氢燃气轮机纯烧天然气弃风弃光只能靠储能和向上级电网售电来缓解。场景二在基准上加入 P2G 和 CCS但不掺氢用来考察“电-气-碳耦合”单独作用的效果。场景三加入燃气掺氢但碳价是固定值而不是阶梯值用来对比阶梯碳价和固定碳价的差别。场景四用完整模型阶梯碳价 P2G CCS 掺氢全上。这种对比矩阵可以把每个模块的边际贡献都拆出来而且在论文里也很容易画成表格或者柱状图。5.2 算例参数建议下面是我常用的一组参数供参考设备/参数数值风电场容量120 MW光伏电站容量80 MW燃气轮机容量60 MW燃气轮机效率45%P2G 电解槽容量50 MWP2G 制氢效率65%储氢罐容量80 MWh能量口径掺氢比上限15%热值比CCS 捕集率90%CCS 耗电系数0.25 MWh/tCO2甲烷化效率70%免费配额强度0.65 t/MWh阶梯区间30 t / 60 t阶梯碳价80 / 120 / 160 元/t注意这里的所有能量单位都统一到 MWh碳单位统一到吨。如果不统一后面的约束很容易出现数量级错误。我自己在调试时吃过亏储氢罐和甲烷化用不同的能量口径导致氢平衡怎么调都不满足最后发现是单位没换算。5.3 从模型机理看的典型结论我按照这套模型跑下来的经验有几个比较稳定的趋势你可以拿自己的算例验证。第一阶梯碳价下燃气轮机的高出力时段会被明显压缩。尤其是傍晚负荷高峰和电价高峰时段燃气轮机本来应该满发但因为碳排放缺口已经逼近第二档甚至第三档继续高发会触发高碳价调度会把一部分负荷转向上级电网购电。这个转变在场景一里完全看不到因为固定碳价不会产生这种累积惩罚效应。第二掺氢比例在负荷高峰时段倾向于走到上限。原因很直接高峰时段燃气轮机出力大天然气消耗量大排放量也大提高掺氢比例可以压低天然气消耗从而压低碳排放。这个收益在阶梯碳价下比固定碳价更明显因为减少的排放处在高价区间的边际价值更高。第三P2G 的运行策略受“电-氢-气”价差驱动。凌晨风电大发、电价低谷时P2G 出力高制出来的氢一部分进储氢罐一部分用于甲烷化白天电价高峰时段储氢罐放氢和天然气掺混发电甲烷化也优先在碳价压力大的时段运行。如果气价和电价没有足够大的峰谷差P2G 的收益率会明显下降这在参数敏感性分析里能看得很清楚。第四CCS 捕集到 CO2 能否被完全利用取决于甲烷化反应器的容量和氢气供应量。算例里经常出现的情况是CCS 在燃气轮机满发时捕集了大量 CO2但因为 P2G 产氢量不够或者反应器容量限制CO2 用不完只能储存在 CO2 储罐或封存。这时如果模型里没有 CO2 储存约束就等效于“免费排放”结果会过于乐观。6. 调试经验与常见坑最后聊聊我在调试此类模型时踩过的一些坑都是真实花费过时间的地方希望能帮你少走弯路。6.1 掺氢口径不一致导致热值错乱这是最隐蔽的坑。工程上燃气掺氢常用体积比比如掺氢 15%意思是混合燃气中氢气体积占 15%。但相同体积下天然气的热值比氢气高不少直接用体积比乘以燃料热值会低估天然气贡献。我的建议是在模型里用热值掺氢比换算公式是β_energy β_vol · LHV_H2 / (β_vol · LHV_H2 (1 - β_vol) · LHV_NG)其中 LHV_NG 约 36 MJ/m³LHV_H2 约 10.8 MJ/m³。做参数标定时先用这个公式把体积比换算成热值比再代入模型。否则碳排放强度、天然气耗量全都会偏离真实值。6.2 CCS 耗电项被漏掉很多初版模型容易忘记把 P_ccs 放进电功率平衡方程。漏掉之后的结果是CCS 仿佛在不耗电地捕碳捕集量再大也不影响电力调度系统会过度使用 CCS碳捕集成本被严重低估。我习惯在每个设备运行约束旁边同步检查它在哪些平衡方程里出现——P2G 出现在电平衡和氢平衡CCS 出现在电平衡和碳平衡燃气轮机出现在电平衡、气平衡和碳平衡。6.3 CO2 物料守恒约束过于松散如果模型里没有约束“甲烷化消耗的 CO2 不能超过 CCS 捕集量和 CO2 储罐放出量之和”求解器会私自“凭空变出 CO2”让甲烷化白白产气结果特别好看但不物理。这个问题的排查方法很简单把结果里的 CO2 利用量拉出来看看它是否等于捕集量如果不等于说明耦合约束没生效。6.4 储氢罐首末状态约束为了避免调度策略把储氢罐当成免费垃圾桶首末状态约束必须加S_h2(1) S_h2(T1)。不加的话模型会在最后一个时段把所有储氢清空把氢全部烧掉或者卖掉人为拉低总成本。同样的问题也会出现在 CO2 储罐和储能电池上建议所有有存储环节的统一加首末平衡约束。6.5 Yalmip 求解失败时的调试顺序遇到 infeasible 时不要慌按下面的顺序排查通常能定位问题。先用 sdpsettings(verbose,1) 看求解器输出的是哪类错误。如果是“infeasible”把目标函数临时设成常数 0先求可行解如果还是不可行说明约束本身就矛盾。然后逐条“注释”约束同时固定变量取值看哪条约束和其余约束冲突。最容易出问题的通常是电平衡和储氢动态方程尤其是下标维度对不上时Yalmip 不会立刻报错但会生成一个形状不对的约束矩阵。另外Yalmip 的 sdpvar 定义中一个常见错误是忘记 full导致变量变成对称矩阵。写代码时建议用 assert 检查每个变量的 sizeassert(isequal(size(P_wt_pre), [1, T]), P_wt_pre dimension mismatch);检查一遍所有预测数据的维度做数据驱动建模时能省很多事。我在实际项目里还有一个习惯把整个调度模型封装成一个 function输入是风光预测、电价、气价、碳价、设备参数输出是各设备的出力计划、总成本和碳排放量。这样跑敏感性分析和场景对比时不需要反复复制代码只需要在循环里改参数调用函数即可效率和可维护性都会高很多。function [result] vpp_schedule(P_wt_pre, P_pv_pre, P_load, price_buy, price_NG, cparam) % 该函数返回调度结果结构体 result包含各出力序列和成本 ... end参数敏感性分析是这个方向绕不开的一步也是论文里最出彩的部分。我通常会对碳价曲线、掺氢上限、CCS 耗电系数三个参数分别做多组扫描观察总成本、碳排放强度、P2G 利用小时数三个指标的变化这些结果比单独的算例更有说服力。跑完模型再回头看这套系统的核心其实不是某个单一设备而是碳价格信号如何同时影响电力调度、燃料配比和碳捕集水平这三层决策。把阶梯碳价这个关键信号建模清楚了P2G、CCS 和掺氢的调度策略就都跟着顺理成章了。