考虑碳捕集与电转气的虚拟电厂优化调度Matlab实现

这个题目在“双碳”背景下算是很典型的组合调度问题：虚拟电厂（VPP）把垃圾焚烧、碳捕集、电转气（P2G）以及常规的风光储聚合成一个整体，最后用Matlab去求解最优调度策略。项目名字看着长，拆开其实就三条主线——VPP负责聚合，碳捕集负责把垃圾焚烧烟气的CO2“扣”下来，P2G负责把电能转成天然气再进网。我做这套模型前后花了两个多月，中间踩了不少坑，这篇就把它背后怎么建模、怎么用Matlab实现、怎么调试的问题一次性说清楚。

这个项目适合两类人看：一类是做虚拟电厂优化调度、综合能源系统方向的研究生或工程师，尤其是想用Matlab把碳捕集、P2G这类模型快速跑起来做算例分析的人；另一类是电力行业里做碳交易、负碳技术评估的从业者，想搞清楚垃圾焚烧电厂配碳捕集之后，运行成本到底怎么算、P2G加进来到底能带来多少协同收益。对有一定Matlab基础但不熟悉优化建模的读者，我也尽量把约束的来龙去脉写明白。

先说结论：这套调度模型最终的落点是混合整数线性规划（MILP），目标是最小化系统总运行成本，包括燃料成本、启停成本、购售电成本、碳交易成本，还要满足功率平衡、机组出力爬坡、碳捕集运行范围、P2G转换效率等一系列约束。实测下来，用YALMIP+Cplex求解24小时调度问题，规模大概几十个连续变量加几十个0-1变量，求解时间通常在几秒到十几秒，完全在可接受范围内。

下面逐个拆开讲。

1. 这种配置的价值：为什么要把碳捕集、电转气、垃圾焚烧绑在一起

1.1 VPP设计视角：聚合产生的协同收益

虚拟电厂本身不改变物理设备的接线方式，它靠的是信息通信和优化调度把分散的分布式能源聚合起来，对外作为一个整体参与电力市场或接受电网调度。传统的VPP模型通常只聚合风电、光伏、储能、燃气轮机，这种组合的问题在于：调节能力有限、碳排放没有被系统性地管起来。

加入垃圾焚烧机组之后，系统的基荷稳定性明显提高，因为垃圾焚烧不像风电光伏那样完全看天吃饭；但它有一个天然的短板——燃烧过程必然产生CO2。如果不对这个排放做处理，整个VPP的碳足迹就会很难看，在碳交易机制下会带来额外的成本。

碳捕集装置的价值在于把垃圾焚烧和系统“负碳化”连在一起。捕集下来的CO2可以直接封存或利用，而如果这个CO2刚好可以被P2G的甲烷化过程当成原料，就形成了一种很有意思的循环：垃圾焚烧产生CO2，碳捕集捕获CO2，P2G用捕集到的CO2加上电解水产生的氢气合成甲烷，再进天然气管网。整个系统不再只是单向的“发电—用电”，而是电、气、碳三网耦合。

1.2 垃圾焚烧加碳捕集的排放核算逻辑

我一开始做碳核算时最容易搞错的一点，是把垃圾焚烧排放直接当成燃煤电厂来算。垃圾焚烧的CO2排放因数其实分两部分：一部分来自化石碳（比如塑料、橡胶），另一部分来自生物质碳（比如厨余、纸张、木材）。在碳交易核算里，生物质碳部分通常被视为碳中和，不纳入净排放，或者按零排放处理。所以垃圾焚烧电厂的碳排放系数比燃煤机组低，但绝不是零。

加入碳捕集装置之后，净排放量的表达式可以写成：

净排放量 = 原始排放量 - 捕集量

如果捕集率足够高，净排放量可能变成负值。在碳交易机制中，负排放意味着你可以多出配额，转手卖给其他排放主体，这属于实实在在的经济收益。反而正是这个收益，才让碳捕集装置在经济上有了落地的空间。

从调度角度看，碳捕集并不是免费的“环保附件”。它的运行需要额外消耗能量——主要是再生塔所需的蒸汽或电力，这部分能耗必须从汽轮机发电出力里扣除。于是在模型里，VPP净出力不再是垃圾焚烧机组原始出力，而是原始出力减去碳捕集自身消耗后的差值。

