尧图网站建设 
这两年做低碳电力系统相关的网格优化最常被问到的一个问题就是虚拟电厂调度里同时塞进电转气、碳捕集和垃圾焚烧到底该怎么建模单看任何一个设备网上的资料都不算少但要在一套调度模型里把“电、气、碳”三条线串起来很多朋友一开始会发懵。这篇文章我不打算写那种“从原理到展望”的综述而是直接以一个我自己用Matlab实现的虚拟电厂优化调度项目为底子把变量怎么定、约束怎么写、求解器怎么调、算例怎么安排以及我实际踩过的坑按顺序摊开来讲。适合正在做综合能源系统优化、低碳调度、碳捕集与封存建模的研究生以及想用Matlab把调度模型落到代码里的工程师。1. 项目背景为什么要把电转气、碳捕集和垃圾焚烧放在同一个算盘里1.1 虚拟电厂不再只是“风光储”的拼盘传统的虚拟电厂调度大概就是聚合分布式光伏、风电、储能和备用柴油机组以最小运行成本为目标去优化各设备的出力。这种模型已经很成熟核心约束也只有功率平衡、机组出力上下限、储能SOC等。可一旦加入垃圾焚烧机组问题立刻不一样。垃圾焚烧机组的特点是“不能想停就停”。焚烧炉需要维持稳定燃烧垃圾仓里的垃圾每天都要处理城市固废的连续清运要求决定了焚烧机组必须连续运行。在夜间的低负荷时段其他机组可以大幅降出力甚至停运垃圾焚烧机组往往只能压到最低稳燃负荷没办法继续压。这个“刚性出力”如果不建模调度方案在工程上根本执行不了。所以在模型里垃圾焚烧机组不是普通燃煤机组它更像是一个“必须运行但出力范围受限”的特殊节点。再加上碳捕集设备系统里的电气关系又多了一层碳捕集本身要耗电捕集量越大虚拟电厂提供给电网的净出力就越低。如果目标函数里只算机组发电收益不把碳捕集能耗和碳排放成本放进去优化器会把碳捕集率直接拉到零看起来“电量更多、收益更高”实际上碳排放超了碳成本一算可能完全得不偿失。1.2 电转气和碳捕集形成的碳-气闭环电转气Power to GasP2G简单说就是用电力电解水制氢再和二氧化碳甲烷化生成合成天然气。这一过程最吸引人的地方在于它能把暂时用不掉的再生能源电力转化成气体储存起来解决风光出力的波动问题。而碳捕集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 决策变量一张表说清楚整个模型的决策变量可以分成三类机组状态类、能量流类和碳流类。变量名含义类型u(k,t)机组k在时刻t的启停状态二进制v(k,t)机组k在时刻t是否从停运转为启动二进制P(k,t)机组k的有功出力连续P_re(t)可再生能源实际出力连续P_p2g(t)电转气消耗的有功功率连续E_cap(t)碳捕集量连续alpha(t)烟气分流比/捕集率连续0~1G_p2g(t)电转气产气量连续G_sto(t)储气罐储气状态SOC连续G_buy(t)从外部网购气量连续注意机组启停状态必须是二进制变量这就决定了整个模型是一个混合整数规划MILP或混合整数二次规划MIQP。如果你有条件用二次目标CPLEX/Gurobi可以直接求解MIQP如果担心求解时间可以先把机组成本函数线性化为分段函数变成MILP之后求解稳定性更高。2.3 核心约束机组、碳捕集和电转气机组模型部分最基本的是出力上下限约束和爬坡约束Pmin(k) * u(k,t) P(k,t) Pmax(k) * u(k,t)这个约束表达的是只有当机组运行时u1出力才能落在上下限内机组停运时u0出力强制为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Σ_{it-MINUP1}^{t} u(i) MINUP * (u(t) - u(t-1)) Σ_{it-MINDOWN1}^{t} (1-u(i)) MINDOWN * (u(t-1) - u(t))这个模型在实际代码里可以简化不一定每个机组都强制建模但垃圾焚烧机组建议至少设一个MINUP24即全天必须保持并网运行。这样最符合垃圾焚烧电厂连续运行的真实特性。碳捕集部分的建模要抓住两个关键量烟气中的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(t1) 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^2bPc 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 核心指标的变化在这样一个算例上典型结果大致是指标场景A无CCS/P2G场景B完整模型日运行成本元约100000约92000弃风弃光电量MWh358系统净碳排放t约260约190燃气机组耗气量MWh4530电转气产气量MWh—22为什么会是这样因为P2G在夜间风电大发、负荷低谷时启动把原本要被弃掉的6到8个小时风电变成天然气储存起来在早高峰时给燃气机组用相当于把弃电“搬了个家”。碳捕集则是在白天垃圾焚烧机组出力较高的时段多捕碳利用捕集后较低的净排放水平减少碳配额购买。虽然碳捕集和P2G本身都消耗电能但在算例的总成本账上省下的碳费用和减少的弃电惩罚超过了它们的运行成本所以整体经济性反而更好。4.3 结果里值得留意的现象第一个现象是垃圾焚烧机组在夜间出力始终压在最低技术出力40 MW附近而燃气机组和备用机组完全停运。这很符合焚烧炉连续运行的物理特点也是模型里MINUP24约束起作用的结果。如果把这个约束删掉优化器可能会给出夜间停运垃圾焚烧机组、白天再启动的“数学最优”方案但现实中根本没法执行。第二个现象是碳捕集率并不像很多人想的“越高越好”。即便碳价给到了60元/吨捕集率也没有无限逼近1而是停在一个相对合理的水平。原因很简单捕集率升高会让碳捕集电耗大幅上升进而占用电负荷空间甚至导致系统需要用高价购电来维持平衡。模型在碳成本和电成本之间自动找平衡点这正是优化调度的意义。5. 实际调参中踩过的坑与排查方法5.1 模型怎么突然无解了这是我遇到最多的问题。无解往往不是因为目标函数写错而是约束之间出现了物理上不可能的矛盾。最常见的有三种现象可能原因排查思路求解器返回Primal infeasible爬坡约束与启停状态矛盾比如机组刚启动就要求出力从0跳到60 MW用sol.info定位逐条注释掉约束找到冲突的锚点储气SOC出现NaN储气初值没有设导致递推式在第1小时就出现未定义变量检查Gsto(1)是否赋值二进制变量维度和出力维度不一致P(:,k)和u(:,k)总是必须在同一时间维度用size()检查每个变量的维度还有一个隐蔽问题垃圾焚烧机组的最低出力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产气储气、燃气机组消耗气体它们之间环环相扣漏掉任何一个环节优化结果就会给出看似省钱、实则无法落地的方案。我建议刚开始接触这类模型的朋友先不要急着堆设备数量把一套最简单的“电-气-碳”闭环跑通再逐步加约束换场景。跑通之后你自然就明白下一步该在哪儿加随机性、在哪儿加设备模型了。 虚拟电厂电-气-碳联合优化调度的Matlab建模实现
这两年做低碳电力系统相关的网格优化最常被问到的一个问题就是虚拟电厂调度里同时塞进电转气、碳捕集和垃圾焚烧到底该怎么建模单看任何一个设备网上的资料都不算少但要在一套调度模型里把“电、气、碳”三条线串起来很多朋友一开始会发懵。这篇文章我不打算写那种“从原理到展望”的综述而是直接以一个我自己用Matlab实现的虚拟电厂优化调度项目为底子把变量怎么定、约束怎么写、求解器怎么调、算例怎么安排以及我实际踩过的坑按顺序摊开来讲。适合正在做综合能源系统优化、低碳调度、碳捕集与封存建模的研究生以及想用Matlab把调度模型落到代码里的工程师。1. 项目背景为什么要把电转气、碳捕集和垃圾焚烧放在同一个算盘里1.1 虚拟电厂不再只是“风光储”的拼盘传统的虚拟电厂调度大概就是聚合分布式光伏、风电、储能和备用柴油机组以最小运行成本为目标去优化各设备的出力。这种模型已经很成熟核心约束也只有功率平衡、机组出力上下限、储能SOC等。可一旦加入垃圾焚烧机组问题立刻不一样。垃圾焚烧机组的特点是“不能想停就停”。焚烧炉需要维持稳定燃烧垃圾仓里的垃圾每天都要处理城市固废的连续清运要求决定了焚烧机组必须连续运行。在夜间的低负荷时段其他机组可以大幅降出力甚至停运垃圾焚烧机组往往只能压到最低稳燃负荷没办法继续压。这个“刚性出力”如果不建模调度方案在工程上根本执行不了。所以在模型里垃圾焚烧机组不是普通燃煤机组它更像是一个“必须运行但出力范围受限”的特殊节点。再加上碳捕集设备系统里的电气关系又多了一层碳捕集本身要耗电捕集量越大虚拟电厂提供给电网的净出力就越低。如果目标函数里只算机组发电收益不把碳捕集能耗和碳排放成本放进去优化器会把碳捕集率直接拉到零看起来“电量更多、收益更高”实际上碳排放超了碳成本一算可能完全得不偿失。1.2 电转气和碳捕集形成的碳-气闭环电转气Power to GasP2G简单说就是用电力电解水制氢再和二氧化碳甲烷化生成合成天然气。这一过程最吸引人的地方在于它能把暂时用不掉的再生能源电力转化成气体储存起来解决风光出力的波动问题。而碳捕集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 决策变量一张表说清楚整个模型的决策变量可以分成三类机组状态类、能量流类和碳流类。变量名含义类型u(k,t)机组k在时刻t的启停状态二进制v(k,t)机组k在时刻t是否从停运转为启动二进制P(k,t)机组k的有功出力连续P_re(t)可再生能源实际出力连续P_p2g(t)电转气消耗的有功功率连续E_cap(t)碳捕集量连续alpha(t)烟气分流比/捕集率连续0~1G_p2g(t)电转气产气量连续G_sto(t)储气罐储气状态SOC连续G_buy(t)从外部网购气量连续注意机组启停状态必须是二进制变量这就决定了整个模型是一个混合整数规划MILP或混合整数二次规划MIQP。如果你有条件用二次目标CPLEX/Gurobi可以直接求解MIQP如果担心求解时间可以先把机组成本函数线性化为分段函数变成MILP之后求解稳定性更高。2.3 核心约束机组、碳捕集和电转气机组模型部分最基本的是出力上下限约束和爬坡约束Pmin(k) * u(k,t) P(k,t) Pmax(k) * u(k,t)这个约束表达的是只有当机组运行时u1出力才能落在上下限内机组停运时u0出力强制为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Σ_{it-MINUP1}^{t} u(i) MINUP * (u(t) - u(t-1)) Σ_{it-MINDOWN1}^{t} (1-u(i)) MINDOWN * (u(t-1) - u(t))这个模型在实际代码里可以简化不一定每个机组都强制建模但垃圾焚烧机组建议至少设一个MINUP24即全天必须保持并网运行。这样最符合垃圾焚烧电厂连续运行的真实特性。碳捕集部分的建模要抓住两个关键量烟气中的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(t1) 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^2bPc 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 核心指标的变化在这样一个算例上典型结果大致是指标场景A无CCS/P2G场景B完整模型日运行成本元约100000约92000弃风弃光电量MWh358系统净碳排放t约260约190燃气机组耗气量MWh4530电转气产气量MWh—22为什么会是这样因为P2G在夜间风电大发、负荷低谷时启动把原本要被弃掉的6到8个小时风电变成天然气储存起来在早高峰时给燃气机组用相当于把弃电“搬了个家”。碳捕集则是在白天垃圾焚烧机组出力较高的时段多捕碳利用捕集后较低的净排放水平减少碳配额购买。虽然碳捕集和P2G本身都消耗电能但在算例的总成本账上省下的碳费用和减少的弃电惩罚超过了它们的运行成本所以整体经济性反而更好。4.3 结果里值得留意的现象第一个现象是垃圾焚烧机组在夜间出力始终压在最低技术出力40 MW附近而燃气机组和备用机组完全停运。这很符合焚烧炉连续运行的物理特点也是模型里MINUP24约束起作用的结果。如果把这个约束删掉优化器可能会给出夜间停运垃圾焚烧机组、白天再启动的“数学最优”方案但现实中根本没法执行。第二个现象是碳捕集率并不像很多人想的“越高越好”。即便碳价给到了60元/吨捕集率也没有无限逼近1而是停在一个相对合理的水平。原因很简单捕集率升高会让碳捕集电耗大幅上升进而占用电负荷空间甚至导致系统需要用高价购电来维持平衡。模型在碳成本和电成本之间自动找平衡点这正是优化调度的意义。5. 实际调参中踩过的坑与排查方法5.1 模型怎么突然无解了这是我遇到最多的问题。无解往往不是因为目标函数写错而是约束之间出现了物理上不可能的矛盾。最常见的有三种现象可能原因排查思路求解器返回Primal infeasible爬坡约束与启停状态矛盾比如机组刚启动就要求出力从0跳到60 MW用sol.info定位逐条注释掉约束找到冲突的锚点储气SOC出现NaN储气初值没有设导致递推式在第1小时就出现未定义变量检查Gsto(1)是否赋值二进制变量维度和出力维度不一致P(:,k)和u(:,k)总是必须在同一时间维度用size()检查每个变量的维度还有一个隐蔽问题垃圾焚烧机组的最低出力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产气储气、燃气机组消耗气体它们之间环环相扣漏掉任何一个环节优化结果就会给出看似省钱、实则无法落地的方案。我建议刚开始接触这类模型的朋友先不要急着堆设备数量把一套最简单的“电-气-碳”闭环跑通再逐步加约束换场景。跑通之后你自然就明白下一步该在哪儿加随机性、在哪儿加设备模型了。