虚拟电厂电-气-碳联合优化调度的Matlab建模实现

这两年做低碳电力系统相关的网格优化，最常被问到的一个问题就是：虚拟电厂调度里同时塞进电转气、碳捕集和垃圾焚烧，到底该怎么建模？单看任何一个设备，网上的资料都不算少，但要在一套调度模型里把“电、气、碳”三条线串起来，很多朋友一开始会发懵。这篇文章我不打算写那种“从原理到展望”的综述，而是直接以一个我自己用Matlab实现的虚拟电厂优化调度项目为底子，把变量怎么定、约束怎么写、求解器怎么调、算例怎么安排，以及我实际踩过的坑，按顺序摊开来讲。适合正在做综合能源系统优化、低碳调度、碳捕集与封存建模的研究生，以及想用Matlab把调度模型落到代码里的工程师。

1. 项目背景：为什么要把电转气、碳捕集和垃圾焚烧放在同一个算盘里

1.1 虚拟电厂不再只是“风光储”的拼盘

传统的虚拟电厂调度，大概就是聚合分布式光伏、风电、储能和备用柴油机组，以最小运行成本为目标去优化各设备的出力。这种模型已经很成熟，核心约束也只有功率平衡、机组出力上下限、储能SOC等。可一旦加入垃圾焚烧机组，问题立刻不一样。

垃圾焚烧机组的特点是“不能想停就停”。焚烧炉需要维持稳定燃烧，垃圾仓里的垃圾每天都要处理，城市固废的连续清运要求决定了焚烧机组必须连续运行。在夜间的低负荷时段，其他机组可以大幅降出力甚至停运，垃圾焚烧机组往往只能压到最低稳燃负荷，没办法继续压。这个“刚性出力”如果不建模，调度方案在工程上根本执行不了。所以在模型里，垃圾焚烧机组不是普通燃煤机组，它更像是一个“必须运行但出力范围受限”的特殊节点。

再加上碳捕集设备，系统里的电气关系又多了一层：碳捕集本身要耗电，捕集量越大，虚拟电厂提供给电网的净出力就越低。如果目标函数里只算机组发电收益，不把碳捕集能耗和碳排放成本放进去，优化器会把碳捕集率直接拉到零，看起来“电量更多、收益更高”，实际上碳排放超了，碳成本一算可能完全得不偿失。

1.2 电转气和碳捕集形成的碳-气闭环

电转气（Power to Gas，P2G）简单说就是用电力电解水制氢，再和二氧化碳甲烷化生成合成天然气。这一过程最吸引人的地方在于：它能把暂时用不掉的再生能源电力转化成气体储存起来，解决风光出力的波动问题。而碳捕集（CCS）捕集下来的高浓度二氧化碳，正是甲烷化反应所需要的原料。两个设备放在同一座虚拟电厂里，天然形成互补关系。

从调度角度看，P2G是一个灵活的可调负荷：电价低、风光大发的时候多用电产气；用电高峰、电价高的时候减少用电，甚至把产出的天然气送去燃气轮机发电或直接外售。碳捕集设备承担了碳减排的任务，捕集的CO2一部分存储或外送，一部分可以供给P2G做甲烷化原料。再加上垃圾焚烧烟气工况相对稳定，二氧化碳浓度不低，为碳捕集提供了相对好的接入条件。三者放在同一个优化问题中，表面上是成本优化，实际上是在做“能量在各个时间断面之间的转移”和“碳在哪一个环节被固定”的路径选择。

1.3 这套虚拟电厂优化调度代码的定位

我想先明确一下，这不是一套给人看“理论推导”的代码，而是一套能跑出结果的调度模型。整体采用确定性日前调度框架，以一天24小时为调度周期，时间尺度是1小时。系统包含垃圾焚烧机组、常规燃气机组、光伏、风电、碳捕集装置、电转气装置以及储气设施。模型最终输出每台机组各时刻出力、启停状态、碳捕集率、电转气功率等决策量，并给出总运行成本和碳排放水平。

代码实现我建议用MATLAB + YALMIP + 商用求解器（CPLEX或Gurobi）的组合。YALMIP的建模语法非常简单，熟悉MATLAB的人很容易上手；底层用MILP求解器处理二进制启停变量，算例规模不大时几分钟就能出结果。如果不想依赖YALMIP，直接用MATLAB自带的intlinprog也可以，只是变量和约束多了之后，代码会变得比较难维护。下面所有讲解都按“能直接复现”的标准来展开。

2. 调度模型怎么搭：目标函数、决策变量和核心约束

2.1 目标函数：不是单纯“发电成本最低”

