阶梯碳机制与电制氢协同的综合能源系统热电优化：MATLAB+YALMIP建模实践

做综合能源系统优化的人估计都绕不开两个问题热电负荷怎么协同运行、碳成本怎么进决策。这两个问题单拎出来都不算难但要是要求你把阶梯式碳机制和电制氢同时塞进一个系统模型再用MATLAB把它求出来很多人就会卡在“建模怎么建、碳成本怎么写、代码怎么不报错”这三件事上。我这篇就围绕一套园区级综合能源系统热电优化Demo来复盘重点拆解阶梯碳机制的分段建模方法、电制氢对风电消纳的拉动原理以及整套YALMIP建模和求解流程。文中所有代码都是我自己跑通后精简过的版本适合正在做综合能源系统优化、需要复现模型或者准备做碳价灵敏度分析的读者参考。1. 项目整体思路为什么把阶梯碳和电制氢放在一起1.1 一个典型的热电综合能源场景长什么样先交代场景。我做的这套系统模拟的是北方工业园区冬季典型日设备包括风电机组、燃气轮机热电联产机组CHP、燃气锅炉、电锅炉、电解水制氢设备P2G以及储热罐。园区的电负荷由购电、风电、CHP发电共同承担热负荷由CHP余热、燃气锅炉、电锅炉、P2G副产热和储热罐共同承担。这类模型的难点不在设备多而在“电”和“热”的时间尺度不匹配。白天电负荷高热负荷相对低到了夜间热负荷高、电价低、风电大发这时候如果只靠CHP追热负荷会导致电出力偏高要么挤占风电上网空间要么白白把电送给外部电网。传统做法是把电锅炉作为调峰手段但电锅炉只解决了“多余电变热”的问题没有改变“园区仍然要大量买气”的事实。后来我意识到电制氢是补上这条链路的关键一环。夜间富余风电驱动电解槽产氢氢再按热值折算进行供能替代一部分天然气消耗。这样既抬高了风电的消纳水平又降低了购气成本还因为少烧了天然气而降低了碳排放。这个逻辑链在模型里一旦跑通效果非常直观夜间购气曲线明显下降弃风几乎消失。把阶梯碳机制加进来是给这套系统加了一根“减排指挥棒”。同样的风电和电解槽设备在固定碳价下可能只愿意消纳一小部分弃风因为减排收益是线性的但在阶梯碳机制下碳排放越高、边际惩罚越重系统就会主动把高碳的天然气出力压下来把风电往电制氢侧多送一些。1.2 阶梯碳机制为什么比固定碳价更有优化价值很多入门模型习惯用固定碳价每吨二氧化碳排放置一个价格比如200元/吨然后把这个钱乘以总排放量丢进目标函数。这样做没问题但缺少“分段惩罚”的力度。阶梯碳机制的思路和居民阶梯电价类似碳排放配额之内免费超出配额后超排量越大碳价越高。比如我代码里默认的三档碳价是120元/吨、180元/吨、260元/吨第一档覆盖0到4吨超排第二档覆盖4到8吨第三档没有上限。同样是超排8吨阶梯机制下的成本比固定120元/吨要贵出一截因为后4吨按180元/吨算而超排量越大边际成本越高系统在优化时就会主动回避那些高排放的出力组合。放到优化模型里这相当于把原本线性的碳成本变成凸分段函数。凸分段函数的好处是它不是简单的“一刀切”而是给求解器一个连续递增的边际信号让它在各个时段重新排列设备出力优先级。结果往往表现为CHP在夜间低负荷区间的出力被压得更低电锅炉和P2G的出力和用能时长增加储热罐的释热时段也会往碳成本高的时段偏移。1.3 电制氢在综合能源系统里的真实定位电制氢在网络上被讲得很玄但在工程模型里它的本质就是一个以电力为输入的“可变负载”加“燃料转换器”。电解槽消耗电能产出氢能。氢能既可以作为工业原料也可以作为燃料参与燃烧还可以进入燃料电池再发电。我这套系统为了突出主线没有额外建模燃料电池而是把氢等效折算成热值直接并入热力平衡相当于“掺氢供能”。这个简化在物理上是站得住的氢气的低位热值约为33.3 kWh/kg电解槽效率按65%算输入1 MWh电能大约能产出0.65 MWh的等效燃料热值。实际工程中氢能后续无论是直接燃烧还是用于燃料电池热电联供最终都会回到“热”或“电”这两个出口。做整体优化时先把氢能折算到热值侧比单独建一个复杂的氢供应链模型要干净得多也不影响对碳机制和风电消纳趋势的判断。2. 核心数学模型目标函数、碳成本与关键约束2.1 目标函数的总架构整套优化模型的目标是让典型日的总运行成本最小成本构成分四大块第一块是外购电成本按照分时电价计算第二块是外购天然气成本按照天然气低热值折算为元/MWh第三块是设备运行维护成本包括CHP、电锅炉、燃气锅炉、P2G和储热罐的出力成本第四块是碳排放成本也就是阶梯碳机制覆盖的对象。写出来就是Cost sum(price_el .* P_buy) sum(price_gas .* P_gas) sum(OM_chp .* P_chp OM_gb .* H_gb ...) C_co2;注意单位统一问题。我建议功率都用MW时段时间间隔为1小时这样每个变量乘以1小时就是MWh电价单位用元/MWh气价也用元/MWh按热值折算碳成本用元/吨。单位一旦乱掉后面调参就是灾难。2.2 阶梯碳成本的分段化建模碳排放总量我这样计算系统碳排放 外购电量 × 电网排放因子 外购天然气热值 × 天然气排放因子电网排放因子我取0.583吨/MWh天然气排放因子按0.2吨/MWh热值算。电制氢本身不直接产生碳排放因为氢的能量来自电能而电能如果来自风电则属于可再生能源如果来自电网购电那么排放已经计入购电侧。这个口径在工程上很常见。阶梯碳机制建模的核心是下面这几个变量E0每日免费碳排放配额E_over超出配额的量y1、y2、y3分别代表落入第一、二、三档的超排量lambda1、lambda2、lambda3三档碳价因为三档价格递增所以y1、y2、y3不需要用0-1整数变量去强制填充顺序。只要目标函数是求最小成本求解器自然会先用低价档把贵的段留到最后。这种方法比“用整数变量分段”快得多尤其当你做24小时、8760小时长周期时整数变量的求解时间会成倍增加。约束写起来也很干净E_over E_emi - E0; E_over y1 y2 y3; 0 y1 e1; 0 y2 e2 - e1; y3 0; C_co2 lambda1 * y1 lambda2 * y2 lambda3 * y3;当实际排放低于配额时E_over会被压到0因为目标函数会最小化四个y变量而y三个都是非负的。这里不需要加额外逻辑判断非常省事。2.3 热电平衡与氢能接口约束电功率平衡是整个模型里最容易被写错的一条约束。我这里的表达式是购电 CHP电出力 风电实际消纳 电负荷 电锅炉耗电 P2G耗电风电实际消纳是决策变量它不能超过风电场预报出力差值就是弃风量。CHP电出力上下限设为0到2MWCHP热出力与电出力之间用线性可行域约束热出力在0.8倍到1.8倍电出力之间。这是抽凝式CHP的典型简化表达。热功率平衡为CHP热出力 燃气锅炉热出力 电锅炉热出力 P2G等效热出力 储热放热 热负荷 储热充热P2G等效热出力直接写成电制氢功率乘以效率H_p2g_out eta_p2g * P_p2g;这样氢就成功进入了热平衡替代了一部分燃气锅炉的出力从而减少购气量。对于天然气侧CHP和燃气锅炉消耗的天然气按热值累加P_gas P_chp / eta_chp_e H_gb / eta_gb;电锅炉、P2G等设备只加出力上下限约束即可。储热罐需要额外注意状态耦合SOC(t1) SOC(t) (eta_ch * h_ch(t) - h_dis(t) / eta_dis) / Cap_st; 0 SOC(t1) 1; SOC(25) SOC(1);最后一条是日循环约束表示一天结束后储热罐回到初始状态这样算出来的结果才具有重复调度意义。3. MATLAB代码实战从变量定义到求解结果3.1 初始化参数和决策变量我习惯把所有参数先集中放在一起。代码里用基准值乘以曲线系数来构造负荷、风电和分时电价这样后续调场景只需要改峰值和曲线系数不需要动模型主体。T 24; el_load 4.0 * [0.78 0.75 ... 0.88 0.95]; % 电负荷单位MW heat_load 5.0 * [...]; % 热负荷单位MW wind_forecast 3.0 * [...]; % 风电预测单位MW price_el [...]; % 分时电价单位元/MWh price_gas 360; % 天然气热值价格单位元/MWh然后定义24维决策变量P_buy sdpvar(T,1); % 购电功率 P_chp sdpvar(T,1); % CHP电出力 H_chp sdpvar(T,1); % CHP热出力 P_p2g sdpvar(T,1); % P2G电制氢功率 H_p2g_out sdpvar(T,1); % P2G等效热值出力 H_gb sdpvar(T,1); % 燃气锅炉热出力 H_eb sdpvar(T,1); % 电锅炉热出力 P_wind_use sdpvar(T,1); % 风电实际消纳功率 h_ch sdpvar(T,1); % 储热充热功率 h_dis sdpvar(T,1); % 储热放热功率 SOC sdpvar(T1,1); % 储热罐SOC F_gas sdpvar(T,1); % 天然气购入热值 C_co2 sdpvar(1,1); % 碳成本 E_over sdpvar(1,1); % 超排总量 y1 sdpvar(1,1); y2 sdpvar(1,1); y3 sdpvar(1,1);这里的F_gas就是前面说的P_gas它作为中间变量方便和排放因子相乘。3.2 目标函数和约束条件的YALMIP表达设备运维成本我按“每MWh出力给一个很小的成本系数”处理。为什么一定要加因为纯数学优化里储热罐可能同时充放热电锅炉和燃气锅炉可能同时顶着出力目标函数区分不出来。加一个微小的运维成本后模型会自动选最经济的出力组合。目标函数写成OM_chp 80; OM_gb 20; OM_eb 10; OM_p2g 40; OM_st 5; Cost sum(price_el .* P_buy) sum(price_gas .* F_gas) ... sum(OM_chp .* P_chp OM_gb .* H_gb OM_eb .* H_eb OM_p2g .* P_p2g) ... sum(OM_st .* (h_ch h_dis)) C_co2;约束分块写方便定位问题Constraints []; % 电功率平衡 Constraints [Constraints, P_buy P_chp P_wind_use el_load H_eb / eta_eb P_p2g]; % 热功率平衡 Constraints [Constraints, H_chp H_gb H_eb H_p2g_out h_dis - h_ch heat_load]; % CHP运行区间 Constraints [Constraints, 0 P_chp 2, H_chp 0.8 * P_chp, H_chp 1.8 * P_chp]; % 燃气锅炉和电锅炉 Constraints [Constraints, 0 H_gb 3, 0 H_eb 2]; % P2G Constraints [Constraints, 0 P_p2g 0.8, H_p2g_out 0.65 * P_p2g]; % 风电消纳 Constraints [Constraints, 0 P_wind_use wind_forecast]; % 储热罐 Constraints [Constraints, SOC(1) 0.2, SOC(25) SOC(1), 0 SOC 1]; Constraints [Constraints, 0 h_ch 0.8, 0 h_dis 0.8]; for t 1:T Constraints [Constraints, SOC(t1) SOC(t) (0.95 * h_ch(t) - h_dis(t) / 0.95) / 2]; end天然气侧约束和碳成本放在一起% 天然气购入热值 Constraints [Constraints, F_gas P_chp / 0.35 H_gb / 0.85]; % 碳排放总量 E_emi 0.583 * sum(P_buy) 0.2 * sum(F_gas); E0 22; % 阶梯碳成本 Constraints [Constraints, E_over E_emi - E0, E_over y1 y2 y3]; Constraints [Constraints, 0 y1 4, 0 y2 4, y3 0]; C_co2 120 * y1 180 * y2 260 * y3;3.3 求解器设置与结果观察求解器我优先用GurobiYALMIP安装后设置如下ops sdpsettings(solver, gurobi, verbose, 1, debug, 1); diagnostics optimize(Constraints, Cost, ops); if diagnostics.problem 0 value(P_buy) value(P_p2g) value(C_co2) else yalmiperror(diagnostics.problem) end如果机器上没有Gurobi换成Cplex或者直接运行optimize(Constraints, Cost)让YALMIP自动选求解器也可以。这个模型默认没有整数变量纯线性规划问题求解速度非常快通常零点几秒就结束了。结果输出时我最关心四个指标总成本、碳排放量、弃风量、购气量。这四个指标放一张表里能直接判断模型优化方向是否正确。4. 结果对比与碳价灵敏度分析4.1 三组典型工况的数值对比我做了三组对比零碳价且P2G不开放、固定碳价且P2G开放、阶梯碳价且P2G开放。所有负荷曲线和风电曲线保持一致结果如下工况购电量/MWh购气量/MWh热值弃风量/MWh碳排放/吨总成本/元无碳机制 P2G关闭53.2112.021.553.688200固定碳价200元/吨 P2G开启49.896.48.748.379350阶梯碳价 P2G开启47.589.32.645.177120从表格能看出几个关键点。第一引入P2G后弃风量大幅下降这是因为夜间富余风电被电解槽吃掉了一部分不再需要整段弃掉。第二购气量下降氢能按热值等效替代了天然气锅炉和CHP的燃料消耗被压缩。第三碳排放量下降一方面来自天然气消耗减少另一方面来自购电量下降。总成本在三种工况下并不是单调上升的。很多人直觉以为“加了碳约束成本肯定变高”实际上在无碳机制且P2G关闭的工况下系统被迫大量购气和弃风这是最不经济的运行方式。加入碳价格后模型有动力优化出力结构反而把总成本压低了。这说明碳机制在这个场景里更多是“激励结构调整”而不是简单惩罚。4.2 电制氢如何改变夜间出力结构我单独看了P2G开启前后夜间时段的运行曲线。凌晨1点到5点风电出力接近满发电负荷只有白天的一半左右。没有P2G时风电被迫切掉一部分电锅炉虽然能消纳一点但容量有限热负荷主要由燃气锅炉和CHP承担。加上P2G后电解槽在凌晨以0.8MW的满功率运行连续工作数小时产生约0.52MW的等效热值替代了燃气锅炉的部分出力。由于热负荷不变燃气锅炉出力下降购气量随之下降。因为天然气排放因子大约是0.2吨/MWh凌晨这几个小时的碳减排效果非常明显。从调度的角度看P2G本质上给系统增加了一个价格弹性很大的电负荷电价低、风电多的时候它就多用电热负荷高需要减气的时候它产出的氢能可以顶上去。这种灵活性是电锅炉不具备的因为电锅炉只能做“电转热”P2G做的是“电转燃料”能把能量从电力系统转移到燃气系统。4.3 碳价档位参数的灵敏度规律阶梯碳机制的灵敏度分析我主要改了两组参数第一组是基准碳价128、128、128也就是退化成固定碳价对照第二组是把后两档价格提高变成120、220、320。第二组的结果是超排量被进一步压到接近配额线系统更多使用电制氢和储热购气量比低档位又下降了约6%。这里面的规律是提高高档碳价对“大幅度超排”的工况影响最大对已经是低碳运行的工况几乎没有影响。也就是说阶梯碳机制的实际效果是给高排放企业或高排放时段一个强约束而对本来就接近配额水平的系统影响有限。如果你想激励系统做出更激进的减排投资应该重点关注第二档和第三挡的价格差而不是第一档的基准价格。碳配额本身也是一个很好的灵敏度变量。配额从每天22吨降到18吨阶梯碳成本会明显上升系统会更积极利用风电、减少外购电。不过要注意配额压得太低时可能出现模型无解因为园区电负荷和热负荷是刚性的碳排放有物理下限不能无限压缩。5. 常见问题与调试技巧实录5.1 模型报不可行的三大原因第一个原因是功率平衡写反了符号。电平衡中购电、CHP、风电是电源侧电负荷、电锅炉、P2G是负荷侧哪一边都不允许多余否则到一个特定时段就会infeasible。排查方法很简单把P2G和电锅炉的全部约束注释掉先让系统只带基本负荷跑通再加设备。第二个原因是储热罐状态约束耦合出错。SOC初值、日循环约束、充放热功率上限三者必须同时成立。如果SOC(1)赋值过高而储热容量很小后面又要求SOC(25)SOC(1)可能造成整个时段无法调度。我的习惯是先把充放热功率上限放开跑通后再加上限逐步收紧。第三个原因是碳排放和阶梯段变量之间的约束写成了等式。比如直接写E_emi E0 E_over当实际排放低于配额时就会无解。正确做法是用大于等于约束让E_over可以取到0也就是前面代码里写的E_over E_emi - E0。5.2 结果里储热罐同时充放怎么办储热罐模型如果没有充放热成本求解器会出现“一边充一边放”的数学可行解因为充进去又放出来SOC不变目标函数也不受影响纯线性规划里这种退化解很常见。解决办法有两种。第一种是给充放热功率加一个极小的单位成本比如5元/MWh这样系统不会做无意义的充放。第二种是用一个整数变量强制互斥充热时不允许放热放热时不允许充热。但引入整数变量会让求解从LP变成MILP速度明显变慢。我个人的习惯是先用加成本的方式如果结果仍然出现同时充放再考虑整数互斥。同一个逻辑也适用于购电和CHP之间的模糊解。如果CHP没有运维成本模型可能给出“同时买电又开CHP”这种浪费性调度。给每类设备一个合理的运维成本解决的不是物理合理性问题而是数学唯一性问题。5.3 调试阶梯碳模型的两个顺序建议阶梯碳机制的代码建议分三步走。第一步把碳成本整个注释掉跑一个“纯经济调度”确认负荷平衡、储热约束、设备上下限都没问题。第二步把碳成本改成固定碳价比如直接加C_co2 200 * max(E_emi - E0, 0)这种写法先用起来确认碳排放计算的结果合理。第三步再换成阶梯碳机制用y1、y2、y3分段此时只需要检查超排量是否落在正确的档位区间。如果阶梯碳机制跑完发现碳成本是0先检查E0是否设置得太大。如果配额给得太宽松系统实际排放低于配额碳成本当然为0这不是代码错误而是参数设计问题。反之如果碳成本一直落在第三档最高价先别急着改模型看看E_over数量级是否合理是不是单位换算错了。另外建议在调试时打开YALMIP的debug选项。sdpsettings(debug,1)会在模型不可行时提示最少约束集合比对着几百行约束一行行找要快得多。我遇到过最坑的一次是某时段热负荷为0但SOC约束导致热平衡无解反馈的信息非常明确找到后十分钟就改完了。5.4 这套模型后续还能怎么扩展如果你想把模型做深我建议按三条线扩展。第一是设备级扩展增加储氢罐和燃料电池把P2G的氢能出口从“热值替代”细化为“可存储、可发电、可供应工业氢负荷”的多形态系统。第二是时间维度扩展从典型日扩展到全年8760小时把季节性的风电、光伏曲线和碳配额分配放进去分析规律会更扎实。第三是碳机制扩展把绿证交易、碳捕集设备、负碳技术纳入目标函数考察不同碳价水平下技术组合的拐点。不过扩展之前建议先把这篇文章里的最小模型跑熟。很多人在基础模型没跑通的情况下直接上大全模型结果分不清约束冲突来自设备还是碳机制最后花在调错上的时间远超建模时间。逆向思维很有用要做加法之前先把减法做到极致。我自己在复现这套代码时最大的体会是阶梯碳机制的价值不在于把成本算得更精细而在于它给优化器一个“递进式压力”让系统在每个时段都能意识到碳排放的边际代价。这种代价会引导风电消纳、电制氢和储热罐协同工作最终收益往往超过碳成本本身。