1.3 电转气发挥作用的两个时机

P2G在调度模型里不是每时每刻都运行的，它主要出现在两个典型时段。第一个是风电光伏大发、电价低谷、系统弃电风险高的时段，此时如果VPP内部把多余的电直接卖给电网，价格很低，不如转成天然气储存或出售；第二个是碳捕集装置捕集到CO2但暂时无法封存或外运的时段，P2G的甲烷化装置反而能把它消纳掉，生成甲烷进燃气网。

P2G的加入，使得调度问题从纯电力平衡扩展到了电气耦合平衡。传统调度模型里，功率平衡约束只有电这一条；现在要多出一条天然气的产出与销售/储存平衡。而天然气的产出量和电力输入量之间的关系，本质上是一个转换效率约束。这就是“计及电转气协同”的核心含义。

2. 系统建模的核心细节

2.1 垃圾焚烧机组建模：先用线性化再谈优化

垃圾焚烧机组和常规火电机组在模型形式上很像，都需要处理出力上下限、爬坡速率、最小启停时间。但垃圾焚烧有一个额外特点：燃料不是可调度的煤或气，而是城市生活垃圾，来料基本是连续的，而且热值波动大。实际园区里通常有垃圾储坑，可以短时调整进料量，但不可能像火电厂那样按需采购燃料。所以我在建模时，把燃料输入近似为与实际出力分段线性相关的变量，而不是完全固定的常数。

用Matlab写约束时，我会把垃圾焚烧机组当作一个带0-1状态变量（启停状态）的常规机组来处理：

• 出力范围：P_waste_min <= P_waste(t) <= P_waste_max
• 爬坡约束：-RD <= P_waste(t) - P_waste(t-1) <= RU
• 最小运行/停运时间约束：用常规的MILP公式
• 启停逻辑约束：u_waste(t) - u_waste(t-1) = start(t) - stop(t)

这里最容易被忽略的是启动时的燃料消耗。垃圾焚烧锅炉冷启动会额外消耗天然气助燃，这部分成本虽然不大，但在24小时尺度下如果启停次数多，累积起来也挺可观。我在目标函数里单独加了启动成本项C_start * start(t)。

2.2 碳捕集装置建模：别把它的能耗想得太复杂

碳捕集装置的物理过程很复杂，涉及吸收塔、再生塔、溶剂循环，但站在调度层面向外看，它的外特性其实可以概括为“捕集多少碳，就消耗多少能量”。省事的做法是把它建模成一个线性能耗模型：

P_capt(t) = P_capt_base + k_capt * Q_capt(t)

其中Q_capt(t)是捕集到的CO2量，P_capt_base是维持装置运行的基本电耗，k_capt是单位捕集量的边际能耗。这个公式的好处是保持了整个优化问题的线性特性，可以直接在MILP框架里求解。

还需要加一个捕集能力上限：

Q_capt(t) <= Q_capt_max

同时，碳捕集装置的能耗会减少VPP的净出力，所以功率平衡约束里的“电源出力”要写成：

P_vpp_power(t) = P_waste(t) + P_wind(t) + P_PV(t) + P_storage(t) - P_capt(t)

写YALMIP约束时，P_capt(t)是一个连续变量，和Q_capt(t)一起作为决策变量放进优化问题。

如果想让模型更贴近实际，还可以用两段式线性函数来逼近碳捕集能耗的非线性曲线，即捕集率从零到某个阶段时，单位能耗较高；达到一定捕集量后，单位能耗趋于平稳。这样既保持线性，又比单一斜率更真实。我最终的算例模型就采用了这种分段处理。

2.3 P2G设备的运行约束

P2G设备本质上是一个“电转气”转换器，输入是电能，输出是天然气（或氢气）。在建模时，把它的工作范围分成三段是有必要的：第一段是待机状态，电耗为0，产气为0；第二段是正常运行状态，输入电功率下限到上限之间，产气量与输入电量近似线性；第三段是满负荷状态，输出达到最大值。

数学表达上，正常段的转换关系很简单：

G_p2g(t) = eta_p2g * P_p2g(t)