很多初学者会把目标函数写成“系统运行成本最小”，然后只加燃料成本。实际做低碳调度时，必须把碳排放的外部性也内部化。我的模型目标函数如下：

min J = Σt Σk [ f_k(P_k,t) * u_k,t + C_{su,k} * v_k,t ] + Σt [ c_p2g * P_p2g,t + c_ccs * E_cap,t ] + Σt [ c_cut * (P_re_avail,t - P_re,t) ] + Σt [ λ_co2 * (E_emit,t - E_cap,t - E_quota,t) ]

逐项解释：

• 第一项是发电成本。f_k(P_k,t)可以简化为二次函数aP² + bP + c，也可以分段线性化；u_k,t是机组启停状态；v_k,t是启动变量，用于计启动成本。
• 第二项是电转气运行成本和碳捕集设备运行成本。电转气成本主要来自电解槽的损耗、水耗和维护；碳捕集成本可以按捕集的CO2吨数折算，也可以按再生塔能耗折算。
• 第三项是弃风弃光惩罚。可再生能源有预测出力序列，如果调度不能全额消纳，就要在目标函数里给一个惩罚，否则优化器为了照顾机组调节能力，很可能主动弃掉便宜的风光。
• 第四项是碳成本。E_emit,t是虚拟电厂实际排放量，E_cap,t是当时间接或直接固定的CO2量，E_quota,t是免费碳配额。如果实际排放超过配额，需要按碳价购买；反之可以出售盈余配额。碳价本身按场景参数设置。

为什么这样设计？因为如果目标函数里没有碳价，那么碳捕集设备和P2G的“碳固定”价值就无法体现。优化器宁可排放CO2，也不愿意花额外电耗去运行碳捕集装置。加入碳价之后，模型才会在“用电捕碳增加的电费”和“减少购买碳配额获得的收益”之间做权衡。

2.2 决策变量：一张表说清楚

整个模型的决策变量可以分成三类：机组状态类、能量流类和碳流类。

注意，机组启停状态必须是二进制变量，这就决定了整个模型是一个混合整数规划（MILP）或混合整数二次规划（MIQP）。如果你有条件用二次目标，CPLEX/Gurobi可以直接求解MIQP；如果担心求解时间，可以先把机组成本函数线性化为分段函数，变成MILP之后求解稳定性更高。

2.3 核心约束：机组、碳捕集和电转气

机组模型部分，最基本的是出力上下限约束和爬坡约束：

Pmin(k) * u(k,t) <= P(k,t) <= Pmax(k) * u(k,t)

这个约束表达的是：只有当机组运行时u=1，出力才能落在上下限内；机组停运时u=0，出力强制为0。

爬坡约束要区分向上爬坡和向下爬坡，用上一时刻的出力做基准：

P(k,t) - P(k,t-1) <= RampUp(k) * (1 - v(k,t)) + Pmax(k) * v(k,t) P(k,t-1) - P(k,t) <= RampDown(k)

启动过程允许一定的宽限，不然机组刚启动那一时刻根本爬不了那么多。这一步是我调参时经常被卡住的地方。

对于垃圾焚烧机组，还要额外加“连续运行”约束。最简单的做法是给一个最小运行时间MINUP和最小停运时间MINDOWN：

Σ_{i=t-MINUP+1}^{t} u(i) >= MINUP * (u(t) - u(t-1)) Σ_{i=t-MINDOWN+1}^{t} (1-u(i)) >= MINDOWN * (u(t-1) - u(t))

这个模型在实际代码里可以简化，不一定每个机组都强制建模，但垃圾焚烧机组建议至少设一个MINUP=24，即全天必须保持并网运行。这样最符合垃圾焚烧电厂连续运行的真实特性。

碳捕集部分的建模要抓住两个关键量：烟气中的CO2总量和捕集能耗。设垃圾焚烧机组发电量为P_rdf(t)，单位发电量对应的原始CO2排放因子为e0，则实际可捕集的CO2量为：

E_raw(t) = e0 * P_rdf(t)

捕集量等于捕集率乘以原始排放量：

E_cap(t) = alpha(t) * E_raw(t)

碳捕集装置的能耗可以简化为单位捕集能耗beta乘以捕集量：

P_ccs(t) = beta * E_cap(t)

这一项必须进入电功率平衡方程。如果丢掉了这一步，模型给出的“零成本捕碳”方案就是空中楼阁。

电转气部分同样要分成电和气两条线。电线上，P2G是一个负荷；气线上，其产气量与输入电功率之间用转换效率eta_p2g连接：

