1. 为什么虚拟电厂开始盯上“氢”和“碳”:问题拆解与项目动机
1.1 单一储能救不了高比例新能源接入的场站
先聊个现实问题。过去几年我做虚拟电厂(VPP)相关项目时,最直观的感受是:如果 VPP 内部只有风电、光伏加电化学储能,优化调度做到最后往往会撞上一堵墙——弃电率压不下去,顶峰能力又不够。原因不复杂,电化学储能的容量和功率是绑定的,几小时尺度上的能量搬移还行,但跨天、跨周的调节需求它真的扛不住;而且一台储能电池的度电成本折算下来,在现货市场里未必跑得赢峰谷价差。
于是行业内开始把目光投向更长周期的能量载体。氢气就是其中一个方向。电解水制氢可以把富余电力转成氢,需要时再用氢发电回补电网,能量存储周期从小时级拉长到天级甚至周级。这就是所谓的 P2G(Power-to-Gas)环节。但问题是,纯 P2G 的往返效率并不理想,电转氢再氢转电,整体效率能到 40% 出头就算不错了,投入产出比看着就让人犹豫。
1.2 CCS 和掺氢为什么总被绑在一起讨论
把碳捕集(CCS)加进来之后,逻辑就变了。P2G 制氢需要大量电力,而 VPP 内部如果有燃气机组,燃气轮机排出的 CO₂ 浓度高、便于捕集,捕下来的 CO₂ 又能和氢气通过甲烷化反应合成天然气(也就是电转气-气转天然气的完整链条),再送回燃气机组烧掉。这一圈下来,碳排放被循环利用了,燃气机组的燃料成本也部分对冲了。
至于燃气掺氢,是另一个角度的改良。现有燃气机组直接烧纯天然气,碳排放强度是固定的;掺入体积比 10% 到 30% 的氢气后,同等发电量下的 CO₂ 排放量实打实下降。对于还不想大规模改造机组、又想快速降低碳强度的场景,掺氢是性价比很高的选项。但掺氢比例不是越高越好,它会影响燃气轮机的热值、火焰速度、燃烧稳定性,甚至机组寿命,所以调度模型里掺氢比必须作为变量、受约束地优化,而不是拍脑袋定一个数。
这个项目的标题之所以把 P2G-CCS 耦合和燃气掺氢放在一起,本质上是在做一个“多能互补、碳氢联动”的 VPP 内部资源整合问题:风电光伏多了,制氢;氢气多了,一部分掺进天然气烧掉发电,一部分和捕集的 CO₂ 合成天然气存起来;碳排放量又被阶梯碳交易成本约束着。整个系统里电、氢、气、碳四条线缠在一起,这就是优化调度模型复杂度和难度的来源。
1.3 阅读理解题:这个项目实际要解决哪几件事
我把标题拆成几个子问题,便于后面展开时对号入座:
- 虚拟电厂聚合了什么:风电、光伏、燃气轮机(含掺氢运行)、电转气装置、碳捕集装置、储氢罐、天然气储罐,以及与外网购售电的交互。
- 阶梯碳交易怎么建模:碳排放权的价格不是线性的,免费配额内低价、超过配额逐级升高,这种分段递增的成本函数会让调度策略显著偏向低碳运行。
- P2G-CCS 耦合意味着什么:P2G 的耗电、产氢,CCS 的捕碳、耗电,以及两者结合后氢和 CO₂ 合成天然气的产量,都存在物理上的量比关系。
- 掺氢燃气机组怎么进模型:燃气机组的出力-燃料消耗、碳排放特性会随掺氢比例变化,天然气管网/储气设施的混氢浓度也需要做物料平衡约束。
下面我把每个子问题的建模思路和 Matlab 实现路径掰开来讲,特别是阶梯碳交易那部分的分段线性化处理、P2G-CCS 单元的物料平衡约束怎么写,以及 Yalmip 环境下怎么把混合整数线性规划(MILP)问题稳定求解出来。
2. 阶梯碳交易机制建模:从阶梯电价逻辑到分段线性化的落地
2.1 阶梯碳价背后的政策含义与优化动机
阶梯碳交易的规则可以这样理解:主管部门给 VPP 分配一定的免费碳排放配额,实际排放量在配额以内不花钱;超出配额的部分,划分成若干个区间,越靠后的区间每一吨 CO₂ 的惩罚价格越高。这与居民阶梯电价是一个逻辑——用得越多,边际成本越高。
在优化调度模型里,这个机制带来的直接结果是:燃气机组每多发一度电,产生的碳排放成本不再是固定的,而是取决于当天总排放量落在哪个区间。这样调度模型在决策时,不仅会考虑“发一度电赚多少钱”,还会考虑“发这一度电是否会把排放量顶进更高价区间”。这就是阶梯碳交易对调度策略产生影响的核心路径——它让碳排放成本的边际值内生化了。
用具体的数来算一笔账。假设碳交易基价为 60 元/t,阶梯递增 20 元/t,排放区间长度 1000 t。免费配额 800 t,实际排放 2200 t。那么超出配额的部分为 1400 t。其中 0-1000 t 那一段,成本为 1000 × 60 = 60000 元;1000-1400 t 那一段,要看清这 400 t 落在第一级增档还是第二级。如果每个阶梯区间长度为 1000 t,则 1400 t 中落在第一级区间(0-1000)内的 1000 t 按 60 元算,剩余 400 t 按第二级 80 元算。但如果实际排放达到 2600 t,超过配额 1800 t,则 1000 t 按 60 元,800 t 按 80 元(这里以一级阶梯为例)。有的规则会设两三个档位,比如 60、80、100,超过越多越贵,这个逻辑是一致的。
这种分段函数的难点在于:它不是一个平滑的凸函数可以直接扔给求解器,需要做分段线性化,引入 0-1 整数变量来标识排放量落在哪个区间。后面第 4 节我会给具体的 Yalmip 代码框架。
2.2 免费配额分配方式与模型表达的取舍
模型里还需要处理免费配额的计算。常见做法有两种:一是按 VPP 内机组容量和基准年排放强度估算一个固定配额;二是按当年实际发电量乘以一个行业基准排放强度来算动态配额。前者在模型中是一个常数,处理起来最简单;后者虽然更贴近实际政策(比如按基准线法分配),但会让模型多一层非线性——配额本身依赖发电量,发电量是决策变量。
从复现论文的角度,我建议先用固定配额建立基线版本,跑通流程后再升级到动态配额。原因很简单:固定配额下,碳交易成本函数只是关于总排放量的分段函数,线性化手段成熟、稳定性好;动态配额会把这个分段函数的断点变成变量,问题规模扩大、求解难度明显上升。论文里画出来的调度结果图,绝大多数都是固定配额假设下得到的,这不影响结论的参考价值。
对于免费配额这个数值本身,建议项目里把它设为一个可调参数。比如说某 VPP 的配额量设置 2000 kg 和 800 kg,得到的调度策略会有明显差异,这正好可以作为灵敏度分析的素材,放到论文或项目报告里也是加分项。
2.3 总碳排放量的构成:别漏了外购电的间接排放
阶梯碳交易建模还有一个容易踩坑的点:排放量统计范围。VPP 内部的碳排放来源不只是燃气机组燃烧,还包括从电网购电对应的间接排放。虽然 VPP 是用户侧主体,但如果项目背景是园区型 VPP,外购电的间接排放常常被纳入核算边界。处理方式很简单:购电量乘以电网平均排放因子,折算成等效 CO₂ 排放量,加进总排放量里一起参与阶梯计价。
需要注意的是,如果模型中允许 VPP 向电网售电,那么售电部分对应的排放不应计入 VPP 的排放账户,否则会出现逻辑矛盾:同一度电,发电侧计了排放、用户侧又计了排放。在代码实现时,只需用“净购入电量 = 购电 - 售电”乘以排放因子,或者干脆只对购电量计排放,售电部分不做重复计算。我在复现一些开源代码时就发现过这个细节没处理好导致结果异常的情况。
3. P2G-CCS 耦合单元建模:物料平衡、能量流与运行约束
3.1 从电解水到甲烷化的全链条物质流动
P2G-CCS 耦合单元是这个 VPP 里的“化学反应核心”。完整链条分三步:
- 电解水制氢:P2G 装置消耗电力,将水电解为氢气和氧气。这里有一个约束:P2G 的耗电功率和产氢速率近似成正比,模型里用一个转换效率 η_el 连接两者即可。
- CCS 捕碳:CCS 装置从燃气轮机排烟中捕集 CO₂,需要消耗一定电能(再生能耗、压缩能耗等)。捕集到的 CO₂ 一部分可以封存或外售,另一部分送入甲烷化反应器。
- 甲烷化合成天然气:氢气和 CO₂ 在反应器内发生 Sabatier 反应:CO₂ + 4H₂ → CH₄ + 2H₂O。按化学反应计量关系,1 mol 甲烷需要 4 mol 氢气加 1 mol 二氧化碳。换算成质量流量,每合成 1 kg 甲烷大约需要 0.5 kg 氢气和 2.75 kg CO₂。
这个计量关系非常关键。很多人在建模时直接把 P2G 产氢量当作“氢负荷”供给燃气轮机掺氢使用,忽略了甲烷化还需要消耗一部分氢;或者捕集了 CO₂ 却不用,白白增加了 CCS 的运行能耗。在耦合系统里,氢气要优先分配给甲烷化反应,剩余氢再进入储氢罐或掺氢燃烧,否则物料不平衡,整个模型跑出来的“优化结果”物理上是不可行的。
3.2 运行域与耦合约束的数学表达
假设一个调度周期为 24 h,时间间隔 1 h,决策变量包括:
- P2G 耗电功率 P_P2G(t),上限为 P_P2G_max,下限为 P_P2G_min(电解槽有最小稳定运行功率,不能低于某个值,否则启动/停机的损耗过大)。
- 产氢流量 F_H2_el(t) = η_el × P_P2G(t) / HHV_H2,其中 HHV_H2 可取 39.4 kWh/kg,用来把电功率换算成氢气质量流量。
- CCS 捕碳量 F_CO2_capture(t) = η_ccs × F_flue(t),其中 F_flue(t) 是燃气轮机排烟 CO₂ 流量,η_ccs 是捕集率,一般取 0.85 到 0.95。CCS 的耗电功率可以近似为 P_CCS(t) = λ_CCS × F_CO2_capture(t),λ_CCS 是单位捕碳能耗,常见范围在 0.2-0.5 kWh/kg CO₂。
- 甲烷化消耗氢气 F_H2_meth(t) = 4/1 × (F_CH4_meth(t)/M_CH4) × M_H2,用摩尔质量换算后,约等于每生产 1 kg CH₄ 消耗 0.5 kg H₂。
- 甲烷化消耗 CO₂ F_CO2_meth(t) = 1/1 × (F_CH4_meth(t)/M_CH4) × M_CO₂,约等于每生产 1 kg CH₄ 消耗 2.75 kg CO₂。
- 甲烷化产品 F_CH4_meth(t) 送入天然气储罐或直接供燃气机组使用。
氢平衡约束可以写为:
F_H2_el(t) + F_H2_storage_out(t) = F_H2_meth(t) + F_H2_blend(t) + F_H2_storage_in(t)
其中 F_H2_blend(t) 是掺入燃气轮机燃料的氢气流量。储氢罐模型用简单的能量/质量状态方程:
S_H2(t+1) = S_H2(t) + F_H2_storage_in(t) × Δt - F_H2_storage_out(t) × Δt
同时限制 S_H2(t) 在上下限之间,以及储氢罐的进出速率不超过最大充放速率。
CCS 捕集量与甲烷化用碳量的平衡则写成:
F_CO2_capture(t) ≥ F_CO2_meth(t)
剩余的 CO₂(捕了但没用掉的)可以封存,有封存成本;也可以外售给其他工业用户,有售碳收益。模型里把这个差值乘以一个单位封存成本或者外售价格即可。
3.3 为什么 CCS 不能随便“满功率”运行
CCS 是耗电大户。捕集率越高,对应的电耗和再生能耗越高。在 VPP 里,这就变成了一个经济权衡:燃气机组发一度电赚的钱,减去这一度电对应的额外碳排放成本,再减去捕碳耗电造成的少发电或者购电成本,剩下的才是净收益。
实际调度结果往往会出现这样的模式:在碳价较低、电价较高的时段,CCS 捕集率会调低甚至停机,因为直接缴碳税更划算;在碳价高、电价低(尤其是光伏大发的中午)时段,CCS 反而会开足马力,因为捕下来的碳不只是减排,还能和光伏制氢合成天然气储存起来。这种“时变运行策略”是耦合系统相比独立运行的最大优势,也是调度模型价值所在。
我第一次跑这种模型时,给 CCS 设置了一个固定捕集率,结果优化器完全不用它,因为捕碳电耗折算成成本后,高于直接买碳配额的成本。后来把捕集率改成 0.5-0.95 的连续变量,结果才合理起来。所以代码里最好把 η_ccs 作为决策变量,而不是常量。
4. 掺氢燃气机组建模与网络约束:热值修正、双燃料平衡和 Weymouth 线性化
4.1 掺氢比是决策变量时,机组模型怎么写
燃气轮机掺氢后,燃料由天然气和氢气组成。掺氢的体积比定义为:
θ_blend(t) = V_H2(t) / (V_CH4(t) + V_H2(t))
由于氢气和天然气的热值差异很大,机组发电量对应的燃料总热量需求应写为:
P_GT(t) / η_GT = F_CH4(t) × LHV_CH4 + F_H2_blend(t) × LHV_H2
其中 η_GT 是机组发电效率,LHV_CH4 取 55.5 MJ/kg,LHV_H2 取 120 MJ/kg(氢的 LHV 数值上是天然气的两倍多,质量基)。这个公式是双燃料热平衡的核心,不能写成“天然气消耗量 + 氢气消耗量”简单相加,热值不同,必须各自乘各自的热值。
掺氢比对机组运行的影响体现在三方面是:氮氧化物排放特征变化、燃烧稳定边界(贫燃熄火极限)、以及机组出力爬坡速率限制。在简化调度模型中,最常见的方式是对掺氢比设置上下限,比如 θ_blend(t) ∈ [0, 0.3],表示体积掺氢比不超过 30%。这个上限的取值依据来自燃气轮机制造商的保证值和实际改造案例的保守估计。
如果要更精细一些,可以定义一个二次效率修正项:η_GT_eff(t) = η_GT_base × (1 - α × θ_blend(t)),其中 α 是掺氢对效率的惩罚系数。但这个式子引入非线性后,求解难度会上升。如果论文目标是 MILP 求解,建议把修正项线性化,或者干脆忽略这个效率变化,只在约束里限制掺氢比上限。我在做工程化项目时更倾向于后者,因为掺氢比不超过 20% 时,效率变化通常小于 1 个百分点,对调度结果的影响远小于碳价波动的影响。
4.2 燃气机组排放核算:掺氢的减排乘数
燃气机组碳排放的函数是:
E_GT(t) = F_CH4(t) × EF_CH4 + F_H2_blend(t) × EF_H2
其中 EF_CH4 是天然气燃烧的 CO₂ 排放因子(kg CO₂/kg CH₄,大约 2.75),EF_H2 是氢气燃烧的 CO₂ 排放因子(0)。这样掺氢后单位发电量的碳排放强度自然下降,减排效果直接体现在排放量中。
注意一个细节:燃气掺氢燃烧虽然降低了机组端的碳排放,但如果氢气来自电解水制氢,而电解水用电可以来自光伏风电,也可以来自电网购电。如果是电网购电制氢,那么制氢过程中的间接排放应该追溯计入。也就是说,氢气并不是“零碳”的,它的碳足迹取决于制氢电力的来源。在 VPP 模型中处理方式是:制氢耗电从 VPP 内部新能源出力中优先供给,不足部分从电网购电补足,这样既有物理意义,又不会把模型复杂化。
4.3 储气系统与管网掺氢约束的处理
VPP 内部的天然气储罐和储氢罐需要分别建模,但如果存在管网级掺氢,还需要对管网中的混氢浓度做约束。在一些精细模型中,管网的动态需要 PDE(偏微分方程组)描述,这会极大增加求解难度;但在日前调度这类离线优化场景,通常的做法是使用稳态 Weymouth 方程加节点混氢浓度约束。
稳态 Weymouth 方程描述了天然气管道流量与两端压力差之间的关系,是非线性的。在 MILP 框架下,可以把它近似为一组分段线性约束,或者干脆在 VPP 模型中把管网简化为“无压降的虚拟集中节点”,只保留节点物料平衡。对绝大多数 VPP 调度研究而言,集中式储气模型就已经足够支撑结论,不需要引入管网动态。如果你的场景特别关注管网约束,我建议单独拆出来做一个小案例验证,不要一开始就放进大模型。
在代码中,储气罐的状态约束建议写成:
S_gas(t+1) = S_gas(t) + F_meth(t) - F_gas_GT(t) - F_gas_sell(t)
其中 F_meth(t) 是甲烷化产气,F_gas_GT(t) 是供燃气机组的天然气流量,F_gas_sell(t) 是外售天然气。天然气外售在 P2G 经济性评估中很重要,因为甲烷化产气成本不低,如果 VPP 有富裕气量,外售可以获得额外收益。这个通道不要漏掉。
5. 虚拟电厂优化调度全模型:目标函数、约束集与决策变量清单
5.1 目标函数拆解:运行收益减综合成本
虚拟电厂调度问题的目标函数一般写成最大化净收益,包括售电收入、售氢收入、售气收入、碳交易成本、运维成本、购电成本等。展开写:
max Σ_t [ Price_elec(t) × P_sell(t) - Price_buy(t) × P_buy(t) + Price_H2 × F_H2_sell(t) + Price_gas × F_gas_sell(t) - C_OM_all(t) - C_carbon(t) ]
其中 C_carbon(t) 是阶梯碳交易成本,按第 2 节的分段函数计算;C_OM_all(t) 是各设备运维成本,可以简化为与出力/产气量成正比的比例系数;Price_elec(t) 是分时电价,需要从项目数据或假设中给定。
在源代码实现时,我建议把收益项和成本项分开累积,方便事后做敏感性分析。比如售电收入、售氢收入、碳成本各自存成一个数组,最后画堆叠图非常直观。这也是审稿人喜欢看的形式,“收益结构分析”差不多是这类论文的标配图表。
5.2 约束集分类:功率平衡、爬坡约束、状态变量边界
约束集按角色可以分为以下五类:
- 电力功率平衡:P_wind(t) + P_pv(t) + P_GT(t) + P_buy(t) + P_storage_discharge(t) = P_load(t) + P_P2G(t) + P_CCS(t) + P_sell(t) + P_storage_charge(t)。这里电储能可选,如果模型里没有电储能,就把相关项删掉。
- 燃气机组运行约束:出力上下限、爬坡约束、最小启停时间(如果引入机组组合变量)。
- P2G-CCS 耦合约束:产氢-耗电关系、捕碳-耗电关系、甲烷化物料平衡(见第 3 节)。
- 储氢/储气动态约束:状态转移方程和容量、充放速率限制。
- 碳排放与阶梯碳交易约束:总排放量计算、分段区间标识变量、碳成本线性化(见第 4 节)。
如果模型进一步考虑备用约束(旋转备用),需要加一条:可调节出力(燃气机组 + 电储能 + P2G 调节能力)在上调方向上大于等于系统备用需求。P2G 和 CCS 在备用方面其实有一定调节作用,但绝大多数论文把它当纯刚性负荷处理,我建议先不加备用约束,等主模型跑通后再扩展。
5.3 变量规模与求解耗时预估
以 24 h、每小时一个点计算,决策变量数量大约为:
- 连续变量:约 14 类 × 24 + 储氢、储气状态 2 × 24,总共约 400 个。
- 整数变量:阶梯碳交易分区指示变量 3 × 24 = 72 个,如果涉及机组启停再加 24 个。
这个规模对 CPLEX 或 Gurobi 来说很小,秒级甚至亚秒级就能求解。真正影响求解速度的不是变量数量,而是分段线性化时引入的大 M 法和整数变量个数。如果阶梯碳交易的档位设了 5 档,每个时段 5 个 0-1 变量,24 h 就是 120 个,也还能接受。但如果你打算做全年 8760 h 的长周期仿真,那就得考虑简化:要么把时段聚合成典型日,要么用启发式算法替代 MILP。我的建议是,先把 24 h 的 MILP 模型跑得稳定透彻,再谈规模化。
6. Matlab 实现核心代码框架:Yalmip + 求解器配置与关键语句解读
6.1 环境准备和求解器选择
Matlab 环境下做数学优化建模,目前最顺手的方式还是Yalmip。它本身不是求解器,只是一个建模层,底层调用 CPLEX、Gurobi 或开源求解器。我个人习惯用 Gurobi,因为它在 MILP 上的表现稳,MIP Gap 收敛快。如果考虑版权,SCIP 也能在 Yalmip 里无缝调用。
安装方面,常规流程是:下载 Yalmip 源码包放入 Matlab 路径 → 下载求解器安装包并配置 license → 在 Matlab 中运行yalmiptest验证。Tutorial 里有大量现成例子,建议先跑一遍optimize基础例程确认建模和求解链路通畅。
6.2 阶梯碳交易分段线性化的 Yalmip 实现
阶梯碳交易成本函数的分段线性化是这块代码里最需要小心的部分。基本思路是引入 0-1 变量标识总排放量所在的阶梯区间,区间内排放量乘以对应价格累加。核心代码框架如下:
% 阶梯碳交易参数 e0 = 800; % 免费配额,kg step_len = 1000; % 阶梯区间长度,kg price0 = 0.06; % 配额内价格(通常为0或很小),元/kg price_step = [0.08, 0.10, 0.12]; % 超出配额后各区间价格,元/kg n_step = length(price_step); E_total = ...; % 总排放量,由其他变量计算得到,连续变量 % 引入分区变量 delta = binvar(n_step, 1); % 每个区间是否启用 E_seg = sdpvar(n_step, 1); % 落在每个区间的排放量 % 约束:总排放 = 配额内部分 + 各区间排放量 E_total = e0 + sum(E_seg); % 区间限制:区间 k 被启用时,E_seg(k) <= step_len;未启用则 E_seg(k) = 0 M = 5000; % 大M值,应大于所有可能的区间排放量 for k = 1:n_step constraints = [constraints, E_seg(k) >= 0]; constraints = [constraints, E_seg(k) <= step_len * delta(k)]; constraints = [constraints, E_seg(k) <= (E_total - e0 - sum(E_seg(1:k-1)))]; constraints = [constraints, E_seg(k) >= (E_total - e0 - sum(E_seg(1:k-1))) - M * (1 - delta(k))]; end这一段是“顺序启用”的分区表达——第 k 个区间只有在第 k-1 个区间被填满后才可能被启用,通过累计排放量约束来保证区间之间的先后顺序。实际实现时,一种更稳妥的写法是用 SOS2 约束,Yalmip 有现成支持:
% 用SOS2实现分段线性函数 E_axis = [0, 1000, 2000, 3000] + e0; price_axis = [0.06, 0.08, 0.10, 0.12]; lambda_ = sdpvar(length(E_axis), 1); constraints = [constraints, E_total == lambda_' * E_axis']; constraints = [constraints, sum(lambda_) == 1, lambda_ >= 0]; constraints = [constraints, carbon_cost == lambda_' * (E_axis .* price_axis)'];SOS2 的好处是不用手动写大 M 约束,求解器原生支持特殊有序集,数值稳定性更好。上面大 M 的写法是为了帮助你理解机制,建议在最终代码中优先用 SOS2 或 Yalmip 自带的piecewise函数。
6.3 燃气掺氢约束与储氢动态的代码示例
燃气机组燃料平衡用sdpvar定义变量后直接写约束即可:
% 燃气机组双燃料热平衡 % F_CH4_GT: 天然气流量, F_H2_GT: 掺氢流量, P_GT: 机组出力 constraints = [constraints, P_GT / eta_GT == F_CH4_GT * LHV_CH4 + F_H2_GT * LHV_H2]; % 掺氢体积比约束(换算成质量比后再约束) % 体积比 = (F_H2_GT / rho_H2) / (F_CH4_GT / rho_CH4 + F_H2_GT / rho_H2) % rho_CH4=0.717 kg/m3, rho_H2=0.09 kg/m3 vol_H2 = F_H2_GT / 0.09; vol_CH4 = F_CH4_GT / 0.717; constraints = [constraints, vol_H2 <= 0.30 * (vol_H2 + vol_CH4)];掺氢约束写的注意点是:题目要求通常是体积比,而模型里计算用的是质量流量,所以必须通过密度的换算或者直接用量热值表达。有的文献直接用“掺氢能量比”,那就变成热值比例,公式会稍微不同,不要混用。
储氢罐的动态用等式约束表达:
% 储氢罐:S_H2下一时刻 = 当前储量 + 充入 - 放出 for t = 1:T-1 constraints = [constraints, S_H2(t+1) == S_H2(t) + F_H2_in(t) - F_H2_out(t)]; end constraints = [constraints, S_H2 >= S_H2_min, S_H2 <= S_H2_max]; constraints = [constraints, F_H2_in >= 0, F_H2_in <= F_H2_in_max]; constraints = [constraints, F_H2_out >= 0, F_H2_out <= F_H2_out_max];注意sdpvar定义的向量索引从 1 开始,而时段编号也从 1 开始,所以 S_H2(1) 对应初始时段开始时的储量,S_H2(T+1) 对应调度周期结束时的储量,需要额外定义一下。这类“时间索引错位”问题在调试时最容易让人困惑,建议在代码开头用注释明确每个向量的时间语义。
6.4 求解器调用及后处理绘图
求解主流程代码如下:
ops = sdpsettings('solver', 'gurobi', 'verbose', 2, 'debug', 1); result = optimize(constraints, -objective, ops); if result.problem == 0 % 求解成功,提取变量 P_GT_value = value(P_GT); F_H2_GT_value = value(F_H2_GT); E_total_value = value(E_total); else disp(['求解失败: ', result.info]); end后处理建议画四类图,基本可以覆盖论文/项目图表的全部需求:
- 电功率平衡堆叠图:风电、光伏、燃气、购电 vs 负荷、P2G、CCS、售电。
- 掺氢比和燃气机组出力时序曲线:看掺氢策略随电价波动的变化。
- 储氢罐/储气罐储量曲线:看能量存储的充放节奏。
- 碳排放量与阶梯碳成本分段展示:柱状图叠加阶梯价格区间,直观展示碳成本构成。
这些图用 Matlab 的plot、area、bar函数就能完成,关键是数据组织得好——建议在求解后立刻把value(...)的结果存成 struct,比如results.P_GT、results.F_H2_GT,然后单独写一个plot_results.m脚本,不要在求解脚本里画图,否则调试过程会非常痛苦。
7. 典型算例设计与结果分析:从数据构造到调度策略解读
7.1 算例参数怎么构造才合理
一个能支撑论文结论的算例,参数不能随意拍脑袋。我的建议是按以下思路构造:
- 新能源出力曲线:用某典型日的风电、光伏归一化出力数据,乘上装机容量。数据可以用 Matlab 自带
load某个示例,或者自己生成一个带时序特征的曲线(光伏中午高、晚上零;风电夜间高、白天波动)。 - 负荷曲线:用园区典型工作日的负荷模式,峰在白天、晚峰,谷在凌晨。
- 分时电价:参考国内某省现货市场出清价或峰谷电价结构,设置 3-4 个时段的价格档。
- 碳排放参数:免费配额、阶梯区间长度、各区间价格,参考试点碳市场的实际参数数量级。
具体数值我列一个建议表:
| 参数 | 数值 | 备注 |
|---|---|---|
| 风电装机 | 100 MW | 典型日出力曲线 |
| 光伏装机 | 80 MW | 典型日出力曲线 |
| 燃气轮机装机 | 60 MW | 掺氢上限 30% |
| P2G 额定功率 | 30 MW | 效率 η_el = 0.75 |
| CCS 捕集能力 | 15 t/h | 捕集率 0.85-0.95 |
| 储氢罐容量 | 5 t | 初始储量 2.5 t |
| 天然气储罐 | 30 t | 初始储量 15 t |
| 免费碳配额 | 800 t/d | 基准线法 |
| 碳交易基价 | 60 元/t | 每档递增 20 元/t |
7.2 对照实验设计:四种场景才好讲故事
只跑一个场景的数据在论文或项目结论里没有说服力,建议至少做四个对照场景:
- 场景 A:无 P2G-CCS、无掺氢、无阶梯碳交易(传统 VPP,固定碳价)。
- 场景 B:有 P2G-CCS,无掺氢,阶梯碳交易。
- 场景 C:有 P2G-CCS、有掺氢,固定碳价(无阶梯)。
- 场景 D:全部技术都加,阶梯碳交易。
这四组对比跑完之后,可以清晰拆解出各项技术的贡献:P2G-CCS 对弃风弃光的抑制效果、掺氢对碳排放强度的降低、阶梯碳价对燃气机组出力的约束。我在项目里跑过类似的对照,结果基本都指向同一个结论:单项技术都有收益,但叠加后存在耦合增益,也就是“1+1>2”的效果。这种结论很有项目报告价值,也容易在答辩时展开。
7.3 调度结果怎么解读:三组典型对比
以我实际跑的算例经验,几个值得关注的规律:
- 午间光伏大发时段:电价低、光伏出力高,此时 P2G 大概率满负荷运行,把剩余电力转化为氢气;CCS 捕集的 CO₂ 与氢气共同合成天然气,为晚高峰储备“燃料”。表现在功率曲线上,午间 P2G 耗电功率和甲烷化产气率都处于高位。
- 晚高峰时段:电价高、负荷高,燃气机组满发,掺氢比可能达到上限 30%(因为高电价下多发电的收益高于掺氢带来的效率损失)。此时储氢罐快速放空,天然气储罐也可能快速出清。
- 碳价敏感时段:当总排放量逼近阶梯分界点时,优化器会主动调整燃气机组出力——可能在某个时段让机组少发甚至停机,转为从电网购电,避免排放量“跨档”导致整体碳成本跳升。这个行为在结果里表现为:机组出力曲线在某个时段出现异常下降,对应碳成本曲线出现拐点。
这种“为控制碳排放总量而牺牲局部发电收益”的策略,在固定碳价模型里是看不到的,也恰恰是阶梯碳交易模型的核心价值所在。如果在结果图里发现不了这类现象,建议回头检查碳交易分段约束是否真的生效、大 M 值是否设置过大导致分区变量失效。
8. 复现这个项目时最常踩的坑:排查链路与代码调试心得
8.1 坑一:阶梯碳交易的“区间未启用但仍有排放量”
这是我在复现类似代码时遇到的第一个大坑。现象是:总排放量明明没超过配额,碳交易成本却出现了负值或非零值;或者是第三个区间从未启用,但 E_seg(3) 依然有数值。
原因是分段约束写得不严谨。经典的错误写法是只约束E_seg(k) <= step_len * delta(k),没有约束未启用区间的变量必须为 0 的下界关系。加上E_seg(k) >= -M * (1 - delta(k))也没用,因为排放量是非负的。正确的做法需要保证:当 delta(k)=0 时,E_seg(k) 必须恰好为 0。这时需要引入激活标识,或者检查大 M 约束里是否存在目标函数对 E_seg 有正向激励导致它“偷跑”的可能。
调试方法:把优化结果打印出来,逐时段检查 delta 和 E_seg 的对应关系。如果发现矛盾,在约束集里临时添加E_seg(k) <= M * delta(k)且E_seg(k) >= 0,并确认目标函数不会“奖励”多余的 E_seg。若目标函数中碳成本是 E_seg 的增函数,理论上优化器不会主动增加 E_seg,但如果约束没有完全锁死自由度,数值求解器可能给出不可行但“看起来正常”的解。所以建议碳成本项写成sum(lambda_ .* price_axis) * E_total而不是sum(E_seg .* price_step),这两种表达式在分段线性化中是等价的,但后者的数值敏感度更高。
8.2 坑二:掺氢约束引入非线性导致求解器报错
用 Yalmip 写掺氢体积比约束时,如果写成:
F_H2_GT / (F_CH4_GT + F_H2_GT) <= 0.3;这就是一个非线性约束,Yalmip 会尝试调用非线性求解器(如 fmincon、ipopt),大概率在小规模问题上能跑,但速度慢、稳定性差,而且很难断言全局最优。
解决思路有两条。一是把约束改写成线性形式:F_H2_GT <= 0.3 * F_CH4_GT + 0.3 * F_H2_GT,移项可得0.7 * F_H2_GT <= 0.3 * F_CH4_GT。如果掺氢比按体积算,还需要通过密度换算,但数学形式同样是线性的。二是在建模前把掺氢比固定为几个离散档位,用整数变量选择,这样模型是 MILP,求解更干净。我的经验是:只要论文结论不依赖掺氢比连续变化,优先用离散档位方案,省心且稳定。
调试时如果 Yalmip 报“Nonlinear constraints detected”,可以在sdpsettings里加'solver', 'gurobi'强制指定,如果不能求解,Yalmip 会直接报错告诉你问题类型不匹配,这正是排查信号。
8.3 坑三:P2G 和 CCS 同时高运行时电力平衡失稳
在午间光伏大发的情况下,优化器倾向于让 P2G 和 CCS 同时高负荷运行。如果模型里没有约束“P2G + CCS + 负荷 <= 新能源 + 购电上限”,求解器会给出一个物理上不可行的方案——比如购电功率超过联络线容量。这个坑在能量平衡约束里通常不会暴露,因为等式约束本身包含了变量耦合,但如果你把 P2G 最大耗电功率设得过大,且购电上限没有正确设置,就会出现“虚拟电厂从电网买几百 MW 电去制氢”的荒唐结果——从经济角度看可能不合理,但从约束角度看它合法。
改进办法是设置联络线功率上下限:
constraints = [constraints, P_buy <= P_line_max, P_sell <= P_line_max]; constraints = [constraints, P_buy + P_sell_source <= P_load + P_P2G + P_CCS + P_sell];为了避免这种结果,还需要在目标函数里为购电设置足够大的成本系数、为售电设置合理的价格。如果模型里多了“多余的电量可以白白送给 P2G 制氢而没有任何电费压力”,那问题一定出在目标函数缺少购电成本惩罚项。检查目标函数中是否遗漏了Price_buy(t) * P_buy(t)这项。
8.4 调试工作流:从不可行解到收敛最优解的排查顺序
碰到求解器报不可行(infeasible),不要急着去翻约束,我按以下顺序排查:
- 用
yalmiptest确认 Yalmip 与求解器连接正常。 - 注释掉目标函数,只求可行性(
optimize(constraints, [], ops)),如果可行,再逐步加目标函数。 - 如果不可行,用
constraints的diagnostic属性定位不可行约束,但更高效的办法是二分排除法:删除一半约束,看是否可行,逐渐缩小范围。 - 重点检查时间索引:状态变量 S(t+1) 与 S(t) 的错位、初值定义、末值约束是否相互矛盾。
- 检查大 M 值:M 太大(超过 1e6)会导致数值病态,M 太小(小于最大可能的区间排放量)会导致有效解被错误排除。
我有个习惯:在首次跑通前,把所有不等式约束的右端项都放宽 20%,先把模型跑通,再逐步收紧到目标值。这样能最快定位是哪类约束导致不可行。
9. 扩展方向:这个模型还能往哪些方向走
如果基本模型已经跑通,想继续深入,可以考虑以下几个方向,每个方向都有明确的论文或工程价值:
- 多时间尺度协同:日前调度 + 日内滚动修正 + 实时反馈。P2G-CCS 和储氢罐的时间常数比电储能大得多,多时间尺度模型能更好处理这种“慢系统”和“快系统”的耦合。
- 源-荷不确定性:风电光伏出力和负荷都有预测误差,引入鲁棒优化或随机规划(场景法)可以让调度策略更实用。阶梯碳交易和鲁棒优化结合时,需要注意分段函数的鲁棒形式化处理。
- 电-碳-绿证市场协同:绿证交易和碳交易存在联动,VPP 的综合收益还应考虑绿证收益。把绿证纳入目标函数后,光伏风电的“环境价值”会进一步凸显,调度结果会更强调新能源消纳。
- 需求响应机制的引入:VPP 内部如果有可调负荷(如电解铝、充电桩),需求响应可以和 P2G 制氢形成竞争关系——可调负荷减少用电时,P2G 获得更多低价绿电;反之亦然。这需要在模型中增加负荷弹性约束。
扩展方向的选择取决于你的项目周期和最终目标。如果只是课程项目或结题报告,建议先把基本模型做扎实,把四组对照实验跑完,加上灵敏度分析(碳价、配额、掺氢上限三个维度),已经足够撑起一篇结构完整的报告。如果是做期刊论文,鲁棒优化或市场协同方向更值得投入。
10. 关于设备建模参数取值的一些经验参考
在复现过程中你会发现,参数取值的合理性往往比求解算法的选择更影响结果的可信度。根据我在实际项目和文献阅读中积累的数据,给出几个常用参考范围:
- 电解槽效率:碱性电解槽 60%-75%,PEM 电解槽 65%-85%。调度模型里取 70%-75% 比较稳妥。
- CCS 捕集能耗:胺法燃烧后捕集,再生能耗约 3.5-4.0 GJ/t CO₂,折合电耗约 0.3-0.5 kWh/kg CO₂。先进溶剂或膜法可以更低,但工业级保守取 0.4。
- 甲烷化效率:热量损失考虑在内,化学计量转化率接近 100%,但整个 P2G 链条效率在 40%-60% 之间。模型里甲烷化直接用量比关系即可,不需要单独设效率系数。
- 燃气轮机掺氢极限:目前主流重型燃机(如 GE、西门子的 H 级机型)掺氢比例可以到 30%-50%,但多数在役机组或改造项目的保守值在 15%-20%。模型里设 20% 作为默认,30% 作为上限情景。
这些参数在论文的“算例设置”部分只需要一句话描述来源,在代码里则是常量定义。建议把它们集中放在一个parameters.m脚本里,方便批量修改和灵敏度分析。
11. 一些个人体会和实操建议
最后聊点代码之外的感受。这类“多技术耦合 + 优化调度”的题目,难点从来不在某个具体技术的建模,而在于多个时间常数、多种能量载体、多个市场机制交织后的系统性协调。P2G-CCS 是一个“吃电、产氢、耗碳、出气”的四端元件,燃气掺氢让本来只有“电-气”的耦合又多了一个维度,阶梯碳价格又把碳排放从一个环境约束变成经济约束——这种多重耦合关系如果只靠脑子想,很容易在某个细节上漏掉关键的物理约束或经济变量。
我的实操建议是:拿到题目后先不急着写代码,先在草稿纸上画出能量流和物料流图,标明每个设备的输入输出变量,逐个写出量纲自洽的方程。这个工作比调代码花费的时间更长,但价值也更高。等所有方程都在纸面上对齐之后,Matlab 建模不过是机械的翻译工作。
另外,养成“每加一组约束就跑一次”的习惯。不要等把所有约束写完了再调试,那会让错误定位变得极其困难。我通常先跑“纯电模型”(新能源 + 燃气 + 电储能 + 碳成本),验证基本结果合理后,再分阶段加入 P2G、CCS、掺氢、储氢罐。每加一个模块就跑通一次,同时检查新模块变量在结果中的取值是否符合物理直觉。
这个项目如果完整复现下来,无论是用于课程设计、毕业论文还是一次内部技术汇报,收获都不会小。毕竟它把电、氢、气、碳四条线非常紧凑地串在了一个虚拟电厂的框架里,你调通它的那天,基本也就理解了综合能源系统优化最有代表性的那一类问题。