1. 从卖环保概念到算清碳账单为什么虚拟电厂必须引入阶梯碳交易这几年做虚拟电厂VPP优化调度的同行应该都有同一个感觉调度模型里最不好处理的已经不是功率平衡而是碳排放那本账。以前大家习惯的做法是给单位发电量乘一个固定碳价算个碳成本加进目标函数就完事了。但实际参与碳交易试点后会发现现实里的碳价根本不是固定的——排放越界越多每吨配额的价格会跳到更高的档位这就是阶梯碳交易tiered carbon trading。这个机制对调度决策的影响非常大。举个直观例子固定碳价下燃气轮机多烧一吨气碳成本是线性增长的优化器会一直算到一个边际平衡点但换成阶梯碳价后一旦总排放量跨过某个档位阈值后续所有排放量都按更高的碳价计费目标函数是分段线性且非凸的。这时候模型就必须配上0-1变量描述当前落在哪个碳价区间问题从LP直接变成MILP复杂度完全不同。与此同时单靠燃气轮机自身减排已经很难同时满足经济性和碳排放目标所以现在主流的低碳改造方案基本是两条线并行走一条是P2G-CCS耦合一条是燃气掺氢。P2G电转气利用低谷电或弃风弃光电解水制氢CCS把燃气轮机尾气里的二氧化碳捕集下来捕集到的CO2和电解产生的氢气送去甲烷化生成合成天然气再回用燃气掺氢则是把一部分电解氢直接混进燃气轮机的燃料里烧。这样电-氢-气-碳四者就在同一个虚拟电厂拓扑里形成了一个闭环调度优化的不是单个设备而是整个碳-氢-气循环的经济性。这篇内容面向正在做综合能源、新能源电力系统优化方向的工程师和研究生尤其在用Matlab做调度仿真的朋友。我会把模型怎么建、约束怎么列、阶梯碳价怎么线性化、Yalmip代码怎么写、算例结果怎么看一次讲清楚。全程基于实际可复现的24小时调度案例看完可以直接照着搭自己的算例。2. 案例模型怎么搭设备拓扑、能量流与调度时间尺度的选择2.1 标准VPP拓扑与24小时调度窗口我这次搭建的虚拟电厂聚合了六类单元风电机组、光伏阵列、燃气轮机支持掺氢、电化学储能、P2G-CCS耦合装置电解槽甲烷化反应器氢气储罐合成天然气储罐、以及与上级电网的交互接口。调度周期取典型日24小时时间步长1小时这是目前论文和工程验证里最常见的设置既能体现日内峰谷差异又不会让变量规模失控。如果做日内滚动优化可以把步长缩到15分钟但核心约束结构不变。能量流方向是这样的风、光优先供给负荷储能在低谷充电、高峰放电燃气轮机在峰时出力燃料来自天然气网购气和甲烷化合成的天然气同时按掺氢比例混入氢气P2G装置在电价低谷或风光过剩时启动电解制氢后兵分两路——一路直接进燃气轮机掺烧另一路和CCS捕集到的CO2进入甲烷化反应器生成合成天然气余电从上级电网购入补足。这里有一个建模时容易忽略的点P2G和CCS不是两个独立设备它们通过甲烷化反应器形成强耦合。换句话说CCS捕集了多少碳甲烷化就需要消耗对应比例的氢氢的来源又是P2G电解水产物。所以调度模型里必须同时描述氢平衡、碳平衡和气平衡三者的时间耦合关系直接影响了低谷时段的购电策略。2.2 风光出力、电价、负荷的典型日曲线处理我用的典型日数据是一组归一化后的模拟值风电在凌晨1点到5点出力较高中午光伏达到峰值负荷则有早晚两个峰。分时电价按峰谷平三段设置峰时1.1元/kWh、平时0.7元/kWh、谷时0.35元/kWh。在Matlab里我习惯把所有外部输入统一成行向量或列向量方便后面索引。比如N_t 24; P_wind [0.52 0.55 0.58 0.60 0.55 0.45 0.38 0.30 0.28 0.25 ... 0.22 0.20 0.18 0.17 0.16 0.18 0.22 0.28 0.35 0.42 ... 0.48 0.50 0.53 0.54] * 120; % 风电出力单位 MW P_pv zeros(1, N_t); P_pv(6:19) [0.05 0.15 0.32 0.48 0.62 0.75 0.85 0.88 ... 0.86 0.78 0.65 0.50 0.32 0.15] * 80; % 光伏出力 P_load [0.55 0.52 0.50 0.48 0.46 0.50 0.58 0.68 0.80 0.85 ... 0.92 0.88 0.78 0.70 0.75 0.82 0.90 0.95 0.98 0.88 ... 0.76 0.68 0.62 0.58] * 150; price [0.35 0.35 0.35 0.35 0.35 0.35 0.50 0.70 1.10 1.10 ... 0.70 0.50 0.50 0.70 0.70 1.10 1.10 1.10 1.10 0.70 ... 0.50 0.50 0.35 0.35];这里乘以的120、80、150是装机容量单位MW。实际用的时候替换成你自己的预测数据和装机参数即可。2.3 掺氢燃气轮机的出力-排放-热值联动关系掺氢燃气轮机的建模比普通燃气轮机多了一层燃料热值的修正。纯甲烷的低位热值约50 MJ/kg氢气的低位热值约120 MJ/kg但体积密度差别很大。为了方便调度模型计算我建议直接按能量占比定义掺氢比例记为k_h2(t)表示t时段燃气轮机输入燃料中以氢气形式提供的能量占比。单位时间内燃气轮机输入的总燃料能量可以写成F_total(t) P_gt(t) / (eta_gt * LHV_total(t))其中LHV_total(t)是掺氢混合燃料的平均低位热值随掺氢比例变化LHV_total(t) k_h2(t) * LHV_h2 (1 - k_h2(t)) * LHV_gas这样写的好处是天然气的购买量只计算非氢部分氢气的消耗量单独从储氢罐里扣除。排放量也按天然气侧的消耗来算掺氢比例越高同出力下的碳排放就越低。实际工程里燃气轮机掺氢比例有上限目前主流机组改造后能到20%至30%的体积掺氢比折算能量占比大约在7%到10%。我算例里取能量占比上限0.15比较稳妥。3. 核心数学模型目标函数、约束方程组与非线性项的线性化处理3.1 目标函数里那笔阶梯碳成本怎么写模型的目标函数是最小化虚拟电厂日运行总成本分六块向上级电网购电成本、购买天然气成本、各设备运维成本、阶梯碳交易成本、弃风弃光惩罚。购电成本好理解就是每个时段的购电功率乘分时电价再累加C_grid sum(price_buy(t) * P_buy(t))购气成本是天然气购气量乘气价注意这里只算天然气不算掺入的氢气因为氢气来自P2G自产已经在P2G运维成本里计过了C_gas price_gas * sum(F_gas_purchase(t))运维成本按各类设备的单位出力和单位充放电功率折算C_om sum(k_om_gt * P_gt(t) k_om_p2g * P_p2g(t) k_om_bess * (P_bch(t) P_bdis(t)))碳交易成本不是线性项它取决于全天累计净排放量。先算总排放量E_total再减去免费碳配额E_quota差值就是需要购买的配额量Q_carbon。阶梯碳交易规则下碳价按Q_carbon落入的档位取值于是C_carbon写成C_carbon lambda_1 * min(Q_carbon, q1) lambda_2 * max(0, min(Q_carbon - q1, q2 - q1)) lambda_3 * max(0, Q_carbon - q2)如果Q_carbon为负数说明配额有盈余可以按回购价出售获得收益。这个分段函数就是我们后面要做0-1线性化的目标。3.2 功率平衡、储能、燃气轮机基础约束功率平衡是每一类调度模型的骨架这里所有电源和储能放电之和等于负荷加上所有用电设备之和P_wind(t) P_pv(t) P_gt(t) P_bdis(t) P_buy(t) P_load(t) P_bch(t) P_p2g(t) P_ccs(t)储能约束包括SOC递推、充放电功率上下限、同一时段不能同时充放三个条件SOC(t1) SOC(t) (eta_ch * P_bch(t) - P_bdis(t) / eta_dis) * dt P_bch(t) P_bch_max * u_bch(t) P_bdis(t) P_bdis_max * u_bdis(t) u_bch(t) u_bdis(t) 1初末SOC我一般设为相同的0.2保证调度结果具有日循环特性。燃气轮机约束除了出力上下限还有爬坡约束和最小启停时间约束。最小启停时间在这个尺度下可以简化但爬坡必须加P_gt(t) - P_gt(t-1) ramp_up P_gt(t-1) - P_gt(t) ramp_down3.3 P2G-CCS耦合约束氢、碳、气的物料守恒这是整个模型最容易被写崩的地方。P2G装置消耗电功率P_p2g(t)按效率eta_p2g转化为氢气的化学能F_h2_prod(t) eta_p2g * P_p2g(t)这里的F_h2_prod用能量单位MW表示后续掺氢和甲烷化的氢气消耗也都统一用能量单位避免摩尔质量换算的麻烦。CCS捕集量正比于燃气轮机实际碳排放量E_ccs(t) beta_ccs * E_gt(t)E_gt(t)根据燃气轮机消耗的天然气量计算已经扣除了掺氢带来的减排效果。甲烷化反应器同时消耗氢气和二氧化碳比例固定。按能量折算我这里用一个简化系数k_r表示单位CO2捕集量对应消耗的氢能F_h2_reac(t) k_r * E_ccs(t) F_sng_prod(t) eta_m * E_ccs(t)F_sng_prod(t)是合成天然气产量同样按能量计。储氢罐和合成天然气储罐的递推约束分别为H2_st(t1) H2_st(t) F_h2_prod(t) - F_h2_gt(t) - F_h2_reac(t) SNG_st(t1) SNG_st(t) F_sng_prod(t) - F_gas_gt(t)其中F_h2_gt(t)是掺氢消耗的氢气能量F_gas_gt(t)是燃气轮机消耗的合成天然气能量。注意到这里天然气网购气和合成天然气都可以供燃气轮机使用天然气网购气直接进入F_gas_gt的补足通道。这个方程组把电、氢、气、碳四条能量流串成了一张网任何一个设备的调度决策都会沿着物料守恒链传导到其他设备。3.4 阶梯碳价的0-1线性化方法附Yalmip代码分段线性函数写进MILP的标准做法是引入二进制变量。我习惯用Yalmip的define和binary变量来手工实现因为这样最可控。基本思路是将Q_carbon可能在的区间编号假设分成三段则引入两个0-1变量z1、z2表示档位再用大M约束把成本表达式限制在对应分段里。% 阶梯碳价参数 lambda [80, 120, 200]; % 各档碳价元/吨 q_th [0, 5000, 10000]; % 档位边界吨 M 100000; % 大M z binvar(2, 1); % 两个0-1变量表示落在哪个分段 % 确保只落在其中一个分段 Q_str sdpvar(1, 1); % 实际净排放含正负 % 分段区间选择约束 constraints [constraints, ... Q_str q_th(1) - M*(1-z(1)), ... Q_str q_th(2) M*z(1), ... Q_str q_th(2) - M*(1-z(2)), ... Q_str q_th(3) M*(1-z(2))];更稳妥的做法是引入分段变量seg1、seg2、seg3把Q_str拆成三段非负部分每段乘对应碳价后求和。这种方式不会出现大M选段时的数值问题推荐优先使用seg sdpvar(3, 1); Q_str sum(seg); constraints [constraints, 0 seg(1) 5000]; constraints [constraints, 0 seg(2) 5000]; constraints [constraints, 0 seg(3) 20000]; C_carbon lambda(1)*seg(1) lambda(2)*seg(2) lambda(3)*seg(3);注意这里如果Q_str是负数配额盈余分段变量法就不能直接用需要额外判断符号。我算例中免费配额取得比较保守全程Q_str都为正所以分段变量法够用如果你要处理配额盈余场景就得加一组符号变量。4. Matlab实现与求解选型从环境配置到跑通24小时调度4.1 环境Matlab Yalmip Gurobi/CplexMatlab版本建议R2021b以上Yalmip用新版求解器优先Gurobi或Cplex两者都可以申请学术授权学生和高校科研用户直接走官方申请流程就行。如果你的机器装的是开源求解器也可以让Yalmip自动调用SCIP但MILP求解速度会慢不少24小时模型问题不大一旦把时间粒度缩到15分钟差距就会非常明显。我第一次跑这个模型时用的是Cplex 12.10算24小时、每小时一个点变量规模大概是连续变量约150个、二进制变量约80个求解时间在十几秒到一分钟之间。同一个模型如果改用PSO这类智能算法不仅每次求解时间翻几倍而且很容易陷入局部最优后面我会专门聊这个问题。4.2 代码骨架变量、约束、求解器调用顺序下面是整个主脚本的结构我按这个顺序写代码基本不会乱%% 初始化 clc; clear; close all; load(typical_day_data.mat); % 你的典型日数据文件 %% 决策变量定义 P_gt sdpvar(N_t, 1); P_p2g sdpvar(N_t, 1); P_buy sdpvar(N_t, 1); P_bch sdpvar(N_t, 1); P_bdis sdpvar(N_t, 1); H2_st sdpvar(N_t1, 1); SNG_st sdpvar(N_t1, 1); F_h2_gt sdpvar(N_t, 1); F_h2_reac sdpvar(N_t, 1); F_sng_prod sdpvar(N_t, 1); E_ccs sdpvar(N_t, 1); u_gt binvar(N_t, 1); % 燃气轮机关停 u_bch binvar(N_t, 1); % 储能充电状态 u_bdis binvar(N_t, 1); % 储能放电状态 %% 约束 constraints []; % 功率平衡、储能、P2G-CCS等约束按前一节写进去 %% 目标函数 objective C_grid C_gas C_om C_carbon C_wind_penalty; %% 求解 ops sdpsettings(solver, gurobi, verbose, 2, showprogress, 1); optimize(constraints, objective, ops);变量初始化这里有个容易犯的错储气罐状态变量长度我写的是N_t1因为末时刻的递推会用到N_t1的值。你可以在约束里把H2_st(1)设为初值并强制H2_st(N_t1)等于初值形成日循环约束。如果不加这个末态约束求解器可能会把储气罐在最后几个时段内存量全放空结果看起来经济性很好实际却无法连续运行第二天。4.3 典型坑二进制变量暴增、求解超时、不可行解的排查这个模型最容易出现的问题是二进制变量失控。如果给燃气轮机加启停状态、给储能加充放状态、给碳阶梯加分段选择再给购电加状态变量二进制变量数量立刻翻倍。Yalmip里看起来只是多了几个变量求解器内部的分支定界树会膨胀得很快。我的建议是第一版模型先不加机组启停状态把燃气轮机设成最小出力可降到10%额定出力当作可调节连续机组处理。储能和阶梯碳价的状态变量是必须的其他能简化的先简化。等第一版跑通、结果合理再逐步加启停约束。不可行解是另一个高发问题。最常见的原因是功率平衡约束里外购电力上限设太小低谷时段负荷低、P2G又启动两边抢电导致平衡被破坏。排查方法很简单把P2G的最小出力设成0把购电上限调大先看约束是否可行如果还是无解就把功率平衡约束的等式改成不等式加松弛变量输出松弛变量的值定位是哪个时段出了问题。我在调试时会在每个约束块后面加一行注释记录设备索引范围方便用check(constraints)定位冲突源。Yalmip的check会返回每个约束的残差极大地方便了找错。5. 仿真结果与对比实验阶梯碳价和掺氢到底改写了哪些决策5.1 基准算例参数表下面是我这次算例用的主要参数你可以直接抄过去替换自己的数据参数数值说明燃气轮机额定功率60 MW最低出力10 MW储能容量120 MWh最大充/放功率30 MWP2G额定功率40 MW制氢效率0.7CCS捕集率0.85捕集燃气轮机碳排放掺氢能量比上限0.15燃气轮机燃料能量占比免费碳配额基准0.3 t/MWh按燃气轮机出力发放阶梯碳价80/120/200元/吨5000/10000分界弃风弃光惩罚300 元/MWh保证最大消纳天然气价格2.8 元/m3按热值折算后约0.9 元/kWh谷/平/峰电价0.35/0.70/1.10元/kWh免费配额确定方式是重点我按燃气轮机每发1 MWh电力发放0.3吨碳配额。这样设置会让燃气轮机在高负荷时段自动产生较大的配额缺口阶梯碳价就真正成了影响出力决策的硬约束。5.2 低碳场景下的各设备出力曲线解读跑完基准算例后我先看燃气轮机和P2G的出力曲线两者呈现明显的互补关系夜间0点到6点风电出力高、电价低燃气轮机接近最小出力P2G满载制氢上午8点以后负荷攀升燃气轮机开始逐步加出力P2G逐渐降负荷晚上18点到21点负荷达到峰值燃气轮机满发储能放电补足缺口。这个结果符合预期低谷电和弃风电驱动P2G运行产出的氢一部分储存起来供白天掺氢使用一部分进甲烷化反应器消化CCS捕集的CO2高峰时段燃气轮机带着掺氢燃料高发碳排放被CCS部分捕集形成一个完整的日内碳循环。储氢罐的SOC曲线也很有意思早上6点左右达到最高点之后随着掺氢消耗逐步下降到晚上21点基本见底次日凌晨又重新补充。这种夜间制氢-白天掺烧的日内节奏正是P2G-CCS耦合系统参与调度的典型特征。5.3 阶梯碳价与线性碳价调度结果对比为了说明阶梯碳价的价值我把碳交易成本从阶梯函数改成固定碳价100元/吨做了对照实验。固定碳价场景下燃气轮机全天总出力比阶梯碳价场景高约8%碳排放量高出约12%。原因在于固定碳价下每个时段的边际排放成本恒定优化器只会做全局平均式的取舍而阶梯碳价下一旦全天累计排放量逼近下一档阈值后续时段的实际边际碳成本就跳到120元甚至200元燃气轮机在峰时段的相对经济性迅速恶化优化器宁可多用一部分外购电和储能放电。这正是阶梯碳交易机制进入调度模型的真正意义它让碳成本从账本里的平均成本变成了每个时段动态变化的边际成本直接影响的是出力优先级排序而不只是总成本数值。5.4 掺氢比例上限对购电策略与碳交易成本的影响掺氢比例上限从0逐步调到0.25观察三个指标天然气购气量、外购电量、碳交易成本。掺氢比例提高后天然气购气量明显下降因为部分燃料能量由自产氢气承担。碳交易成本同步下降因为单位出力的碳排放系数降低了。比较意外的是外购电量不降反升原因是P2G为了制出更多氢气增加了低谷期用电负荷这部分电量一部分来自风电一部分来自夜间低价外购电。从结果看掺氢比从0提到0.15时系统总运行成本下降了约4%碳排放下降约9%但从0.15提到0.25时成本下降幅度明显收窄只有1%左右碳排放继续下降但幅度也在变小。这说明掺氢存在一个经济性和减排效果的边际拐点并不是掺得越多越好。6. 代码复用经验与下一步扩展建议6.1 适合直接改的参数向量和可替换模块整套代码里最核心的可复用部分就是约束构建段。你拿到别人的代码后优先看三处一是目标函数里各项成本前的系数二是P2G-CCS耦合约束里的k_r和eta_m这两个联动系数三是阶梯碳价的分段阈值和碳价向量。h2_ratio这一行是扩展掺氢相关约束的入口你可以直接把k_h2(t)从常数上限改成随时间变化的决策变量实现掺氢比的实时优化。储气罐容量参数也很关键。很多算例默认储气罐无限大P2G夜间制氢全部存起来这样看起来低碳指标很好但实际工程里储氢罐投资巨大。你可以把H2_st上限从500 MWh缩到200 MWh再跑一次观察碳排放和购气量的变化这个结果可以直接用来论证储氢容量配置方案。6.2 从08:00-20:00峰时挖数据的小技巧分析结果时别只盯着全天总量我把每个时段的购电功率、燃气轮机出力、P2G耗电量都导出成表格然后用峰时8点到20点和谷时20点到次日8点分别做平均能看出更细的调度规律。比如峰时段储能的放电深度、掺氢比例的利用率、CCS捕集量与燃气轮机出力的比例关系这些指标是论文章节里最有说服力的数据。我建议在代码末尾加一个结果汇总矩阵% 结果汇总矩阵按行时段/总购电/燃机出力/P2G耗电/碳捕集量/碳交易成本 result_table table((1:24), value(P_buy), value(P_gt), ... value(P_p2g), value(E_ccs), value(C_carbon)); writetable(result_table, dispatch_result.csv);有了这个CSV文件后续画图、敏感性分析、横向对比都不用来回重跑模型。6.3 换成智能算法求解需要注意什么如果你的场景必须用智能算法比如要处理燃气轮机启停次数等整数变量特别多的组合优化问题我的建议是不要直接拿PSO或GA跑原始MILP模型。把二进制变量尽量控制在10个以内剩余连续性约束全部保留然后用罚函数处理可行性。但这个做法有代价处理等式约束的时候罚函数权重很难调权重太小约束不满足权重太大目标函数被淹没。我个人的经验是只要模型规模没有到几千个变量级别优先用Gurobi/Cplex配合Yalmip把MILP交给专业求解器处理。只有在需要对比智能算法性能、或者论文要求展示元启发式方法时才考虑把阶梯碳价部分用惩罚函数近似。在实际操作中我最常踩的坑是终于跑出结果后发现储能初末SOC没闭合整体曲线看起来没问题但第二天无法连续运行。现在我的习惯是把储气罐初末状态约束单独列一个段落每次调参后先跑一遍可行性检查再去看经济性和碳排放指标。做这类多能源耦合调度模型能不能连续运行、碳平衡能不能闭合比成本数字本身更能说明模型质量。