G_p2g(t) = eta_p2g * P_p2g(t)

储气罐的连续状态约束为：

G_sto(t+1) = G_sto(t) + G_p2g(t) - G_fuel(t) - G_sell(t)

其中G_fuel(t)是燃气机组消耗的气量，G_sell(t)是外售气量。储气罐有容量上下限和注入/采出速率限制，这和电池储能SOC的建模思路完全一样。

2.4 电、气、碳三条平衡如何闭环

系统最终必须满足电功率平衡。所有电源出力加上从外部购电，等于所有负荷：

Σk P(k,t) + P_re(t) + P_grid_buy(t) + P_dis(t) = P_load(t) + P_p2g(t) + P_ccs(t) + P_ch(t)

其中P_grid_buy是虚拟电厂与外网交换的功率，允许购电和售电两种状态。垃圾焚烧机组、燃气机组本身的电力必须进入平衡，碳捕集和P2G的消耗也要出现在右侧。

气平衡则是：P2G产气加上外部购气，等于燃气机组消耗加外售气量。碳平衡则简化成一个等式：

E_net(t) = E_emit(t) - E_cap(t)

这个净排放量就是参与碳交易计算的基准。捕集下来的CO2一部分直接存储，另一部分送去P2G甲烷化，这部分在实际模型里可以作为P2G产气量的一个原料约束，比如甲烷化所需要的CO2量与产气量成正比，可以写成：

G_p2g(t) <= mu * E_cap(t)

这里mu是单位产气量对应的CO2消耗系数。这个约束把“碳”和“气”真正绑在一起，模型会知道P2G不是无中生有，而是要消耗碳捕集资源。

3. Matlab代码实现：从空模型到可复现算例

3.1 代码结构：先把文件分清楚

写这种优化模型，最忌讳的就是把所有东西堆在同一个main.m里。我的习惯是把代码拆成四块：

• main.m：主程序，负责设置参数、调用模型、输出结果。
• init_data.m：填负荷数据、风光预测、机组参数、碳价、效率系数等。
• build_model.m：用YALMIP定义变量、写约束和目标函数。
• plot_result.m：绘制各机组出力曲线、碳捕集率曲线、电转气功率曲线、功率平衡图。

这样后续做场景对比时，只需要改init_data.m，不用反复动模型主体。

3.2 用YALMIP定义决策变量

直接进入代码。下面这段是一个可以运行的模型骨架。首先是初始化：

%% 初始化 T = 24; % 24小时 G = 3; % 机组数量：1垃圾焚烧 2燃气 3备用 Pmin = [40; 10; 0]; % MW Pmax = [120; 60; 30]; Ramp = [30; 30; 20]; % 机组成本系数 aP^2+bP+c a = [0.01; 0.015; 0.02]; b = [12; 14; 18]; c = [50; 30; 20]; % 污染物/碳相关 e0 = 0.45; % t/MWh 垃圾焚烧原始CO2排放系数 beta = 0.2; % MWh/tCO2 捕集单位CO2对应的电耗 eta_p2g = 0.6; % P2G综合效率 quota = 0.05; % tCO2/MWh 免费配额折算系数 lam_co2 = 30; % 元/tCO2

然后是YALMIP变量：

%% 定义YALMIP变量 u = binvar(T, G, 'full'); % 机组启停 v = binvar(T, G, 'full'); % 启动动作 P = sdpvar(T, G, 'full'); % 机组出力 P_re = sdpvar(T, 1); % 可再生实际出力 P_re_avail = ... % 风光预测数据(读入) alpha = sdpvar(T, 1); % 捕集率 Pp2g = sdpvar(T, 1); % P2G功率 Gp2g = sdpvar(T, 1); % P2G产气量 Gsto = sdpvar(T, 1); % 储气罐SOC Gbuy = sdpvar(T, 1); % 购气量 Pgrid = sdpvar(T, 1); % 外购电，正数购电，负数售电 P_ccs = sdpvar(T, 1); % 碳捕集电耗 Ecap = sdpvar(T, 1); % 碳捕集量 Eemit = sdpvar(T, 1); % 实际排放量

定义一个变量集合便于后续循环添加约束。用YALMIP的好处是约束可以直接用>=、<=这样的运算符写。

3.3 核心约束怎么写进模型

下面把机组约束写成一个循环：