其中eta_p2g是电转气综合效率。注意这里的效率不是单一电解槽效率，还包括了压缩、甲烷化等环节，整体通常取50%-60%。如果只考虑制氢，效率可以到70%-80%，但考虑到后续要合成甲烷，我按55%来取。

约束条件至少包括：

• 0 <= P_p2g(t) <= P_p2g_max
• G_p2g_min <= G_p2g(t) <= G_p2g_max
• 天然气产量需要满足当地气网销售或储气罐约束

如果模型里有储气罐，还要加入储气量平衡方程：

S_gas(t+1) = S_gas(t) + G_p2g(t) - G_sell(t)

这个方程在Matlab里表现为差分约束，YALMIP可以直接支持带时移的变量表达式，写起来很顺手。

我建议把P2G的最小技术出力设成非零值。因为现实中电解槽在极低负荷下运行效率很差，而且频繁启停会缩短寿命。模型里如果允许P2G在0到满负荷之间任意取值，最优解往往会利用这个“连续可调”的便利，给出一些不符合物理实际的调度方案。

2.4 系统功率平衡与电网交互约束

VPP整体与上级电网的交互，用联络线约束来刻画。我通常会设置一个最大购电功率和最大售电功率，同时为了避免“又买又卖”这种明显背离交易逻辑的结果，还会加一个购售电互斥约束。

功率平衡约束是每个时刻必须满足的硬约束：

P_waste(t) + P_wind(t) + P_PV(t) + P_disch(t) + P_buy(t) = P_load(t) + P_charge(t) + P_capt(t) + P_p2g(t) + P_sell(t)

把等式右侧的P_capt和P_p2g都当作负荷来看待，碳捕集装置、P2G设备就相当于两个“柔性负荷”，这也是它们能参与优化调度的本质原因。如果题目里还需要计及旋转备用，可再加一条：

P_reserve(t) >= Reserve_req(t)

备用约束我在算例里没有细究，但实际工程中它是必须的，尤其当风光渗透率较高时，备用不足会直接影响调度方案的可行性。

3. Matlab实现与求解器配置

3.1 环境搭建和数据准备

我用的是Matlab + YALMIP + Cplex的组合。YALMIP是一个建模层的工具箱，它对变量类型和约束条件的表达非常灵活，尤其适合快速搭建MILP模型；Cplex负责底层求解。如果你没有Cplex授权，也可以换成Gurobi，或者干脆用Matlab自带的intlinprog，但YALMIP的建模效率能帮你省下大量手写稀疏矩阵的时间。

开始写代码前的第一件事是统一单位。我习惯把功率单位统一为MW，电量单位统一为MWh，时间尺度用小时，电价单位用元/MWh。为什么强调这个？因为一旦功率和电量混用、万和亿混用，优化求解器会很容易出现数值病态，表现为“明明模型没问题，但求解器报infeasible”或者迭代不动。

数据准备部分，典型代码如下：

% 系统数据 T = 24; % 调度周期，小时 load_data = data.Load; % 1x24 向量，单位MW wind_data = data.Wind; % 1x24 向量，单位MW pv_data = data.PV; % 1x24 向量，单位MW price_buy = data.PriceBuy; % 购电价，元/MWh price_sell = data.PriceSell; % 售电价，元/MWh % 机组参数 waste_min = 10; % MW waste_max = 50; % MW RU_waste = 15; % MW/h RD_waste = 15; capt_base = 1.5; % 碳捕集基础电耗 capt_coef = 0.08; % MWh/tCO2 capt_max = 120; % t/h eta_p2g = 0.55; % 电转气效率 p2g_max = 15; % MW

3.2 YALMIP变量定义与目标函数

变量定义的核心原则是：一个变量对应一个物理量，不要为了省代码把多个不同物理含义的东西合并到一个变量里。

