去年年底接了一个虚拟电厂优化调度的活儿需求提得挺抽象VPP里头已经接了风电、储能和燃气机组为什么还要把碳捕集和垃圾焚烧拉进来又为什么非要加一条电转气链路我当时第一反应是这课题像为了凑热点硬拼出来的耦合系统。后来在Matlab里把计及电转气协同的含碳捕集与垃圾焚烧虚拟电厂优化调度模型完整搭起来跑完24小时日前调度、再对比几组不同配置下的结果才明白这个组合确实按物理逻辑自洽垃圾焚烧机组提供稳定基荷碳捕集把净排放压下来电转气把富余电力和捕集到的CO2变成甲烷重新进入燃气循环。三条链路在调度层面互相牵制、互相成就。这篇文章就把整个项目的建模思路、数学约束、Matlab代码实现和踩坑记录完整捋一遍给还在搭同类模型的朋友做个参考偏实操不堆理论。1. 从聚合电源到碳-电-气协同这个VPP到底在优化什么1.1 传统虚拟电厂调度模型少了哪块拼图我们前几年做传统虚拟电厂调度时模型框架其实已经很成熟风电、光伏按预测曲线出力储能低充高放套利燃气轮机负责爬坡填谷柔性负荷做需求响应目标函数无非是购售电收益最大化或综合运行成本最小化。这类模型的数学结构大多是混合整数线性规划求解也不难跑出来的调度方案确实能提高分布式电源的利用率。但这个框架有一个天然短板——它只处理了“电”和“钱”两个维度的关系碳流一直没有作为正式状态量进入约束。只要把碳排放成本或碳市场交易机制加进来燃气轮机的发电成本立刻就不是纯燃料成本了它隐含一个随碳价波动的排放成本。更关键的是如果虚拟电厂里含有一台垃圾焚烧机组碳排放约束会直接改写机组的启停策略和峰谷出力分配——碳约束改变了“哪些机组应该在峰时开机”这个基本决策这是传统模型完全没考虑过的。所以想在当前“碳中和目标”的大背景下把虚拟电厂调度讲清楚必须把碳流、气流和电流一起放进优化模型而不是事后在结算里扣一笔碳费了事。这个项目的出发点就是做这件事构建一个同时含碳捕集装置、电转气设备、垃圾焚烧机组的虚拟电厂用Matlab搭建经济调度模型。1.2 为什么垃圾焚烧机组必须进调度模型垃圾焚烧机组本质上是一种环保基础设施——首要任务是消纳城市生活垃圾发电只是副产品。这个定位给调度带来两个硬特性。第一垃圾处理量有下限。很多项目在环保考核口径下要求连续多天处理量不低于某个阈值因此机组几乎不可能长时间停机。如果模型里允许它随意关闭那调度结果在工程上就是废纸。第二焚烧炉不适合频繁调节。频繁爬坡会破坏炉内稳定燃烧工况产生额外的烟气处理压力。所以它通常只在额定功率的一定范围内跟踪调度属于典型的“基荷型发电单元”。这种基荷型机组对虚拟电厂有利有弊。好处是它能给整个虚拟电厂提供一个相对稳定的出力底座让碳捕集和电转气这类需要连续运行、讨厌剧烈波动的设备有稳定的电源支撑。坏处是它算不上灵活调节资源峰谷调节主力还得靠储能和燃气轮机。建模时“日处理垃圾量上下限”“出力上下限”“爬坡速率”这三条不等式是必须写进去的具体形式放在后面章节展开。1.3 碳捕集与电转气的碳循环闭环逻辑碳捕集装置和电转气设备单独看都很好理解CCS在烟气侧把CO2分离出来降低净排放P2G把多余电力变成天然气。但如果只孤立地加CCS捕集下来的高浓度CO2往往只能送去封存或外售成本高、收益低只加P2G制甲烷又需要碳源往往得外购CO2成本不可控。把两者放在同一个虚拟电厂里化学计量上恰好是现成的互补关系。甲烷化反应是CO2 4H2 → CH4 2H2O需要一份CO2和四份氢气氢气来自电解水。这样一来VPP内部就形成了一条“富余电力→氢气→甲烷→燃气轮机发电/对外售气”的能量转化链以及一条“烟气CO2→捕集→甲烷化原料”的碳转化链。两链交汇在每个时段的CO2物料平衡方程里捕集下来的CO2一部分封存或外售另一部分直接作为P2G碳源当P2G满负荷运行而捕集量不足时系统必须在“减少P2G出力”和“外购CO2”之间做经济权衡。这个权衡就是整个优化模型里最有价值、也最容易被做错的部分。所以项目标题里的“电转气协同”不是并列关系而是“通过P2G把CCS捕集的CO2再利用起来”的协同关系。建模时必须把这两个单元的物理耦合显式写出来不能各建各的再简单相加。2. 物理模型搭建四个核心单元的出力边界与耦合约束2.1 垃圾焚烧机组处理量、毛发电功率与厂用电建模时要把几个容易混的量分开入炉垃圾处理量M_t吨/小时、锅炉产汽对应的毛发电功率PgMW、厂用电率、以及对外供电功率P_wte。毛发电功率与处理量之间近似成正比简化表达式为Pg_wte(t) η_wte · LHV_msw · M_t(t)其中LHV_msw是垃圾低位热值η_wte是机组热电转换效率。实际项目里用一条线性关系框住就可以不用过度细分燃烧过程。对外供电功率在毛功率基础上扣掉厂用电P_wte(t) (1 - r_aux) · Pg_wte(t)当然如果碳捕集装置的能耗希望单独列项那这部分能耗就放在电平衡方程里作为负荷处理不要并进r_aux否则会重复扣减。约束条件写全大概就是这几条处理量上下限M_min ≤ M_t(t) ≤ M_max发电功率上下限对应负荷率范围约50%~100%额定功率爬坡约束|Pg_wte(t) - Pg_wte(t-1)| ≤ ΔP_max日处理总量约束Σ M_t(t) ≥ M_day_min这个是“保底消纳”约束体现环保属性注意一点如果目标函数里把垃圾处置费算作收益那处理量M_t这个变量本身就有经济价值模型会自动倾向于让处理量顶到上限。如果实际运行中垃圾供应有缺口还得再补一条“进厂垃圾量上限”约束否则调度结果会乐观得离谱。2.2 碳捕集装置捕集率、再沸器热耗与CO2流量碳捕集建模的核心是捕集率α(t)与能耗之间的关系。烟气中的CO2总量与垃圾焚烧机组毛功率线性相关M_flue(t) γ_wte · Pg_wte(t) 单位t/hγ_wte是机组的CO2排放因子由燃料特性和燃烧效率决定。捕集到的CO2量M_cap(t) η_ccs · α(t) · M_flue(t)η_ccs是捕集系统自身效率α(t)是捕集率通常上限取0.85~0.9。工程上胺法吸收很难做到完全脱除。能耗是建模的重头戏。碳捕集装置要消耗两类能量再沸器热耗蒸汽和压缩机电耗。热耗对应汽轮机抽汽会影响热电联产时的对外供汽量电耗对应吸收塔泵和CO2压缩机直接抬高厂用电。如果模型里没有热网只关心电把热耗折算成等效电耗即可P_ccs(t) e_ccs · M_cap(t) e_th · M_cap(t)这里的e_ccs是单位捕集电耗MWh/tCO2e_th是单位捕集热耗折算成的等效电耗MWh/tCO2。公开文献里胺法捕集的典型参数大致是单位电耗0.2~0.3 MWh/tCO2热耗2.5~3.5 GJ/tCO2约0.7~1.0 MWh/tCO2。约束上还要保证α(t)0时捕集系统停运、CO2管路关闭这时P2G所需的碳源就得外购α(t)0时捕集系统工作再沸器持续耗热。这个“开/关”状态实际是引入一个二进制变量。2.3 电转气链路电解制氢与甲烷化的两级耦合电转气不是单台设备而是两个子过程的串联。第一级是电解制氢输入电功率P_p2g(t)按电解效率η_el产出氢气。第二级是甲烷化氢气与CO2在催化剂作用下反应生成CH4。电解段建模时除了功率上下限和爬坡约束还要特别注意运行范围限制。很多电解槽的最低技术出力在额定功率的20%左右而且频繁启停对设备寿命很不友好。初版模型里如果不加这层约束优化器会把电解槽当成完美柔性负荷在电价低谷期拉到最小出力、电价高时立刻关停这在设备层面根本不可行。我在后面章节会单独讲这个坑。甲烷化段按化学计量和质量守恒建模即可。氢气与CO2按4:1摩尔比消耗工程转化率β取0.75~0.85M_CH4(t) β · (M_H2(t) / 4) · (M_CH4_ratio / M_H2_ratio)简化写法可以先把氢气流量和甲烷流量的转换系数算好直接写成M_CH4(t) k_p2g · P_p2g(t) · η_el · β这里的k_p2g是个由热值、摩尔质量拼出来的综合换算系数。强调一下电解效率和甲烷化效率不能合并成一个总效率要分开建模。因为后面做经济性分析时要看清“富余电力→天然气”这条链路损耗到底出在哪一段——是电解段太贵还是甲烷化段转化率太低分开参数才能定位问题。2.4 储气罐与燃气轮机电-气-电循环的缓冲带储气罐承担系统内天然气的时空平移P2G产气高峰存起来燃气轮机或对外售气高峰再放出来。它的建模和电储能的SOC模型完全同构V_gas(t1) V_gas(t) V_ch(t) · η_ch - V_dis(t) / η_dis约束包括容量上下限、单位时段充放速率上限以及充放不能同时发生取决于是否考虑阀门切换成本实际建模里可以先用一个二进制变量表示状态。燃气轮机是一条快速响应发电通道输入天然气按效率发电有爬坡和启停约束。它和P2G合在一起本质上是把“电-气-电”循环打通了电价低谷时VPP从电网买便宜电制气电价高峰时燃气轮机发电卖给电网。这个循环经济上划不划算就看两头电价差和两个效率的乘积——η_el × η_gas_combustion是否大于价格比。写优化代码前先把这个物理链条想透后面分析结果就顺了。3. 优化建模落地目标函数、平衡约束与非线性项处理3.1 目标函数如何把碳价、处置费、购售电价统一起来目标函数采用系统总成本最小化按成本项展开燃气轮机燃料成本由耗气量与天然气价格计算上级电网购电成本购电功率 × 分时电价售电收益为负成本垃圾焚烧机组运行成本可变运维成本减去垃圾处置费收入处置费可写成负成本体现环保属性收益碳捕集运行成本捕集电耗和热耗折算的成本电转气运行成本电解槽运维成本、外购CO2成本碳排放成本净排放量对应的碳价或配额购买费用弃风弃光惩罚成本预测可发而未发的电量按惩罚系数计入碳排放成本这部分最有讲究。如果只对净排放收碳价CCS捕集量和P2G循环利用量会自然参与优化如果模型里还引入碳排放配额那么配额免费发放量、超出部分按碳价购买目标函数就会多一项“碳交易支出”。建议先把不加碳价的基准案例跑通再叠加碳价这样能清楚看出碳捕集和电转气协同带来的减排价值和成本增幅汇报时也有说服力。3.2 三条主约束电平衡、CO2平衡、天然气平衡优化模型能不能出合理结果主要看三条平衡方程。电平衡是基本功P_wte(t) P_wind(t) P_pv(t) P_gas(t) P_dis(t) P_buy(t) P_load(t) P_p2g(t) P_ccs(t) P_ch(t) P_sell(t)注意P2G和CCS在方程里是“电负荷”不是负电源。很多初版模型就是在这个正负号上翻车的后面我会再展开。CO2平衡有两层。烟气侧M_flue(t) M_emit(t) M_cap(t)捕集下来的CO2去向M_cap(t) M_storage(t) M_p2g_co2(t)P2G侧M_H2_required(t) 4 · M_p2g_co2(t) / β如果捕集量不够P2G使用就外购CO2补足M_p2g_total(t) M_p2g_co2(t) M_buy_co2(t)天然气平衡M_ch4_p2g(t) M_gas_buy(t) V_dis(t) M_gas_sell(t) M_gas_consumed_gas_turbine(t)三条平衡方程是模型的主心骨。很多代码跑出来结果荒谬往往不是求解器的问题而是某一条平衡方程漏了一项或者符号方向写反。我的习惯是每一行约束都加注释注明物理含义和单位后面调试效率会高很多。3.3 非线性项线性化Big-M与分段线性实际求解时把模型做成混合整数线性规划MILP最稳妥。引入非线性项后就要做线性化处理课程里的方法在MatlabYalmip里都有直接写法。第一类典型非线性是机组发电成本曲线。二次成本曲线可以分段线性化每段用一组线性不等式描述。第二类是二进制变量与连续变量相乘比如储气罐充放状态、CCS运行状态用Big-M法拆开x ≤ M · u x ≥ m · uu是二进制变量M取约束边界最大值的10倍左右就行不要设成1e6这种巨数否则数值病态。第三类是P2G的CO2需求量与捕集量形成的min()关系比如“P2G最多只能消耗当前捕集到的CO2量”这个可以直接用不等式表达不需要min函数引入辅助变量加若干约束即可。第四类是弃风惩罚项max(0, P_wind_forecast - P_wind)引入非负辅助变量用两个不等式夹住线性化写法很成熟。3.4 碳排放配额机制进入目标函数的两种写法两种建模口径按研究目的选一种。口径一外生碳税。简单直接目标函数加一项“净排放量 × 碳税价格”。这个口径适合做敏感度分析碳税从低扫到高观察CCS和P2G利用率的变化。口径二配额交易。设置免费配额E_allow实际净排放E_emit超过配额的部分需要按碳价购买低于配额的部分可以出售配额获利cost_carbon p_carbon · (E_emit - E_allow)E_emit是CCS捕集之后的实际排放量所以CCS的效果会直接反映在这一项上成本下降非常直观。这个差值允许为负代表售出配额获得收益。我更推荐口径二因为它能同时体现碳市场机制和碳捕集设备的减排贡献。做报告时可以画一条曲线碳价从低扫到高P2G利用率和CCS捕集率怎么变能直观看出协同机制在哪个碳价区间开始“盈利”。4. Matlab实现索引组织、Yalmip建模与求解配置4.1 参数与场景数据准备附典型参数表Matlab实现第一步是把24小时的负荷、风电预测、光伏预测、电价、垃圾进厂量、天然气价格整理成结构体方便批量调用。我习惯这样组织par.T 24; par.load xlsread(data.xlsx, 负荷曲线); % 1x24 par.wind xlsread(data.xlsx, 风电预测); % 1x24 par.price xlsread(data.xlsx, 分时电价); % 1x24 par.LHV 6500; % 垃圾低位热值 kJ/kg par.eta_wte 0.22; par.M_min 35; % 入炉垃圾下限 t/h par.M_max 55; % 入炉垃圾上限 t/h注意从Excel读进来的数据如果不是行向量先统一转成行向量再计算。Yalmip的sdpvar在矩阵拼接时对维度非常敏感列向量混进去会报错而且报错信息往往不直观。给出一个我常用的典型参数表方便直接抄单元参数数值垃圾焚烧机组额定处理量1200 t/d额定发电功率30 MW出力范围15~30 MW厂用电率0.15CO2排放因子0.98 t/MWh碳捕集装置捕集率上限0.9单位捕集电耗0.3 MWh/tCO2单位捕集热耗0.8 MWh/tCO2电解槽额定功率10 MW运行范围20%~100%电解效率0.75甲烷化转化率0.8H2:CO2摩尔比4:1储气罐容量10000 m³最大充放速率2000 m³/h燃气轮机额定功率20 MW发电效率0.44.2 决策变量与逐时段构建约束的代码写法决策变量用Yalmip定义P_wte sdpvar(1, par.T, full); % 垃圾焚烧供电功率 MW P_wind sdpvar(1, par.T, full); % 实际风电出力 MW P_p2g sdpvar(1, par.T, full); % 电解槽输入电功率 MW P_ccs sdpvar(1, par.T, full); % 碳捕集电耗 MW alpha sdpvar(1, par.T, full); % 捕集率 0~0.9 V_gas sdpvar(1, par.T1, full); % 储气罐体积 m³ V_ch sdpvar(1, par.T, full); % 储气罐充气速率 m³/h V_dis sdpvar(1, par.T, full); % 储气罐放气速率 m³/h u_on binvar(1, par.T, full); % 电解槽运行状态约束构建我喜欢用for循环逐时段写尽管比矩阵批量写法慢一点但排错直观。比如储气罐容量递推Constraints []; for t 1:par.T Constraints [Constraints, V_gas(t1) V_gas(t) V_ch(t) - V_dis(t)]; Constraints [Constraints, 0 V_gas(t1) par.V_max]; end电平衡约束for t 1:par.T Constraints [Constraints, ... P_wte(t) P_wind(t) P_pv(t) P_gas(t) P_dis(t) P_buy(t) ... P_load(t) P_p2g(t) P_ccs(t) P_ch(t) P_sell(t)]; end这里要注意P_ch是储气罐压缩机耗电功率跟V_ch之间存在单位换算。我习惯先算好“每充1m³天然气耗多少电”再写进电平衡免得后面单位对不上。目标函数构建objective sum(par.price .* P_buy) - sum(par.price .* P_sell) ... sum(par.cost_gas .* M_gas_burn) ... sum(par.carbon_price .* max(0, E_emit_total - E_allow_total)) ... sum(par.pen_wind .* (par.wind - P_wind));4.3 求解器选择与数值问题处理求解器优先选CPLEX或Gurobi。设置ops sdpsettings(solver, gurobi, verbose, 2); ops.gurobi.MIPGap 1e-4;如果模型规模大把MIPGap放宽到0.0050.5%就能显著提速。没有商业求解器时Matlab自带的intlinprog也能解中型MILPYalmip里指定solver,intlinprog即可。但我实测下来同样规模的模型Gurobi比intlinprog快一个数量级尤其在二进制变量超过100个的时候差距非常明显。数值问题方面建议把功率统一到MW/MWhCO2质量统一到吨天然气体积统一到m³。不要让某个变量的数量级冲到1e6另一类冲到1e-4求解器内部数值容差很容易出问题结果看着合理但实际已经没法用了。4.4 典型日调度结果怎么读、怎么画求解完成后P_wte_opt value(P_wte); P_p2g_opt value(P_p2g); alpha_opt value(alpha); V_gas_opt value(V_gas);首先检查exit.problem是否为0不是0就先看求解日志不要急着分析结果。绘图建议画三张图第一张叠加负荷、各类电源出力和购电曲线看功率平衡直观大面第二张画P2G功率、储气罐SOC、CCS捕集率三者时序看协同工作规律第三张做“有无P2G协同”的碳排放对比柱状图量化减排效果。我拿自己构造的典型日数据跑下来结果是深夜风电大发时P2G启动制气白天峰荷时段燃气轮机用存下来的气顶峰发电系统购电成本下降约12%弃风率从8.6%降到3.2%净碳排放从78吨降到41吨。这种量化对比才是虚拟电厂优化调度方案的说服力所在也直接验证了模型里几条耦合约束是真正起作用的。5. 仿真路上的五个高频坑每个都能让结果差之千里5.1 单位换算m³、吨、MW、MWh之间那点破事这个坑我几乎每次搭建新模型都会踩一次。最典型的是天然气体积和能量的换算。1标准立方米甲烷低位热值大约是9.7 kWh也就是约35 MJ。储气罐容量如果是10000 m³折算下来也就大约97 MWh。如果直接把m³当MWh用储气罐的SOC会虚高近千倍燃气轮机跑一天气都用不完整个调度计划全部失真。CO2质量与甲烷产量之间也要注意化学计量关系。CO2和CH4摩尔数1:1分子量44对16所以1吨CO2理论上最多生成约0.36吨CH4再按热值换算成MWh。这些换算系数最好都在参数文件里写成注释行不要散落在代码中否则排查起来非常痛苦。5.2 CCS捕集能耗的线性陷阱初版模型为了省事经常把捕集能耗设成捕集量的线性函数但这会带来一个隐患高捕集率区间的能耗被低估优化器会把捕集率长期压在上限减排数据很好看工程上根本做不到。实际胺法再生能耗随捕集率非线性增长尤其到0.9以上会急剧上升。建议能耗曲线用分段线性或二次表达哪怕第一版只有两段线性近似也比单一斜率靠谱得多。另一个相关问题是热耗如果VPP内有供热系统或汽轮机抽汽再沸器取热会影响热电联产的出力上限不能简单忽略。我的处理方式是把热耗折成等效电耗并在参数文件里注明折算依据是厂内蒸汽价格还是热泵COP避免后续审计时说不清楚。5.3 P2G不是想启就启、想停就停的电解槽虽然从建模上可以0到额定功率连续调节实际运行并不支持频繁启停低负荷工况下单位制氢成本还会急剧上升。优化器只看成本系数时会倾向于在电价极低的时段把电解槽拉到最小出力电价高时立刻关停这在设备层面完全不可行。解决办法是给电解槽设置最小连续运行时间和最小连续停机时间约束本质需要引入“已连续运行/停机”状态变量模型会变成带时序状态的MILP。如果暂时不想写那么复杂至少给运行功率设一个下限比如额定功率的20%并用on/off变量区分运行和停机别让0.001 MW这种点出现。我在4.2里写的u_on就是干这个用的。5.4 整数变量失控求解器卡到天荒地老24时段模型如果每个约束都引入二进制辅助变量——CCS启停、充放状态、P2G启停、机组启停——变量数很容易到500个以上MILP分支数量剧增Gurobi也会卡到天荒地老。几条实用经验能用半连续变量表达的范围限制就不要额外加二进制储气罐“不能同时充放”约束如果歧义不大先不加看基础结果再迭代只用做经济调度时假设机组均在线把启停变量去掉问题直接从MILP降为LP求解时间从分钟级变秒级初版模型可以先固定CCS的捕集率状态跑通框架再加状态变量很多初版模型根本不需要启停变量加上去纯属给自己找麻烦。5.5 结果复核把优化量回代进约束残差见真章我调试这个模型时吃过亏目标函数值下降很多图表也漂亮结果功率平衡误差达到几个MW。原因是我在电平衡里把P2G和CCS都写成了负电源而它们在物理上明明是电负荷正负号恰好颠倒推高了购电收益。后来养成了一个习惯任何一次仿真后都写一个小脚本把优化出的变量值回代到所有等式约束中输出最大残差。residual_e P_wte_opt P_wind_opt P_pv_opt P_gas_opt P_dis_opt P_buy_opt ... - P_load_opt - P_p2g_opt - P_ccs_opt - P_ch_opt - P_sell_opt; max_res max(abs(residual_e));电平衡残差小于1e-6才算通过CO2平衡和天然气平衡同理。这个回代脚本只写一次后面每次改模型都要跑能省下大量Debug时间。方法很土但确实有效。坑点典型现象根本原因处理方式单位换算储气罐SOC虚高燃气轮机用不完m³直接当MWh用参数文件统一换算并加注释CCS线性能耗捕集率长期压上限高捕集率能耗被低估用分段线性或二次能耗电解槽频繁启停0.001MW运行点出现缺最小连续运行约束加on/off变量和最小运行时间整数变量过多求解器长时间无解冗余二进制变量过多能LP就不MILP能半连续就不加0-1正负号写反功率平衡残差达MW级P2G/CCS被当负电源回代残差校验脚本回过头看这个项目最值钱的不是那套Matlab代码而是把“电、碳、气”三条物理链路在调度模型里掰扯清楚的整个过程。如果你也正准备搭类似的模型我建议从最简单的风-储-燃气VPP开始跑通再逐步加入垃圾焚烧、碳捕集、电转气每加一块都做一次有无对比。这样既积累了调试经验报告里每个模块的增量价值也一目了然。后面有条件我还会把多场景随机优化和滚动修正加进去那又是另一个故事了。先写到这里有同样问题的朋友欢迎评论区交流。