%% 约束集合 Ccon = []; for k = 1:G % 出力上下限 Ccon = [Ccon, Pmin(k)*u(:,k) <= P(:,k) <= Pmax(k)*u(:,k)]; % 启停逻辑 for t = 2:T Ccon = [Ccon, v(t,k) >= u(t,k) - u(t-1,k)]; end Ccon = [Ccon, v(:,k) >= 0]; % 爬坡简化：必要时用大M法 for t = 2:T Ccon = [Ccon, P(t,k) - P(t-1,k) <= Ramp(k)*(1-v(t,k)) + Pmax(k)*v(t,k)]; Ccon = [Ccon, P(t-1,k) - P(t,k) <= Ramp(k)]; end end

这里用了一个比较粗暴的大M启动爬坡处理，实际工程中还可以更精细，比如把启动爬坡和正常运行爬坡分开。但这已经足够跑通一个调度算例。

然后是碳捕集和电转气约束：

%% 碳捕集关系 % 原始CO2排放：只计垃圾焚烧部分 Eraw = e0 * P(:,1); Ccon = [Ccon, Ecap == alpha .* Eraw]; Ccon = [Ccon, P_ccs == beta * Ecap]; Ccon = [Ccon, 0 <= alpha <= 1]; % 实际净排放 Ccon = [Ccon, Eemit == Eraw - Ecap]; %% P2G和储气 Ccon = [Ccon, Gp2g == eta_p2g * Pp2g]; Ccon = [Ccon, 0 <= Pp2g <= 60]; % P2G容量限制 Ccon = [Ccon, 0 <= Gsto <= 80]; % 储气罐容量MWh % 储能连续方程 for t = 2:T Ccon = [Ccon, Gsto(t) == Gsto(t-1) + Gp2g(t) - 5]; % 5假设燃气消耗 end Ccon = [Ccon, Gsto(1) == 20];

储气这一段我故意把燃气消耗写成了常数5，真实模型中应该用燃气机组出力换算。这里作为骨架简化即可。

3.4 目标函数和求解设置

目标函数设置如下：

%% 目标函数 Objective = 0; for k = 1:G for t = 1:T Objective = Objective + (a(k)*P(t,k)^2 + b(k)*P(t,k) + c(k)*u(t,k)); end end Objective = Objective + 0.5*sum(Pp2g) + 5*sum(Ecap) + 100*sum(P_re_avail - P_re); Objective = Objective + lam_co2 * sum(Eemit - quota * sum(P,2));

然后求解：

%% 求解 ops = sdpsettings('solver','cplex','verbose',2,'usex0',1); sol = optimize(Ccon, Objective, ops); if sol.problem == 0 P_opt = value(P); Ecap_opt = value(Ecap); Pp2g_opt = value(Pp2g); alpha_opt = value(alpha); else disp('求解失败：' + sol.info); end

这里sol.problem == 0来自YALMIP的返回状态，0表示最优解。如果模型无解，需要回到约束里排查。

实际项目中，目标函数的二次项a(k)*P(t,k)^2会让模型变成MIQP。如果数据量大导致MIQP求解慢，可以提前把二次成本分段线性化。比如把每台机组的出力范围切成三段，每段的边际成本近似为常数，然后引入分段权重变量。这样可以保证大型算例的求解稳定性。

4. 典型算例与结果分析

4.1 场景设置：有碳捕集P2G和无碳捕集P2G

为了看出这套模型的价值，我通常会做两个场景对比：

• 场景A：基准场景，不含碳捕集、不含电转气，只有常规机组和垃圾焚烧。
• 场景B：完整场景，含碳捕集、电转气、储气罐和碳市场约束。

两个场景用同一天的负荷和风光预测数据。具体参数是：垃圾焚烧机组额定容量120 MW，最低出力40 MW；燃气机组60 MW，备用机组30 MW；风电预测容量约80 MW，光伏约50 MW；日最大负荷220 MW。

4.2 核心指标的变化

在这样一个算例上，典型结果大致是：

为什么会是这样？因为P2G在夜间风电大发、负荷低谷时启动，把原本要被弃掉的6到8个小时风电变成天然气储存起来，在早高峰时给燃气机组用，相当于把弃电“搬了个家”。碳捕集则是在白天垃圾焚烧机组出力较高的时段多捕碳，利用捕集后较低的净排放水平减少碳配额购买。虽然碳捕集和P2G本身都消耗电能，但在算例的总成本账上，省下的碳费用和减少的弃电惩罚超过了它们的运行成本，所以整体经济性反而更好。

4.3 结果里值得留意的现象

第一个现象是：垃圾焚烧机组在夜间出力始终压在最低技术出力40 MW附近，而燃气机组和备用机组完全停运。这很符合焚烧炉连续运行的物理特点，也是模型里MINUP=24约束起作用的结果。如果把这个约束删掉，优化器可能会给出夜间停运垃圾焚烧机组、白天再启动的“数学最优”方案，但现实中根本没法执行。