% 决策变量 P_waste = sdpvar(T, 1); % 垃圾焚烧发电出力，MW u_waste = binvar(T, 1); % 垃圾焚烧启停状态，0-1 v_start = binvar(T, 1); % 启动动作 v_stop = binvar(T, 1); % 停机动作 P_capt = sdpvar(T, 1); % 碳捕集装置总电耗，MW Q_capt = sdpvar(T, 1); % 捕集CO2量，t/h G_p2g = sdpvar(T, 1); % P2G产气量，MW(以天然气热值计) G_sell = sdpvar(T, 1); % 对外出售天然气量，MW P_buy = sdpvar(T, 1); % 从上级电网购电，MW P_sell = sdpvar(T, 1); % 向电网售电，MW S_gas = sdpvar(T+1, 1); % 储气罐储量，MW

目标函数里，我把碳交易成本放在最后一项。注意碳交易成本和碳捕集的能耗成本是两回事，前者和净碳排放量挂钩，后者已经体现在P_capt对功率平衡的影响里了。如果调研报告给出的碳价是price_carbon，配额数据是预先给定的，那么碳交易成本可以写成：

% 目标函数：总运行成本最小 cost_fuel = c_fuel * P_waste; % 垃圾燃料成本 cost_start = c_start * v_start; % 启动成本 cost_buy = price_buy' * P_buy; % 购电成本 profit_sell= price_sell' * P_sell; % 售电收益 carbon_emission = waste_emission_factor .* P_waste; % 原始排放 carbon_net = carbon_emission - Q_capt; % 净排放 cost_carbon = price_carbon * (carbon_net - carbon_quota_alloc); Objective = sum(cost_fuel) + sum(cost_start) + sum(cost_buy) ... - sum(profit_sell) + sum(cost_carbon);

如果配额分配量大于净排放，cost_carbon就会是负值，相当于收益。这个信号在优化里会自动引导系统多开碳捕集装置，因为每捕集一吨CO2，就是真金白银的收益。

3.3 约束条件写入与求解

约束条件写起来最关键的一点，是不要忘记给0-1变量加逻辑关系：u_waste(t) - u_waste(t-1) = v_start(t) - v_stop(t)

这个约束把所有开关状态和启停动作串起来。如果你漏掉它，会出现一种诡异的结果：机组状态变量在连续几个时间段都为1，但同时出现启动和停机动作都等于1的荒谬解。

接下来是YALMIP约束的核心写法：

Constraints = []; % 垃圾焚烧出力与爬坡 Constraints = [Constraints, waste_min .* u_waste <= P_waste <= waste_max .* u_waste]; for t = 2:T Constraints = [Constraints, -RD_waste <= P_waste(t)-P_waste(t-1) <= RU_waste]; end % 启停逻辑 for t = 2:T Constraints = [Constraints, u_waste(t) - u_waste(t-1) == v_start(t) - v_stop(t)]; end % 碳捕集变量关系 Constraints = [Constraints, P_capt == capt_base + capt_coef * Q_capt]; Constraints = [Constraints, 0 <= Q_capt <= capt_max]; % P2G转换关系 Constraints = [Constraints, G_p2g == eta_p2g * P_p2g]; Constraints = [Constraints, 0 <= P_p2g <= p2g_max]; % 储气罐动态 Constraints = [Constraints, S_gas(t+1) == S_gas(t) + G_p2g - G_sell]; Constraints = [Constraints, 0 <= S_gas <= S_gas_max]; % 功率平衡 Constraints = [Constraints, P_waste + wind_data + pv_data + P_disch ... + P_buy == load_data + P_charge + P_capt + P_p2g + P_sell]; % 求解 ops = sdpsettings('solver', 'cplex', 'verbose', 2); result = optimize(Constraints, Objective, ops);

这一段代码看似短，但里面其实隐含了一个细节：P_waste和u_waste相乘，在YALMIP里会被自动处理成线性化的大M约束。如果你不用YALMIP而是自己手写线性化，就要手动引入大M和辅助变量。对大M的取值，我的建议是取机组最大出力的1.1倍，太大太小都会出问题，取值过大容易让求解器数值不稳定，取值过小又会错误地削减可行域。

3.4 结果后处理：画图前的几个检查

求解完成后，不是直接plot就完事了。我一般会先做三项检查。

第一项是检查求解状态。result.problem必须等于0，表示求解成功。如果等于1，说明不可行；如果等于2，说明遇到数值问题；如果是其他值，要一个个查清楚，不要稀里糊涂拿着一个不可行的解去画图。

