最近在复现虚拟电厂优化调度方向的EI论文拿到这个题目时我第一反应不是急着找代码而是先把P2G、CCS、燃气掺氢和阶梯碳交易这条线在脑子里过了一遍。这四个东西任何单独一个拿出来都能撑起一篇调度优化文章现在被“虚拟电厂”这个框架串在一起本质上考的已经不是“会不会列约束”而是能不能把一个电-气-碳强耦合的系统用一组线性约束完整表达出来再扔给求解器算出一个有物理意义的结果。这篇记录的是我完整复现这套模型的过程物理链条怎么理解、阶梯碳成本怎么线性化、目标函数和约束怎么写、MatlabYalmip代码骨架怎么搭以及我在调试中踩过的几个坑。内容面向正在做VPP调度、碳交易机制建模、P2G-CCS耦合方向复现或毕设的同学也适合刚入门优化调度的研究生看完至少能知道这类代码的“最小可运行版本”长什么样。1. 先看懂P2G-CCS耦合和燃气掺氢的物理链条1.1 三个单元单独看分别做什么P2G全称Power to Gas本质是“电转气”。主流技术路线是电解水制氢产物是氢气和氧气如果进一步和CO2反应还可以通过甲烷化把氢气合成为甲烷。这一步的化学方程式很关键电解水2H2O → 2H2 O2甲烷化CO2 4H2 → CH4 2H2O从调度角度看P2G最吸引人的属性是“耗电灵活”。传统负荷基本不可控但P2G可以跟随风电光伏的波动在弃风弃光时段增大用电把本来要扔掉的新能源电量转化为氢能或天然气储存起来。也就是说P2G既是负荷又是能源转换枢纽。CCS是碳捕集与封存/利用主要目标是从燃气轮机排出的烟气里把CO2“抓”回来。这里必须强调一个容易被忽略的点CCS本身要耗电。捕集每吨CO2大约消耗0.3~0.6 MWh电能这部分电在虚拟电厂的功率平衡里必须单独列一项否则算出来的系统“净收益”全是虚假的。CCS捕集到的CO2有两条去向一是封存二是作为原料送入P2G的甲烷化单元。标题里的“P2G-CCS耦合”用的就是第二条路。燃气掺氢则相对直观就是把一定比例的氢气混入天然气送入燃气轮机燃烧。氢气燃烧不产生CO2所以掺氢能直接降低燃机碳排放强度。但问题是燃气轮机不是随便就能高比例掺氢的火焰速度、回火风险、NOx排放、材料氢脆都会限制掺氢比例工程上体积掺氢比例通常在10%~30%之间。折算到热值比例会更低因为氢气低位热值约为10.8 MJ/Nm3天然气约为36 MJ/Nm3体积掺氢20%时热量比例大约只有7%。所以建模时一定要先确认掺氢比例究竟是“体积比”还是“热值比”很多复现文章含糊其辞最后对不上结果。1.2 耦合之后物理上形成了一条“碳循环”把三个单元放在同一个虚拟电厂里真正的难点不是单个设备约束而是它们之间的物料耦合关系。P2G的甲烷化需要CO2如果没有CCS就要外购CO2成本和碳账都不划算。CCS从燃机烟气里捕集CO2刚好可以作为甲烷化的原料。而P2G产生的氢气一部分可以送到燃气轮机掺烧一部分进入甲烷化反应器合成甲烷合成的甲烷再进入燃气轮机当燃料。这样一来虚拟电厂内部形成了一条“可再生能源发电 → 电解水制氢 → 氢气掺烧 → 燃气轮机发电 → CO2捕集 → 甲烷化 → 再发电”的循环链。这条循环链在数学上会形成几个硬耦合约束送到甲烷化的CO2不能超过CCS捕集到的CO2总量P2G产出的氢气要在“直接掺烧”和“甲烷化”之间分配甲烷化产生的甲烷要进入燃气轮机的天然气消耗平衡CCS的耗电会增加电负荷P2G的耗电也会增加电负荷这两部分都要进入电功率平衡。为什么要强调这些耦合关系因为很多复现版本只是把P2G、CCS和掺氢三套约束“拼”在一起却没有让它们之间产生物料流动。那样的模型虽然也能跑但完全体现不出标题里“耦合”二字的含义结果自然也不对。我复现时是先把下面四条物质流写在纸上再开始写代码的电风电/光伏/燃气轮机/购电 → 负荷/P2G/CCS气甲烷外购天然气 P2G甲烷化产气 → 燃气轮机消耗氢P2G产氢 → 直接掺烧 / 甲烷化 / 储氢罐碳外购天然气燃烧排放CO2 → CCS捕集 → 甲烷化原料 / 封存。2. 阶梯碳交易模型分段碳价不是罚函数是报价曲线2.1 免费配额和实际排放的口径怎么定碳交易模型的第一步是明确“配额”和“实际排放”分别怎么算。我复现时采用最常见EI论文口径免费碳排放配额按燃气轮机发电量乘以配额强度给出也就是E_quota quota_intensity × P_gt这里有朋友会问虚拟电厂里还包括风电和光伏为什么不给新能源配额原因很简单新能源本身不产生碳排放碳配额机制主要约束化石能源机组所以按燃气轮机出力分配配额是合理的。实际碳排放的计算则要小心。如果不考虑CCS实际排放就是燃气轮机消耗的外购天然气量乘以排放因子。但考虑到CCS捕集的CO2有一部分被送回P2G甲烷化相当于这部分碳被捕集后固定到了甲烷里不进入当期大气排放所以实际排放表达为E_real sum( F_CH4_gt × EF_gas ) - sum( CO2_to_meth )其中F_CH4_gt是燃气轮机消耗的总天然气燃料能量EF_gas是天然气排放因子CO2_to_meth是每时段送往甲烷化单元的CO2量。这里需要确保E_real非负也就是扣除的CO2不能超过总排放否则模型会“凭空造出负碳排放”结果非常荒谬这也是我后面调试时踩过的坑之一。2.2 阶梯碳价的分段线性化方法阶梯碳交易机制的核心特征是碳排放超额量越大超额部分适用的碳价越高。常见的三阶梯模型可以写成超额排放量区间碳价[0, E1]p1[E1, E2]p2[E2, ∞)p3每个区间碳价满足p1 p2 p3。比如E1取100 tE2取200 tp1取20元/tp2取40元/tp3取80元/t意思是前100吨超额按20元/吨计算第101到200吨按40元/吨200吨以上的部分按80元/吨。这个分段函数最直接的写法是C_carbon p1 × min(E_extra, E1) p2 × max(0, min(E_extra, E2) - E1) p3 × max(0, E_extra - E2)但max/min这种形式没法直接扔给线性规划求解器。更通用的做法是引入分段变量把超额排放量E_extra拆成三段E1_seg E2_seg E3_seg E_extra 0 ≤ E1_seg ≤ E1 0 ≤ E2_seg ≤ E2 - E1 E3_seg ≥ 0碳交易成本变为C_carbon p1 × E1_seg p2 × E2_seg p3 × E3_seg这里有一个很多人没想明白的点为什么不需要额外加0-1变量来保证“先填满第一段再进入第二段”因为在p1 p2 p3且目标函数要求总成本最小的情况下求解器会自动优先使用低价段。也就是说最优解里E1_seg一定会先被填满E2_seg才会开始增加。如果碳价序列不是严格递增的或者目标函数里存在其他非线性项破坏这种单调性才需要用大M法引入0-1变量做严格分段。我们在复现中先按价格递增的常规情况处理模型规模更小求解速度也更快。还有一个细节超额排放量E_extra本身是max(E_real - E_quota, 0)。为了不引入max非线性项我采用两个约束实现E_extra ≥ E_real - E_quota E_extra ≥ 0由于碳交易成本在目标函数中是需要最小化的正数求解器会自动把E_extra压到max(E_real - E_quota, 0)。这比显式写max要稳妥得多。3. 把调度问题写成数学优化目标函数与约束方程3.1 目标函数总成本最小化虚拟电厂作为运营商目标是让整个系统在满足负荷需求的前提下总运行成本最低。我复现时的目标函数包含六项min C_total C_buy C_gas C_p2g C_ccs C_carbon C_curtail拆开解释C_buy是向上级电网购电的费用形式为分时电价乘以购电功率再乘以时段长度C_gas是外购天然气费用只计算真正外购的天然气量P2G甲烷化产出的甲烷不重复计费C_p2g是P2G设备运行维护成本按耗电量线性计算C_ccs是CCS运行成本按CO2捕集量线性计算包含溶剂损耗、压缩能耗、设备维护等折算C_carbon是上一节算出的阶梯碳交易成本C_curtail是弃风弃光惩罚用来驱动模型优先消纳新能源。弃风弃光惩罚项尤其重要。如果没有这项模型很可能为了降低运行成本而主动放弃一部分新能源这不符合“优先消纳新能源”的基本前提。惩罚系数一般要设置得比购电价和碳成本都高才能保证模型在能消纳时不弃。3.2 电功率平衡所有电源侧等于所有负荷侧对每个时段t电功率平衡方程为P_use_wind(t) P_use_pv(t) P_gt(t) P_buy(t) P_load(t) P_p2g(t) P_ccs(t)其中P_use_wind和P_use_pv是可再生能源实际出力它们还要满足上限约束0 ≤ P_use_wind(t) ≤ P_wind_avail(t) 0 ≤ P_use_pv(t) ≤ P_pv_avail(t)未被使用的部分就是弃风弃光量记为P_curtail(t) P_wind_avail(t) - P_use_wind(t) P_pv_avail(t) - P_use_pv(t)这里有个建模习惯问题不要在电平衡里直接写“P_wind P_pv ... ... P_curtail”那样会把弃风弃光错误地变成“负荷”。正确的做法是让可再生出力作为变量弃风弃光作为可再生预测值与实际使用值之差。燃气轮机的电出力约束包括上下限和爬坡约束P_gt_min × u_gt(t) ≤ P_gt(t) ≤ P_gt_max × u_gt(t)ramp ≤ P_gt(t) - P_gt(t-1) ≤ ramp其中u_gt是表示启停状态的0-1变量P_gt_min和P_gt_max是机组最小/最大技术出力ramp是爬坡速率上限。加入启停变量之后这个模型就正式变成MILP混合整数线性规划需要用支持整数变量的求解器比如Gurobi或CPLEX。3.3 气-氢-碳耦合约束线性化是核心气侧平衡要处理“甲烷从哪里来、到哪里去”F_CH4_ext(t) CH4_p2g(t) F_CH4_gt(t)其中F_CH4_ext是外购天然气燃料功率CH4_p2g是P2G甲烷化产生的甲烷功率F_CH4_gt是燃气轮机消耗的总天然气功率。这个约束把“外购天然气”和“P2G产甲烷”统一到一个渠道里。燃气轮机的出力与燃料消耗关系为P_gt(t) eta_gt × ( F_CH4_gt(t) F_H2_fuel(t) )这里为了简化单位我把F_CH4_gt和F_H2_fuel都定义成“输入燃料功率MW”eta_gt是燃气轮机效率。掺氢比例按热量比定义为alpha(t) F_H2_fuel(t) / ( F_CH4_gt(t) F_H2_fuel(t) )这个表达式因为有除法直接写进约束里会产生非线性项。解决办法是改写成线性不等式F_H2_fuel(t) ≤ alpha_max / (1 - alpha_max) × F_CH4_gt(t)如果alpha_max取0.15那么alpha_max/(1-alpha_max)就是0.1765约束变成“氢气燃料功率不能超过天然气燃料功率的17.65%”非常清爽。P2G制氢环节在大时间尺度下可以简化成线性效率模型H2_prod(t) eta_p2g × P_p2g(t) × dt其中H2_prod是该时段产氢能量MWhP_p2g是P2G耗电功率MWdt是时段长度小时eta_p2g是电解槽的电转氢效率。氢气产量需要在“直接掺烧F_H2_fuel”“甲烷化H2_meth”和“储氢罐充注”之间分配Soc_h2(t1) Soc_h2(t) H2_prod(t) - H2_fuel(t) - H2_meth(t)储氢罐容量约束和首末周期约束为0 ≤ Soc_h2(t) ≤ H2_cap Soc_h2(1) Soc_h2(T1)甲烷化单元消耗氢气和CO2产出甲烷。我复现时采用的折算关系是1 MWh氢气按低位热值33.3 MWh/t计算完全甲烷化大约需要5.5 t CO2因此CO2_to_meth(t) K_CO2_H2 × H2_meth(t)其中K_CO2_H2 5.5 / 33.3约等于0.165 tCO2/MWh氢气。同时考虑甲烷化效率eta_meth产甲烷功率为CH4_p2g(t) eta_meth × H2_meth(t)CCS捕集量上限由燃气轮机燃烧天然气的总排放量决定并受捕集率和捕集容量约束CO2_cap(t) ≤ capture_rate × F_CH4_gt(t) × EF_gas × dt CO2_cap(t) ≤ CCS_capacityCCS耗电为P_ccs(t) k_ccs_elec × CO2_cap(t)最后送进甲烷化的CO2不能超过CCS实际捕集的CO2CO2_to_meth(t) ≤ CO2_cap(t)这套约束写下来模型就真正成了“电-气-碳”三个物理量相互咬合的整体。任何一个环节的数据口径不一致最优解都会偏离物理实际。4. MatlabYalmip代码实现从参数声明到求解器落地4.1 算例参数怎么设复现论文类题目一定不要一上来就写代码先把算例参数表列清楚。我采用24时段调度时间间隔dt1小时系统包含一台100 MW燃气轮机、50 MW风电、30 MW光伏、50 MW P2G设备、40 t/h的CCS装置和30 MWh储氢罐。关键参数如下参数数值时段数 T24燃气轮机容量100 MW燃气轮机效率0.45掺氢比例上限热值比15%P2G容量50 MW电解槽效率0.70甲烷化效率0.83CCS捕集率0.90CCS单位电耗0.40 MWh/t天然气排放因子0.20 tCO2/MWh免费配额强度0.05 t/MWh碳价分段 [p1,p2,p3][20, 40, 80] 元/t储氢罐容量30 MWh注意这套参数不是从天上掉下来的而是参考了多篇同类EI论文的典型设定再结合自己算例调整得到的。复现时千万不要原封不动照搬任何一篇文章的参数一定要先理解每个参数的含义再根据自己系统的容量规模换算。4.2 核心代码骨架下面是我调通后的MatlabYalmip代码骨架去掉了数据加载部分保留了变量声明、约束构建、目标函数和求解设置。代码里单位统一为“功率MW、能量MWh、CO2质量t”。%% 基本数据 T 24; dt 1; % 以示例向量代替实际预测数据实际运行时替换为真实曲线即可 P_wind_avail [ ... 1x24 风电预测 ... ]; P_pv_avail [ ... 1x24 光伏预测 ... ]; P_load [ ... 1x24 负荷预测 ... ]; price_buy [ ... 1x24 分时购电价 ... ]; price_gas 350; % 天然气价格元/MWh eta_gt 0.45; % 燃气轮机效率 P_gt_max 100; P_gt_min 20; ramp 60; eta_p2g 0.70; P_p2g_max 50; eta_meth 0.83; K_CO2_H2 5.5 / 33.3; % tCO2/MWh氢气 EF_gas 0.20; % tCO2/MWh天然气 capture_rate 0.90; CCS_capacity 40; % t/h k_ccs_elec 0.40; % MWh/t alpha_max 0.15; % 掺氢比例上限热值比 H2_cap 30; % MWh %% 碳交易参数 quota_intensity 0.05; % tCO2/MWh E1 100; % 阶梯分段长度 E2 200; p1 20; p2 40; p3 80; c_cur 500; % 弃风弃光惩罚元/MWh %% 变量声明 P_use_wind sdpvar(1, T); P_use_pv sdpvar(1, T); P_gt sdpvar(1, T); u_gt binvar(1, T); % 燃机启停0-1变量 P_buy sdpvar(1, T); P_p2g sdpvar(1, T); P_ccs sdpvar(1, T); F_CH4_gt sdpvar(1, T); % 燃机消耗的总天然气燃料功率MW F_CH4_ext sdpvar(1, T); % 外购天然气燃料功率MW CH4_p2g sdpvar(1, T); % P2G甲烷化产甲烷功率MW F_H2_fuel sdpvar(1, T); % 掺烧氢气燃料功率MW H2_prod sdpvar(1, T); % 产氢能量MWh H2_meth sdpvar(1, T); % 甲烷化消耗氢能量MWh Soc_h2 sdpvar(1, T1); % 储氢罐状态MWh CO2_cap sdpvar(1, T); % CCS捕集CO2量t CO2_to_meth sdpvar(1, T); % 送往甲烷化的CO2量t E_extra sdpvar(1); E1_seg sdpvar(1); E2_seg sdpvar(1); E3_seg sdpvar(1); C []; % 约束集合 %% 电功率平衡 C [C, P_use_wind P_use_pv P_gt P_buy P_load P_p2g P_ccs]; C [C, 0 P_use_wind P_wind_avail]; C [C, 0 P_use_pv P_pv_avail]; C [C, 0 P_buy 200]; % 购电上限可根据外部电网容量调整 %% 燃气轮机约束 C [C, P_gt_min * u_gt P_gt P_gt_max * u_gt]; C [C, P_gt eta_gt * (F_CH4_gt F_H2_fuel)]; C [C, F_H2_fuel (alpha_max / (1 - alpha_max)) * F_CH4_gt]; for t 2:T C [C, -ramp P_gt(t) - P_gt(t-1) ramp]; end %% 天然气平衡 C [C, F_CH4_ext CH4_p2g F_CH4_gt]; C [C, F_CH4_ext 0]; %% P2G与储氢约束 C [C, 0 P_p2g P_p2g_max]; C [C, H2_prod eta_p2g * P_p2g * dt]; C [C, Soc_h2(1) 0.5 * H2_cap]; C [C, Soc_h2(T1) Soc_h2(1)]; for t 1:T C [C, Soc_h2(t1) Soc_h2(t) H2_prod(t) - F_H2_fuel(t) * dt - H2_meth(t)]; C [C, 0 Soc_h2(t1) H2_cap]; C [C, 0 H2_meth(t)]; end %% 甲烷化与CCS耦合 C [C, CH4_p2g eta_meth * H2_meth]; C [C, CO2_to_meth K_CO2_H2 * H2_meth]; C [C, CO2_to_meth CO2_cap]; C [C, CO2_cap capture_rate * F_CH4_gt .* EF_gas * dt]; C [C, CO2_cap CCS_capacity * dt]; C [C, P_ccs k_ccs_elec * CO2_cap]; %% 碳排放与阶梯碳交易 E_quota quota_intensity * sum(P_gt .* dt); E_real sum(F_CH4_gt .* EF_gas .* dt) - sum(CO2_to_meth); C [C, E_real 0]; C [C, E_extra E_real - E_quota]; C [C, E_extra 0]; C [C, E1_seg E2_seg E3_seg E_extra]; C [C, 0 E1_seg E1]; C [C, 0 E2_seg E2 - E1]; C [C, E3_seg 0]; C_carbon p1 * E1_seg p2 * E2_seg p3 * E3_seg; %% 目标函数 C_curtail c_cur * ( sum(P_wind_avail - P_use_wind) sum(P_pv_avail - P_use_pv) ) * dt; objective sum(price_buy .* P_buy) * dt ... price_gas * sum(F_CH4_ext) * dt ... 30 * sum(P_p2g) * dt ... % P2G运维成本 25 * sum(CO2_cap) ... % CCS运行成本 C_carbon ... C_curtail; %% 求解 ops sdpsettings(solver, gurobi, verbose, 2, mipgap, 0.001); diagnosis optimize(C, objective, ops); if diagnosis.problem 0 disp(求解成功); P_gt_opt value(P_gt); P_p2g_opt value(P_p2g); CO2_cap_opt value(CO2_cap); C_carbon_opt value(C_carbon); else disp(diagnosis.info); end这段代码虽然离“可直接投稿”还有距离但已经具备完整的MILP骨架。把预测数据换进去再补上结果画图就是一个能跑通的最小复现系统。4.3 结果输出与可视化建议求解完成后我一般会先输出四张图24小时电功率平衡堆叠图观察风电、光伏、燃机、购电如何共同满足负荷P2G耗电、制氢量、掺氢量和甲烷化量的时间曲线判断氢气的分配是否合理燃气轮机出力与CCS捕集量对照图看碳排放和碳捕集的时序关系阶梯碳交易成本构成饼图看哪一段碳价区间贡献了主要成本。可视化的最大价值不是“好看”而是能快速发现物理不合理的解。比如某小时P2G大量产氢但储氢罐SOC曲线没有同步变化那一定是平衡方程写错了再比如CCS捕集量在燃机停机时仍然大于0那一定是CO2_cap的上限约束没写全。5. 复现结果里几个值得细看的调度现象5.1 P2G-CCS耦合把弃风时段变成了“燃料生产时段”在基准场景下最容易观察到的现象是夜间负荷下降、风电出力较大且电价处于谷段时燃气轮机出力减少P2G开始增大耗电把多余风电转化成氢气。此时CCS如果还在运行捕集到的CO2不会白白封存而是送入甲烷化反应器与氢气合成甲烷。这些甲烷会存在气路里或直接补充到天然气管道等到早上负荷高峰时再被燃气轮机烧掉。这其实说明了一个很核心的调度逻辑P2G-CCS耦合不是简单“为了环保而环保”它是在电价低谷时用低成本电能合成高品位燃料把“弃风时段”变成了虚拟电厂自己的“燃料生产时段”。从经济账上看它同时节省了外购天然气费用也减少了弃风惩罚还降低了碳交易成本。5.2 阶梯碳价会同时改变燃机出力和CCS运行策略碳交易机制对调度的引导作用不是“一刀切”地减少燃机出力而是会重构整个系统的最优运行点。碳价提高后我观察到的变化包括燃气轮机在负荷峰段的出力适当下调部分电量由购电替代但燃机不会完全停掉因为P2G需要CO2来合成甲烷而CO2的主要来源就是燃机烟气CCS捕集率会向设备上限靠近因为捕集的CO2在系统中具有“原料价值”和“碳减排价值”双重收益风电商用空间增大弃风率明显降低。这一点对于理解“碳交易为什么能促进新能源消纳”很有帮助。碳价本质上给化石能源发电增加了一个额外成本而这个成本会沿着“燃机出力 → 碳排放 → 碳交易成本 → 总成本”的路径传导到每一度电上最终迫使系统寻找更便宜的替代方案。5.3 掺氢比例不是越高越好存在经济最优区间很多复现者做完基础算例就开始狂拉掺氢比例上限以为掺得越多越好。实际上在P2G产氢量受制于电解槽容量和弃风资源的条件下盲目调大掺氢上限通常会得到“掺氢比例约束根本没被触碰到”的结果。从成本结构看掺氢需要消耗P2G产出的氢气而氢气如果用于甲烷化可以替代外购天然气如果用于直接掺烧则替代的也是天然气但减排效果更直接。两者存在竞争关系。随着掺氢比例上限提高模型会在这两种氢气去向之间做权衡。当掺氢比例较低时提高上限能有效降低天然气成本和碳交易成本但当比例高到一定程度后继续提高只会把氢气从甲烷化那边抢过来反而可能减少甲烷产量导致外购天然气增加总成本不降反升。所以“最优掺氢区间”是系统运行的经济结果不是越大越好。5.4 对比场景设计一组能证明“耦合与碳交易都起了作用”的算例复现类文章最忌讳只有一个全局最优解。我建议至少跑下面四个场景场景是否启用P2G-CCS耦合是否启用阶梯碳交易表现结论S1否否基准系统弃风率最高S2是否消纳改善但低碳价值未被价格化S3否是碳成本高但无法回收CO2S4是是总成本与碳排放综合最优S2和S3的对比尤其重要。只看S2会以为“耦合本身已经很好了”只看S3又会以为“碳交易只是增加负担”。只有S4才能展现出“耦合为碳交易提供减碳手段碳交易为耦合提供经济回报”这个正反馈机制。这不只是论文里的“敏感性分析”更是理解标题里几个关键词之间关系的钥匙。6. 复现时踩过的五个坑与排查思路6.1 量纲不统一结果凭空多出两个数量级第一次跑出结果时我心里隐约觉得不对劲燃气轮机一天的发电量只有几百MWh但CO2捕集量居然算出了上万吨。排查到最后根源是我把“燃料功率MW”和“燃料能量MWh”混用了。F_CH4_gt这个变量我一边在燃气轮机出力方程里当成MW用一边又在碳排放公式里乘了EF_gas后当成MWh用量纲差了一倍结果自然乱套。排查量纲问题的方法很简单从头到尾检查每个约束的左右单位是否一致。功率项统一用MW能量项统一用MWh质量项统一用t。时间粒度为1小时时MW乘以1小时就是MWht/h乘以1小时就是t。把功夫花在单位表上永远不亏。6.2 掺氢比例约束写成了除法模型直接变成非线性最初我图省事把掺氢比例约束直接写成alpha(t) F_H2_fuel(t) / (F_CH4_gt(t) F_H2_fuel(t));Yalmip立刻提示Nonlinear constraints求解器也换成了非线性求解器不仅慢还不一定能找到全局最优。后来才意识到这个约束可以巧妙地线性化成F_H2_fuel(t) ≤ (alpha_max / (1 - alpha_max)) × F_CH4_gt(t)这其实是很多新手容易忽略的地方。看到除法和乘法出现在两个变量之间就应当马上思考能不能通过移项改写为线性形式。如果能尽量用线性约束替代否则求解效率和稳定性都会变差。6.3 碳交易分段变量上界漏写出现负碳价或“免费排放”还有一次碳交易成本算出来是负数导致系统拼命增加碳排放来“赚钱”。查完之后发现E1_seg和E2_seg的上界约束漏写了求解器把E_extra拆成了“第一段负数、第二段正数”的形式再用不同碳价一加就凭空造出了负成本。分段碳成本模型必须写全三个约束0 ≤ E1_seg ≤ E1 0 ≤ E2_seg ≤ E2 - E1 E3_seg ≥ 0一个都不能少。上界是物理分段边界下界保证碳价权重不会出现负值。这里没有捷径检查时务必逐行核对。6.4 储氢罐SOC没有闭环约束系统开始“凭空产氢”储氢罐的初始SOC和末端SOC必须要固定或者至少要让整个周期的充放电量满足闭环约束否则模型会把储氢罐当成一个无限能量的“魔法罐”。我最初只加了SOC上下限约束忘记加Soc_h2(T1) Soc_h2(1)结果最优解里储氢罐每天从0充到30 MWh又全部放空成本被严重低估看起来像是“每天都在白赚30 MWh氢气”。修复方式很简单在约束里加上Soc_h2(1) Soc_h2(T1)如果是调度周期重复执行的情况还可以把初始SOC作为优化变量但至少要在两个周期之间保持能量守恒不能出现能量凭空产生。6.5 调试策略先锁定耦合项再逐层放开最后分享一个非常实用的调试顺序。拿到这套模型后不要一上来就把P2G、CCS、掺氢、碳交易全开着跑。我的做法是分四步走先固定P2G和CCS出力为常数只优化燃气轮机、购电和新能源消纳验证电功率平衡与机组约束是否正确再放开P2G让制氢量成为变量检查氢气平衡和储氢罐SOC是否闭合再放开CCS和甲烷化让CO2捕集量与CO2_to_meth咬合观察碳平衡最后才把阶梯碳交易和掺氢比例上限全部加入。每一步都对比前后结果如果某一步出现明显不合理比如弃风率异常、某设备出力恒为0错误大概率就在刚放开的那一层里。这个顺序帮我省下了至少一个晚上的调试时间也建议你复现任何耦合系统时照做。在我个人看来这类EI论文复现真正有价值的收获不是最后那张漂亮的调度曲线图而是从头到尾把“物理过程 → 数学模型 → 代码约束 → 求解结果”这条链路走通的过程。只要单位统一、边界条件完整、耦合关系不丢哪怕不追求和原论文一模一样的结果这套模型本身也已经能作为后续研究的可靠基线了。 虚拟电厂优化调度中的P2G-CCS耦合与阶梯碳交易建模复现
最近在复现虚拟电厂优化调度方向的EI论文拿到这个题目时我第一反应不是急着找代码而是先把P2G、CCS、燃气掺氢和阶梯碳交易这条线在脑子里过了一遍。这四个东西任何单独一个拿出来都能撑起一篇调度优化文章现在被“虚拟电厂”这个框架串在一起本质上考的已经不是“会不会列约束”而是能不能把一个电-气-碳强耦合的系统用一组线性约束完整表达出来再扔给求解器算出一个有物理意义的结果。这篇记录的是我完整复现这套模型的过程物理链条怎么理解、阶梯碳成本怎么线性化、目标函数和约束怎么写、MatlabYalmip代码骨架怎么搭以及我在调试中踩过的几个坑。内容面向正在做VPP调度、碳交易机制建模、P2G-CCS耦合方向复现或毕设的同学也适合刚入门优化调度的研究生看完至少能知道这类代码的“最小可运行版本”长什么样。1. 先看懂P2G-CCS耦合和燃气掺氢的物理链条1.1 三个单元单独看分别做什么P2G全称Power to Gas本质是“电转气”。主流技术路线是电解水制氢产物是氢气和氧气如果进一步和CO2反应还可以通过甲烷化把氢气合成为甲烷。这一步的化学方程式很关键电解水2H2O → 2H2 O2甲烷化CO2 4H2 → CH4 2H2O从调度角度看P2G最吸引人的属性是“耗电灵活”。传统负荷基本不可控但P2G可以跟随风电光伏的波动在弃风弃光时段增大用电把本来要扔掉的新能源电量转化为氢能或天然气储存起来。也就是说P2G既是负荷又是能源转换枢纽。CCS是碳捕集与封存/利用主要目标是从燃气轮机排出的烟气里把CO2“抓”回来。这里必须强调一个容易被忽略的点CCS本身要耗电。捕集每吨CO2大约消耗0.3~0.6 MWh电能这部分电在虚拟电厂的功率平衡里必须单独列一项否则算出来的系统“净收益”全是虚假的。CCS捕集到的CO2有两条去向一是封存二是作为原料送入P2G的甲烷化单元。标题里的“P2G-CCS耦合”用的就是第二条路。燃气掺氢则相对直观就是把一定比例的氢气混入天然气送入燃气轮机燃烧。氢气燃烧不产生CO2所以掺氢能直接降低燃机碳排放强度。但问题是燃气轮机不是随便就能高比例掺氢的火焰速度、回火风险、NOx排放、材料氢脆都会限制掺氢比例工程上体积掺氢比例通常在10%~30%之间。折算到热值比例会更低因为氢气低位热值约为10.8 MJ/Nm3天然气约为36 MJ/Nm3体积掺氢20%时热量比例大约只有7%。所以建模时一定要先确认掺氢比例究竟是“体积比”还是“热值比”很多复现文章含糊其辞最后对不上结果。1.2 耦合之后物理上形成了一条“碳循环”把三个单元放在同一个虚拟电厂里真正的难点不是单个设备约束而是它们之间的物料耦合关系。P2G的甲烷化需要CO2如果没有CCS就要外购CO2成本和碳账都不划算。CCS从燃机烟气里捕集CO2刚好可以作为甲烷化的原料。而P2G产生的氢气一部分可以送到燃气轮机掺烧一部分进入甲烷化反应器合成甲烷合成的甲烷再进入燃气轮机当燃料。这样一来虚拟电厂内部形成了一条“可再生能源发电 → 电解水制氢 → 氢气掺烧 → 燃气轮机发电 → CO2捕集 → 甲烷化 → 再发电”的循环链。这条循环链在数学上会形成几个硬耦合约束送到甲烷化的CO2不能超过CCS捕集到的CO2总量P2G产出的氢气要在“直接掺烧”和“甲烷化”之间分配甲烷化产生的甲烷要进入燃气轮机的天然气消耗平衡CCS的耗电会增加电负荷P2G的耗电也会增加电负荷这两部分都要进入电功率平衡。为什么要强调这些耦合关系因为很多复现版本只是把P2G、CCS和掺氢三套约束“拼”在一起却没有让它们之间产生物料流动。那样的模型虽然也能跑但完全体现不出标题里“耦合”二字的含义结果自然也不对。我复现时是先把下面四条物质流写在纸上再开始写代码的电风电/光伏/燃气轮机/购电 → 负荷/P2G/CCS气甲烷外购天然气 P2G甲烷化产气 → 燃气轮机消耗氢P2G产氢 → 直接掺烧 / 甲烷化 / 储氢罐碳外购天然气燃烧排放CO2 → CCS捕集 → 甲烷化原料 / 封存。2. 阶梯碳交易模型分段碳价不是罚函数是报价曲线2.1 免费配额和实际排放的口径怎么定碳交易模型的第一步是明确“配额”和“实际排放”分别怎么算。我复现时采用最常见EI论文口径免费碳排放配额按燃气轮机发电量乘以配额强度给出也就是E_quota quota_intensity × P_gt这里有朋友会问虚拟电厂里还包括风电和光伏为什么不给新能源配额原因很简单新能源本身不产生碳排放碳配额机制主要约束化石能源机组所以按燃气轮机出力分配配额是合理的。实际碳排放的计算则要小心。如果不考虑CCS实际排放就是燃气轮机消耗的外购天然气量乘以排放因子。但考虑到CCS捕集的CO2有一部分被送回P2G甲烷化相当于这部分碳被捕集后固定到了甲烷里不进入当期大气排放所以实际排放表达为E_real sum( F_CH4_gt × EF_gas ) - sum( CO2_to_meth )其中F_CH4_gt是燃气轮机消耗的总天然气燃料能量EF_gas是天然气排放因子CO2_to_meth是每时段送往甲烷化单元的CO2量。这里需要确保E_real非负也就是扣除的CO2不能超过总排放否则模型会“凭空造出负碳排放”结果非常荒谬这也是我后面调试时踩过的坑之一。2.2 阶梯碳价的分段线性化方法阶梯碳交易机制的核心特征是碳排放超额量越大超额部分适用的碳价越高。常见的三阶梯模型可以写成超额排放量区间碳价[0, E1]p1[E1, E2]p2[E2, ∞)p3每个区间碳价满足p1 p2 p3。比如E1取100 tE2取200 tp1取20元/tp2取40元/tp3取80元/t意思是前100吨超额按20元/吨计算第101到200吨按40元/吨200吨以上的部分按80元/吨。这个分段函数最直接的写法是C_carbon p1 × min(E_extra, E1) p2 × max(0, min(E_extra, E2) - E1) p3 × max(0, E_extra - E2)但max/min这种形式没法直接扔给线性规划求解器。更通用的做法是引入分段变量把超额排放量E_extra拆成三段E1_seg E2_seg E3_seg E_extra 0 ≤ E1_seg ≤ E1 0 ≤ E2_seg ≤ E2 - E1 E3_seg ≥ 0碳交易成本变为C_carbon p1 × E1_seg p2 × E2_seg p3 × E3_seg这里有一个很多人没想明白的点为什么不需要额外加0-1变量来保证“先填满第一段再进入第二段”因为在p1 p2 p3且目标函数要求总成本最小的情况下求解器会自动优先使用低价段。也就是说最优解里E1_seg一定会先被填满E2_seg才会开始增加。如果碳价序列不是严格递增的或者目标函数里存在其他非线性项破坏这种单调性才需要用大M法引入0-1变量做严格分段。我们在复现中先按价格递增的常规情况处理模型规模更小求解速度也更快。还有一个细节超额排放量E_extra本身是max(E_real - E_quota, 0)。为了不引入max非线性项我采用两个约束实现E_extra ≥ E_real - E_quota E_extra ≥ 0由于碳交易成本在目标函数中是需要最小化的正数求解器会自动把E_extra压到max(E_real - E_quota, 0)。这比显式写max要稳妥得多。3. 把调度问题写成数学优化目标函数与约束方程3.1 目标函数总成本最小化虚拟电厂作为运营商目标是让整个系统在满足负荷需求的前提下总运行成本最低。我复现时的目标函数包含六项min C_total C_buy C_gas C_p2g C_ccs C_carbon C_curtail拆开解释C_buy是向上级电网购电的费用形式为分时电价乘以购电功率再乘以时段长度C_gas是外购天然气费用只计算真正外购的天然气量P2G甲烷化产出的甲烷不重复计费C_p2g是P2G设备运行维护成本按耗电量线性计算C_ccs是CCS运行成本按CO2捕集量线性计算包含溶剂损耗、压缩能耗、设备维护等折算C_carbon是上一节算出的阶梯碳交易成本C_curtail是弃风弃光惩罚用来驱动模型优先消纳新能源。弃风弃光惩罚项尤其重要。如果没有这项模型很可能为了降低运行成本而主动放弃一部分新能源这不符合“优先消纳新能源”的基本前提。惩罚系数一般要设置得比购电价和碳成本都高才能保证模型在能消纳时不弃。3.2 电功率平衡所有电源侧等于所有负荷侧对每个时段t电功率平衡方程为P_use_wind(t) P_use_pv(t) P_gt(t) P_buy(t) P_load(t) P_p2g(t) P_ccs(t)其中P_use_wind和P_use_pv是可再生能源实际出力它们还要满足上限约束0 ≤ P_use_wind(t) ≤ P_wind_avail(t) 0 ≤ P_use_pv(t) ≤ P_pv_avail(t)未被使用的部分就是弃风弃光量记为P_curtail(t) P_wind_avail(t) - P_use_wind(t) P_pv_avail(t) - P_use_pv(t)这里有个建模习惯问题不要在电平衡里直接写“P_wind P_pv ... ... P_curtail”那样会把弃风弃光错误地变成“负荷”。正确的做法是让可再生出力作为变量弃风弃光作为可再生预测值与实际使用值之差。燃气轮机的电出力约束包括上下限和爬坡约束P_gt_min × u_gt(t) ≤ P_gt(t) ≤ P_gt_max × u_gt(t)ramp ≤ P_gt(t) - P_gt(t-1) ≤ ramp其中u_gt是表示启停状态的0-1变量P_gt_min和P_gt_max是机组最小/最大技术出力ramp是爬坡速率上限。加入启停变量之后这个模型就正式变成MILP混合整数线性规划需要用支持整数变量的求解器比如Gurobi或CPLEX。3.3 气-氢-碳耦合约束线性化是核心气侧平衡要处理“甲烷从哪里来、到哪里去”F_CH4_ext(t) CH4_p2g(t) F_CH4_gt(t)其中F_CH4_ext是外购天然气燃料功率CH4_p2g是P2G甲烷化产生的甲烷功率F_CH4_gt是燃气轮机消耗的总天然气功率。这个约束把“外购天然气”和“P2G产甲烷”统一到一个渠道里。燃气轮机的出力与燃料消耗关系为P_gt(t) eta_gt × ( F_CH4_gt(t) F_H2_fuel(t) )这里为了简化单位我把F_CH4_gt和F_H2_fuel都定义成“输入燃料功率MW”eta_gt是燃气轮机效率。掺氢比例按热量比定义为alpha(t) F_H2_fuel(t) / ( F_CH4_gt(t) F_H2_fuel(t) )这个表达式因为有除法直接写进约束里会产生非线性项。解决办法是改写成线性不等式F_H2_fuel(t) ≤ alpha_max / (1 - alpha_max) × F_CH4_gt(t)如果alpha_max取0.15那么alpha_max/(1-alpha_max)就是0.1765约束变成“氢气燃料功率不能超过天然气燃料功率的17.65%”非常清爽。P2G制氢环节在大时间尺度下可以简化成线性效率模型H2_prod(t) eta_p2g × P_p2g(t) × dt其中H2_prod是该时段产氢能量MWhP_p2g是P2G耗电功率MWdt是时段长度小时eta_p2g是电解槽的电转氢效率。氢气产量需要在“直接掺烧F_H2_fuel”“甲烷化H2_meth”和“储氢罐充注”之间分配Soc_h2(t1) Soc_h2(t) H2_prod(t) - H2_fuel(t) - H2_meth(t)储氢罐容量约束和首末周期约束为0 ≤ Soc_h2(t) ≤ H2_cap Soc_h2(1) Soc_h2(T1)甲烷化单元消耗氢气和CO2产出甲烷。我复现时采用的折算关系是1 MWh氢气按低位热值33.3 MWh/t计算完全甲烷化大约需要5.5 t CO2因此CO2_to_meth(t) K_CO2_H2 × H2_meth(t)其中K_CO2_H2 5.5 / 33.3约等于0.165 tCO2/MWh氢气。同时考虑甲烷化效率eta_meth产甲烷功率为CH4_p2g(t) eta_meth × H2_meth(t)CCS捕集量上限由燃气轮机燃烧天然气的总排放量决定并受捕集率和捕集容量约束CO2_cap(t) ≤ capture_rate × F_CH4_gt(t) × EF_gas × dt CO2_cap(t) ≤ CCS_capacityCCS耗电为P_ccs(t) k_ccs_elec × CO2_cap(t)最后送进甲烷化的CO2不能超过CCS实际捕集的CO2CO2_to_meth(t) ≤ CO2_cap(t)这套约束写下来模型就真正成了“电-气-碳”三个物理量相互咬合的整体。任何一个环节的数据口径不一致最优解都会偏离物理实际。4. MatlabYalmip代码实现从参数声明到求解器落地4.1 算例参数怎么设复现论文类题目一定不要一上来就写代码先把算例参数表列清楚。我采用24时段调度时间间隔dt1小时系统包含一台100 MW燃气轮机、50 MW风电、30 MW光伏、50 MW P2G设备、40 t/h的CCS装置和30 MWh储氢罐。关键参数如下参数数值时段数 T24燃气轮机容量100 MW燃气轮机效率0.45掺氢比例上限热值比15%P2G容量50 MW电解槽效率0.70甲烷化效率0.83CCS捕集率0.90CCS单位电耗0.40 MWh/t天然气排放因子0.20 tCO2/MWh免费配额强度0.05 t/MWh碳价分段 [p1,p2,p3][20, 40, 80] 元/t储氢罐容量30 MWh注意这套参数不是从天上掉下来的而是参考了多篇同类EI论文的典型设定再结合自己算例调整得到的。复现时千万不要原封不动照搬任何一篇文章的参数一定要先理解每个参数的含义再根据自己系统的容量规模换算。4.2 核心代码骨架下面是我调通后的MatlabYalmip代码骨架去掉了数据加载部分保留了变量声明、约束构建、目标函数和求解设置。代码里单位统一为“功率MW、能量MWh、CO2质量t”。%% 基本数据 T 24; dt 1; % 以示例向量代替实际预测数据实际运行时替换为真实曲线即可 P_wind_avail [ ... 1x24 风电预测 ... ]; P_pv_avail [ ... 1x24 光伏预测 ... ]; P_load [ ... 1x24 负荷预测 ... ]; price_buy [ ... 1x24 分时购电价 ... ]; price_gas 350; % 天然气价格元/MWh eta_gt 0.45; % 燃气轮机效率 P_gt_max 100; P_gt_min 20; ramp 60; eta_p2g 0.70; P_p2g_max 50; eta_meth 0.83; K_CO2_H2 5.5 / 33.3; % tCO2/MWh氢气 EF_gas 0.20; % tCO2/MWh天然气 capture_rate 0.90; CCS_capacity 40; % t/h k_ccs_elec 0.40; % MWh/t alpha_max 0.15; % 掺氢比例上限热值比 H2_cap 30; % MWh %% 碳交易参数 quota_intensity 0.05; % tCO2/MWh E1 100; % 阶梯分段长度 E2 200; p1 20; p2 40; p3 80; c_cur 500; % 弃风弃光惩罚元/MWh %% 变量声明 P_use_wind sdpvar(1, T); P_use_pv sdpvar(1, T); P_gt sdpvar(1, T); u_gt binvar(1, T); % 燃机启停0-1变量 P_buy sdpvar(1, T); P_p2g sdpvar(1, T); P_ccs sdpvar(1, T); F_CH4_gt sdpvar(1, T); % 燃机消耗的总天然气燃料功率MW F_CH4_ext sdpvar(1, T); % 外购天然气燃料功率MW CH4_p2g sdpvar(1, T); % P2G甲烷化产甲烷功率MW F_H2_fuel sdpvar(1, T); % 掺烧氢气燃料功率MW H2_prod sdpvar(1, T); % 产氢能量MWh H2_meth sdpvar(1, T); % 甲烷化消耗氢能量MWh Soc_h2 sdpvar(1, T1); % 储氢罐状态MWh CO2_cap sdpvar(1, T); % CCS捕集CO2量t CO2_to_meth sdpvar(1, T); % 送往甲烷化的CO2量t E_extra sdpvar(1); E1_seg sdpvar(1); E2_seg sdpvar(1); E3_seg sdpvar(1); C []; % 约束集合 %% 电功率平衡 C [C, P_use_wind P_use_pv P_gt P_buy P_load P_p2g P_ccs]; C [C, 0 P_use_wind P_wind_avail]; C [C, 0 P_use_pv P_pv_avail]; C [C, 0 P_buy 200]; % 购电上限可根据外部电网容量调整 %% 燃气轮机约束 C [C, P_gt_min * u_gt P_gt P_gt_max * u_gt]; C [C, P_gt eta_gt * (F_CH4_gt F_H2_fuel)]; C [C, F_H2_fuel (alpha_max / (1 - alpha_max)) * F_CH4_gt]; for t 2:T C [C, -ramp P_gt(t) - P_gt(t-1) ramp]; end %% 天然气平衡 C [C, F_CH4_ext CH4_p2g F_CH4_gt]; C [C, F_CH4_ext 0]; %% P2G与储氢约束 C [C, 0 P_p2g P_p2g_max]; C [C, H2_prod eta_p2g * P_p2g * dt]; C [C, Soc_h2(1) 0.5 * H2_cap]; C [C, Soc_h2(T1) Soc_h2(1)]; for t 1:T C [C, Soc_h2(t1) Soc_h2(t) H2_prod(t) - F_H2_fuel(t) * dt - H2_meth(t)]; C [C, 0 Soc_h2(t1) H2_cap]; C [C, 0 H2_meth(t)]; end %% 甲烷化与CCS耦合 C [C, CH4_p2g eta_meth * H2_meth]; C [C, CO2_to_meth K_CO2_H2 * H2_meth]; C [C, CO2_to_meth CO2_cap]; C [C, CO2_cap capture_rate * F_CH4_gt .* EF_gas * dt]; C [C, CO2_cap CCS_capacity * dt]; C [C, P_ccs k_ccs_elec * CO2_cap]; %% 碳排放与阶梯碳交易 E_quota quota_intensity * sum(P_gt .* dt); E_real sum(F_CH4_gt .* EF_gas .* dt) - sum(CO2_to_meth); C [C, E_real 0]; C [C, E_extra E_real - E_quota]; C [C, E_extra 0]; C [C, E1_seg E2_seg E3_seg E_extra]; C [C, 0 E1_seg E1]; C [C, 0 E2_seg E2 - E1]; C [C, E3_seg 0]; C_carbon p1 * E1_seg p2 * E2_seg p3 * E3_seg; %% 目标函数 C_curtail c_cur * ( sum(P_wind_avail - P_use_wind) sum(P_pv_avail - P_use_pv) ) * dt; objective sum(price_buy .* P_buy) * dt ... price_gas * sum(F_CH4_ext) * dt ... 30 * sum(P_p2g) * dt ... % P2G运维成本 25 * sum(CO2_cap) ... % CCS运行成本 C_carbon ... C_curtail; %% 求解 ops sdpsettings(solver, gurobi, verbose, 2, mipgap, 0.001); diagnosis optimize(C, objective, ops); if diagnosis.problem 0 disp(求解成功); P_gt_opt value(P_gt); P_p2g_opt value(P_p2g); CO2_cap_opt value(CO2_cap); C_carbon_opt value(C_carbon); else disp(diagnosis.info); end这段代码虽然离“可直接投稿”还有距离但已经具备完整的MILP骨架。把预测数据换进去再补上结果画图就是一个能跑通的最小复现系统。4.3 结果输出与可视化建议求解完成后我一般会先输出四张图24小时电功率平衡堆叠图观察风电、光伏、燃机、购电如何共同满足负荷P2G耗电、制氢量、掺氢量和甲烷化量的时间曲线判断氢气的分配是否合理燃气轮机出力与CCS捕集量对照图看碳排放和碳捕集的时序关系阶梯碳交易成本构成饼图看哪一段碳价区间贡献了主要成本。可视化的最大价值不是“好看”而是能快速发现物理不合理的解。比如某小时P2G大量产氢但储氢罐SOC曲线没有同步变化那一定是平衡方程写错了再比如CCS捕集量在燃机停机时仍然大于0那一定是CO2_cap的上限约束没写全。5. 复现结果里几个值得细看的调度现象5.1 P2G-CCS耦合把弃风时段变成了“燃料生产时段”在基准场景下最容易观察到的现象是夜间负荷下降、风电出力较大且电价处于谷段时燃气轮机出力减少P2G开始增大耗电把多余风电转化成氢气。此时CCS如果还在运行捕集到的CO2不会白白封存而是送入甲烷化反应器与氢气合成甲烷。这些甲烷会存在气路里或直接补充到天然气管道等到早上负荷高峰时再被燃气轮机烧掉。这其实说明了一个很核心的调度逻辑P2G-CCS耦合不是简单“为了环保而环保”它是在电价低谷时用低成本电能合成高品位燃料把“弃风时段”变成了虚拟电厂自己的“燃料生产时段”。从经济账上看它同时节省了外购天然气费用也减少了弃风惩罚还降低了碳交易成本。5.2 阶梯碳价会同时改变燃机出力和CCS运行策略碳交易机制对调度的引导作用不是“一刀切”地减少燃机出力而是会重构整个系统的最优运行点。碳价提高后我观察到的变化包括燃气轮机在负荷峰段的出力适当下调部分电量由购电替代但燃机不会完全停掉因为P2G需要CO2来合成甲烷而CO2的主要来源就是燃机烟气CCS捕集率会向设备上限靠近因为捕集的CO2在系统中具有“原料价值”和“碳减排价值”双重收益风电商用空间增大弃风率明显降低。这一点对于理解“碳交易为什么能促进新能源消纳”很有帮助。碳价本质上给化石能源发电增加了一个额外成本而这个成本会沿着“燃机出力 → 碳排放 → 碳交易成本 → 总成本”的路径传导到每一度电上最终迫使系统寻找更便宜的替代方案。5.3 掺氢比例不是越高越好存在经济最优区间很多复现者做完基础算例就开始狂拉掺氢比例上限以为掺得越多越好。实际上在P2G产氢量受制于电解槽容量和弃风资源的条件下盲目调大掺氢上限通常会得到“掺氢比例约束根本没被触碰到”的结果。从成本结构看掺氢需要消耗P2G产出的氢气而氢气如果用于甲烷化可以替代外购天然气如果用于直接掺烧则替代的也是天然气但减排效果更直接。两者存在竞争关系。随着掺氢比例上限提高模型会在这两种氢气去向之间做权衡。当掺氢比例较低时提高上限能有效降低天然气成本和碳交易成本但当比例高到一定程度后继续提高只会把氢气从甲烷化那边抢过来反而可能减少甲烷产量导致外购天然气增加总成本不降反升。所以“最优掺氢区间”是系统运行的经济结果不是越大越好。5.4 对比场景设计一组能证明“耦合与碳交易都起了作用”的算例复现类文章最忌讳只有一个全局最优解。我建议至少跑下面四个场景场景是否启用P2G-CCS耦合是否启用阶梯碳交易表现结论S1否否基准系统弃风率最高S2是否消纳改善但低碳价值未被价格化S3否是碳成本高但无法回收CO2S4是是总成本与碳排放综合最优S2和S3的对比尤其重要。只看S2会以为“耦合本身已经很好了”只看S3又会以为“碳交易只是增加负担”。只有S4才能展现出“耦合为碳交易提供减碳手段碳交易为耦合提供经济回报”这个正反馈机制。这不只是论文里的“敏感性分析”更是理解标题里几个关键词之间关系的钥匙。6. 复现时踩过的五个坑与排查思路6.1 量纲不统一结果凭空多出两个数量级第一次跑出结果时我心里隐约觉得不对劲燃气轮机一天的发电量只有几百MWh但CO2捕集量居然算出了上万吨。排查到最后根源是我把“燃料功率MW”和“燃料能量MWh”混用了。F_CH4_gt这个变量我一边在燃气轮机出力方程里当成MW用一边又在碳排放公式里乘了EF_gas后当成MWh用量纲差了一倍结果自然乱套。排查量纲问题的方法很简单从头到尾检查每个约束的左右单位是否一致。功率项统一用MW能量项统一用MWh质量项统一用t。时间粒度为1小时时MW乘以1小时就是MWht/h乘以1小时就是t。把功夫花在单位表上永远不亏。6.2 掺氢比例约束写成了除法模型直接变成非线性最初我图省事把掺氢比例约束直接写成alpha(t) F_H2_fuel(t) / (F_CH4_gt(t) F_H2_fuel(t));Yalmip立刻提示Nonlinear constraints求解器也换成了非线性求解器不仅慢还不一定能找到全局最优。后来才意识到这个约束可以巧妙地线性化成F_H2_fuel(t) ≤ (alpha_max / (1 - alpha_max)) × F_CH4_gt(t)这其实是很多新手容易忽略的地方。看到除法和乘法出现在两个变量之间就应当马上思考能不能通过移项改写为线性形式。如果能尽量用线性约束替代否则求解效率和稳定性都会变差。6.3 碳交易分段变量上界漏写出现负碳价或“免费排放”还有一次碳交易成本算出来是负数导致系统拼命增加碳排放来“赚钱”。查完之后发现E1_seg和E2_seg的上界约束漏写了求解器把E_extra拆成了“第一段负数、第二段正数”的形式再用不同碳价一加就凭空造出了负成本。分段碳成本模型必须写全三个约束0 ≤ E1_seg ≤ E1 0 ≤ E2_seg ≤ E2 - E1 E3_seg ≥ 0一个都不能少。上界是物理分段边界下界保证碳价权重不会出现负值。这里没有捷径检查时务必逐行核对。6.4 储氢罐SOC没有闭环约束系统开始“凭空产氢”储氢罐的初始SOC和末端SOC必须要固定或者至少要让整个周期的充放电量满足闭环约束否则模型会把储氢罐当成一个无限能量的“魔法罐”。我最初只加了SOC上下限约束忘记加Soc_h2(T1) Soc_h2(1)结果最优解里储氢罐每天从0充到30 MWh又全部放空成本被严重低估看起来像是“每天都在白赚30 MWh氢气”。修复方式很简单在约束里加上Soc_h2(1) Soc_h2(T1)如果是调度周期重复执行的情况还可以把初始SOC作为优化变量但至少要在两个周期之间保持能量守恒不能出现能量凭空产生。6.5 调试策略先锁定耦合项再逐层放开最后分享一个非常实用的调试顺序。拿到这套模型后不要一上来就把P2G、CCS、掺氢、碳交易全开着跑。我的做法是分四步走先固定P2G和CCS出力为常数只优化燃气轮机、购电和新能源消纳验证电功率平衡与机组约束是否正确再放开P2G让制氢量成为变量检查氢气平衡和储氢罐SOC是否闭合再放开CCS和甲烷化让CO2捕集量与CO2_to_meth咬合观察碳平衡最后才把阶梯碳交易和掺氢比例上限全部加入。每一步都对比前后结果如果某一步出现明显不合理比如弃风率异常、某设备出力恒为0错误大概率就在刚放开的那一层里。这个顺序帮我省下了至少一个晚上的调试时间也建议你复现任何耦合系统时照做。在我个人看来这类EI论文复现真正有价值的收获不是最后那张漂亮的调度曲线图而是从头到尾把“物理过程 → 数学模型 → 代码约束 → 求解结果”这条链路走通的过程。只要单位统一、边界条件完整、耦合关系不丢哪怕不追求和原论文一模一样的结果这套模型本身也已经能作为后续研究的可靠基线了。