第二个现象是：碳捕集率并不像很多人想的“越高越好”。即便碳价给到了60元/吨，捕集率也没有无限逼近1，而是停在一个相对合理的水平。原因很简单：捕集率升高会让碳捕集电耗大幅上升，进而占用电负荷空间，甚至导致系统需要用高价购电来维持平衡。模型在碳成本和电成本之间自动找平衡点，这正是优化调度的意义。

5. 实际调参中踩过的坑与排查方法

5.1 模型怎么突然无解了

这是我遇到最多的问题。无解往往不是因为目标函数写错，而是约束之间出现了物理上不可能的矛盾。最常见的有三种：

还有一个隐蔽问题：垃圾焚烧机组的最低出力40 MW，而风光大发时的总负荷只有100 MW，再加上P2G和CCS消耗也无法把120 MW的净负荷顶住时，系统只能被迫弃风。很多初学者这时候把弃风量作为硬性约束要求“必须全额消纳”，结果无解。正确做法是：弃风不做约束，而作为目标函数里的惩罚项，让模型自己去权衡。

5.2 求解太慢怎么办

如果算例只有24个时段和3台机组，其实很快。一旦把时间尺度拉长到8760个小时，或者机组数量增加到10台以上，MIQP的求解时间会指数级上涨。我常用的处理办法有三个。

第一，把二次目标函数分段线性化，把MIQP变成MILP。CPLEX和Gurobi对线性MIP的求解效率远高于二次MIP。第二，设置MIP gap。sdpsettings('cplex.mip.tolerances.mipgap', 0.01)表示允许1%的次优偏差，对工程调度完全够用，速度能快好几倍。第三，使用热启动。先用简化模型算一组结果，把机组启停状态传给完整模型作为初始解，能显著减少分支定界树的搜索范围。

5.3 碳排放结果出现负值

有些朋友在结果里看到“净排放量为负”，觉得很兴奋，认为是实现了负碳。但大多数情况下不是拐点，而是模型漏了约束。负排放往往是因为碳捕集量Ecap的来源是垃圾焚烧产生的CO2，但实际焚烧环节可能已经将部分CO2视为“生物源碳排放”，在碳核算里常被豁免。模型中如果不区分化石源碳和生物源碳，就会发生“一边捕集本来就不用付费的生物源碳、一边卖出碳配额”的不合理行为。

我的建议是：在模型里单独设置一个化石源碳排放比例，假设垃圾焚烧总CO2排放中只有r_fossil比例来自化石成分，比如塑料、橡胶等，剩余部分视为生物源。碳成本计算时只对化石源部分收碳价，碳捕集收益也只针对化石源部分。否则结果会很漂亮，但经不起碳核查。

6. 后续还能怎么扩展：让调度模型更贴合真实挑战

6.1 从确定性调度走向不确定性优化

上面的模型用的是确定性日前预测数据，实际上风电、光伏、负荷都有预测误差。进一步扩展时，可以用场景法或者鲁棒优化把所有不确定量纳入模型。最简单的场景法就是生成多组风光预测场景，然后在目标函数中对所有场景求期望成本，同时保证每个场景都能满足约束。此时P2G和储气的价值会更明显：因为储气罐可以在多个场景之间发挥平衡作用，比只针对一条预测曲线更贴近真实运行。

6.2 把碳捕集的时序特性做得更细

化工过程是有惯性的，碳捕集装置从低负荷切到高负荷，并不是瞬间完成，而是需要分钟到小时级的过渡时间。如果调度周期是1小时，这还可以接受；如果要做分钟级滚动调度，就要给碳捕集设备加上爬坡约束和最小连续运行时间约束。同理，电解槽启动次数和功率波动范围也应当受寿命限制，否则模型会让电解槽频繁启停，实际设备根本受不了。

6.3 一点个人体会

这套代码做完之后，我最大的感受是：调度模型的价值不在于把目标函数写得多复杂，而在于每个约束是否对上了设备真实的物理特性。垃圾焚烧连续运行、碳捕集耗电、P2G产气储气、燃气机组消耗气体，它们之间环环相扣，漏掉任何一个环节，优化结果就会给出看似省钱、实则无法落地的方案。我建议刚开始接触这类模型的朋友，先不要急着堆设备数量，把一套最简单的“电-气-碳”闭环跑通，再逐步加约束换场景。跑通之后，你自然就明白下一步该在哪儿加随机性、在哪儿加设备模型了。