第二项是核对功率平衡是否真正满足。我会把求解得到的变量代回到平衡约束左边减右边，逐时差值的绝对值应当小于一个很小的容差，比如1e-4。如果偏差明显，说明数值精度不够，可能需要降低求解器内部的tolint或tolfcast阈值。

第三项是检查碳捕集量和P2G产气量是不是被“过度使用”。有时候优化解会把P2G开满，但产出的天然气却卖不出去，因为储气罐满了。这种情况通常源于约束不完整，比如没有设置储气罐上限，或者没有设置天然气销售价格下限。排查顺序是：先看储气罐储量曲线，再看天然气的卖出价格是否合理，最后检查气体平衡方程。

画图时，我推荐至少画三张图：

• 第一张：各电源出力与负荷曲线叠加图，包括垃圾焚烧、风电、光伏、储能、购电，能直观看出系统如何平衡负荷。
• 第二张：碳捕集量与P2G运行功率对比图，这两个变量是把握电气碳协同的关键。
• 第三张：碳交易成本累计曲线，能看出哪个时段因为碳排放增加了成本、哪个时段因为捕集在盈利。

保存图片时记得用：

figure(1); plot(1:T, P_waste_value, 'LineWidth', 1.6); grid on; xlabel('时间/h'); ylabel('功率/MW'); title('垃圾焚烧机组出力'); saveas(gcf, 'result_waste.png');

如果你在写论文时还要导数据，建议用writematrix(P_waste_value, 'P_waste.xlsx')直接保存成Excel，方便后续做敏感性分析。直接用save P_waste.mat也可以，但Excel在跨软件交流时更友好。

4. 实测中踩过的坑与解决办法

4.1 一起手就是Infeasible：先检查约束之间的矛盾

我最开始搭出来的模型，十个案例里有四五个直接报infeasible。后来排查下来，发现大多是类似这种问题：垃圾焚烧机组的最小出力是10MW，但某些时段风电+光伏加上最小购电功率已经超过负荷，而且储能又恰好被填满，导致多余的电无处可去。这种时候系统当然无解。

解决办法是重新审视功率平衡约束里各物理量的取值范围。最直接的做法是给P2G设备更高的最大输入功率，或者降低储能约束的上下限，甚至允许弃风弃光。加弃风弃光变量本身就是一种经典处理方式：

P_curtail = sdpvar(T, 1); Constraints = [Constraints, 0 <= P_curtail <= wind_data]; Constraints = [Constraints, P_wind_used == wind_data - P_curtail];

在实际项目中，弃风弃光并不是丢人的事情，它本身就是调度决策的一部分，只是在目标函数里通常要用一个很小的惩罚系数来避免不必要的弃用。

另一个常见原因是负荷数据选取出了问题。有些公开数据集里的负荷是“全年最大负荷”，而算例里用的系统容量对应的是“典型日负荷”，两者数量级差好几倍。这种错误非常隐蔽，但后果很直接——最优解会拼命售电，反而导致购电售电互斥约束失效。

4.2 0-1变量太多导致求解过慢

如果你的系统里有多台机组、多套碳捕集装置、多套P2G，每个时间段都要引入一批0-1变量，那么24小时尺度上MILP的求解规模会指数级增长。我遇到过一个包含5台机组、3套储能、4套P2G的场景，YALMIP生成的约束矩阵直接让Cplex跑了将近二十分钟还没收敛。

解决思路有几个。第一，尽量合并相同特性的机组，用“聚合机组”代替独立机组，把5台同类垃圾焚烧机组聚合成一台容量5倍的大机组。第二，对P2G和碳捕集装置这类运行特性简单的对象，可以先去掉它们的启停0-1变量，只保留连续运行范围，看优化结果是否显著变化；如果变化不大，说明0-1变量在这个场景下是多余的。第三，设置合理的MIP间隙：

ops = sdpsettings('solver', 'cplex', 'cplex.mip.tolerances.mipgap', 0.01);

牺牲0.1%-1%的精度，能换回往往几十倍的提速。工程应用中，1%的MIP间隙带来的误差完全在接受范围内。

4.3 同一个模型换一批数据，结果就不靠谱了

这种问题多半出在参数的量纲和数量级上。比如碳捕集的能耗系数是0.08 MWh/tCO2，而垃圾焚烧的电价成本是几百元/MWh，碳价可能是几十元/吨，三个量的数量级差10倍以上。求解器内部处理这种混合尺度问题，很容易在数值稳定性上翻车。

我的习惯是给碳价、电价、煤耗费用统一折算到“元/MWh”或“元/吨”的标准单位，并检查每个变量的典型取值都在1e-2到1e+4的范围内。如果某个参数过大或过小，就要考虑归一化。YALMIP里有rescale命令，可以自动对约束矩阵做尺度变换，遇到反复报数值问题的时候可以试试。

4.4 碳交易价格设定影响结果的方式

如果你把碳价设得特别高，优化器自然倾向于把碳捕集开满，甚至可能出现捕集量大于产生量的荒谬情况。这不是求解器错了，而是你忘了加一个物理上限约束：捕集量不能超过烟气里的CO2总量。

Q_capt(t) <= Q_capt_max(t)只约束了捕集装置本身的能力上下限，并没有约束它与垃圾焚烧原始排放之间的逻辑关系。更严谨的写法是：

Constraints = [Constraints, Q_capt <= waste_emission_factor .* P_waste];

这条约束的意义在于，碳捕集装置再怎么捕集，也不能捕集到垃圾焚烧机组尚未产生的CO2。这在实际工程中很好理解，但在模型里经常被漏掉。

一旦漏掉，高碳价场景下求解器就会把Q_capt顶到捕集上限，造成净排放量严重偏低，碳交易收益虚增，最终给出一套在现实中根本不可能实现的调度方案。

4.5 表格速查：常见报错与处理方案

这里再补充一个比较隐蔽的坑：YALMIP里的binvar变量在约束里乘以连续变量之后，如果系数矩阵过大，有些免费求解器会直接拒绝求解。如果只能用免费求解器，建议改用intlinprog手动建模，或者把问题拆小。付费的Cplex/Gurobi对这类线性化约束的预处理能力更强，大M约束写进去几乎没有压力。

5. 扩展方向与个人体会

5.1 怎么把这段代码迁移到其他场景

这套模型最大的价值不在那几行代码本身，而在它的结构是可替换的。如果你想把垃圾焚烧换成生物质发电、把碳捕集换成其他负碳技术、把P2G换成储能，大部分约束框架都可以直接套用，改的只是参数和个别变量之间的耦合关系。

比如，如果场景变为“含电转气与储氢的虚拟电厂”，那就要把P2G的产物从天然气改成氢气，同时增加储氢罐的储量方程；如果场景变为“含碳捕集的燃气机组参与调峰”，那就要把垃圾焚烧机组的燃料成本曲线换成天然气的梯级气价。不要推倒重来，先改参数、再改耦合项，复杂度会低很多。

另外，模型的时间尺度也可以从24小时扩展到168小时，也就是一周。只需要把数据维度从24改成168，同时补上储能、储气罐跨日约束即可。不过要注意，调度周期变长之后，垃圾焚烧机组的最小起停时间约束会明显增加MILP求解难度，建议先用MIPGap手段压住求解时间。

5.2 几句实在话

回头来看，这个项目最难的部分其实不是数学建模，也不是Matlab代码本身，而是“你能不能用一个线性模型把真实物理过程描述到足够准确”。碳捕集装置、P2G设备、垃圾焚烧机组，每一项单拎出来都有复杂的非线性特性，但工程调度模型要的是“可求解”和“够准确”之间的平衡。

我一直坚持一个做法：先跑通最简单的小算例，确认约束框架和求解流程没问题，再一步步加入复杂性。不要一上来就把整篇论文里的所有设备、所有约束都塞进模型，那样你根本不知道问题是出在数据上、建模上还是求解器上。

最后再分享一个小技巧：做敏感性分析的时候，把碳价从20元/吨一直扫到200元/吨，画出“碳价—系统总成本—P2G运行小时数”的关系曲线，往往一眼就能看出系统在什么价格区间会发生“质变”——比如P2G从偶尔运行变成全天运行，碳捕集从部分捕集变成满捕集。这条曲线是写论文、做方案汇报时最有说服力的材料，比贴一大堆公式好得多。