先说个现象。我这两年做综合能源系统优化方向的课题最常被问的不是“怎么把模型写出来”而是“为什么我做出来的调度结果没什么意思”。尤其是带热电联产的园区燃气轮机一台电和热就彻底绑死了热负荷一起来机组就得多发电赶上电价低谷时段多余的电要么低价卖掉要么干脆弃掉。后来我把阶梯式碳交易机制和电制氢塞进同一个MATLAB优化模型里调度结果突然就“活了”——碳价把多发电变成实打实的成本压力电制氢又给多余的电找到了消纳出路一压一拉原本僵死的热电耦合被撬开了。这篇文章就从一个小而完整的算例展开讲清楚三件事第一阶梯式碳交易机制为什么比传统单碳价更贴近现实它的分段成本函数怎么建模第二电制氢在综合能源系统里到底承担了什么角色设备模型怎么抽象第三用MATLABYALMIP把这个热电优化MILP问题写出来整个代码骨架和求解思路是怎么组织的。适合正在写论文需要算例支撑、或者刚接触IES优化调度但不想看一堆复杂推导的读者。1. 为什么是“阶梯式碳交易”单一碳价、配额差与分段惩罚的建模演进1.1 单一碳价为什么不够用综合能源系统碳排放核算的最基础套路是把所有碳排放源折算成一个总排放量E给一个无偿配额E0如果排超了就按固定碳价c购买配额没排超就按同一价格卖出。这个模型写起来非常干净C_carbon c * (E - E0)但真拿去调度就会发现一个问题当碳价是常数时碳排放成本本质上只是一个“影子价格”对系统整体的减排压力是线性的、均匀的。也就是说不管超排10吨还是100吨单位成本都一样。现实中的碳交易市场基本不会这么温柔排放越多监管压力越大超排部分通常是分段累进计价类似阶梯电价基础档便宜超标越多单价越贵。这样的非线性惩罚才能把“高碳路径”真正挤出去逼着调度去寻找燃气轮机以外的灵活性资源。1.2 阶梯式碳交易的成本函数表达式我在模型里采用最常见的三段式阶梯碳价设免费配额为E0实际排放为E净买入配额Q E - E0。当Q 0时分成三档第一档0 ≤ Q1 ≤ L1碳价c1第二档L1 ≤ Q2 ≤ L2碳价c2且c2 c1第三档Q3 ≥ L2碳价c3且c3 c2。如果Q 0说明配额有富余允许以较低的价格c_sell出售c_sell一般低于买入档位。于是总碳交易成本写成C_carbon c1 * Q1 c2 * Q2 c3 * Q3 - c_sell * Q_sell这个分段函数在优化模型里属于经典的分段线性函数。它带来的直接效果是当系统排放逼近某个阶梯边界时继续增加排放的边际成本突然跳升调度策略就会发生不连续的切换。这个切换在结果上非常好看也很有论文故事感——它反映了碳市场对“临界排放行为”的抑制。1.3 阶梯区间怎么设为什么是三档三档不是拍脑袋定的。档位太少比如只有两档对中等超排段的约束不够细腻档位太多会增加整数变量直接拉长MATLAB的求解时间。实际算例里我常用L1 100 tCO2L2 300 tCO2c1 30 元/tc2 60 元/tc3 100 元/t当然这些参数和碳市场实际价格、免费配额核发比例都要做灵敏度分析。你完全可以把档位换成4档或5档模型结构不变只是多两个连续变量和区间边界。阶梯碳价还有一个隐藏优势它能天然形成对“配额过度买入”的预防因为低价档容量有限系统不会无脑买配额会更愿意把剩余电量拿去电解制氢从源头降低排放。2. 电制氢在系统里到底扮演什么角色热电耦合的软刀子2.1 以热定电问题的根源燃气轮机热电联产机组有一个绕不开的运行约束在抽气式或背压式结构下发电功率和供热功率存在强耦合。简化建模时常用固定热电比H_chp k * P_chp这个式子看起来简单但它蕴含了一个残酷事实热负荷确定了电出力基本也被锁死。夜间热负荷高、电价低、风电也高的时候CHP机组仍然被迫满发风电就被挤掉。传统系统只能靠电锅炉储热或者直接弃风来缓解这两种做法都不够理想——电锅炉本质上是把电能变成低品味热能量等级掉了弃风则是干脆扔掉零碳电量。2.2 电解槽加储氢给系统装一个灵活负荷槽电制氢设备的本质就是给系统增加了一个“可调节大功率负荷”。电解槽可以在电价低、弃风多的时段吸收多余电能把电能转化为氢的化学能储存起来。在模型里我通常把电解槽功率设为连续变量P_el ∈ [0, P_el_max]电解槽产氢量按效率折算这里我统一用能量单位GJ处理避免kg/h和MWh之间来回换算H2_prod η_el * P_el储氢罐的时序约束是典型的能量库存约束S_h2(t1) S_h2(t) H2_prod(t) - H2_cons(t)并加上容量边界0 ≤ S_h2(t) ≤ S_h2_max这样电解槽就是一个可上可下的“软件负荷”夜间多充电白天少充电完全由优化器说了算。对系统而言它不像电锅炉那样只能服务供热而是把余电变成了一种可跨时段调配的二次能源对碳视角而言消纳风电和光伏制氢本身就是零碳路径。2.3 氢的三条能量出路储氢罐里的氢不能只停留在罐子里模型里要给它设计出口氢燃料电池发电P_fc H2_cons * η_fc_e补峰时段输出电能缓解CHP的电出力压力。氢燃料电池余热回收H_fc P_fc * r_h2h供一部分热负荷直接替代部分CHP供热降低燃料消耗和碳排放。作为氢负荷外售如果园区有加氢站或工业氢负荷可以直接把氢作为产品输出这部分在目标函数里体现为售氢收益。有了这三条出路氢系统就不再是“摆着看的设备”而是真正打通了电、热、气三个能源网络。我自己的习惯是至少保留前两条路否则电解槽消纳的氢没有去处模型很容易出现“制了氢但没用”的冗余设备问题。3. 综合能源系统热电优化的数学模型从能量流到目标函数3.1 先画清楚能流拓扑写代码之前一定要把系统拓扑用表格列出来否则MATLAB里变量一多必乱。我的算例采用以下设备集合设备/环节输入输出说明燃气轮机CHP天然气电、热固定热电比核心产热设备燃气锅炉天然气热备用供热启停灵活电锅炉电热电转热消纳低谷电电解槽电氢电制氢灵活负荷储氢罐氢氢跨时段转移氢能氢燃料电池氢电、热氢的综合利用出口储热罐热热解耦供热与发电风光机组自然电零碳电源按预测曲线给定这个配置不算复杂但足够展示阶梯碳价和电制氢的互动CHP负责基础荷电锅炉和氢燃料电池负责调峰电解槽负责吃多余电。3.2 设备约束与能量平衡约束首先写电功率平衡Pbuy(t) P_chp(t) P_fc(t) P_wind(t) P_pv(t) P_load(t) P_eb(t) P_el(t)每一项都必须量化纲统一我统一用MW。热功率平衡H_chp(t) H_eb(t) H_fc(t) H_tes(t) H_load(t)其中H_tes(t)是储热罐的放热功率正值放热、负值蓄热。CHP运行约束除了热电比还要限制出力上下限和爬坡速率P_chp ∈ [P_chp_min, P_chp_max]-ramp_down ≤ P_chp(t1) - P_chp(t) ≤ ramp_up燃气锅炉同理。电解槽和燃料电池也有出力上下限。最后是储热罐S_tes(t1) S_tes(t) η_tes_c * H_tes_in(t) - H_tes_out(t) / η_tes_dS_tes ∈ [0, S_tes_max]这些约束本身不难难的是别漏掉“设备在同一个时刻既蓄又放”的无效循环。我建议给储热罐的进、出热单独设置二元变量互斥或者在目标函数里加很小的运行维护成本让优化器自动避免无意义循环。3.3 目标函数成本怎么拆目标是最小化系统总运行成本我把它拆成六块购电成本sum(price_buy(t) * Pbuy(t))购气成本sum(price_gas * V_gas(t))其中V_gas(t)是CHP和燃气锅炉的天然气总消耗量单位统一折算为MWh热值运维成本各设备出力乘以单位运维系数碳交易成本上一章写的阶梯碳价函数弃风弃光惩罚penalty * (P_wind_max - P_wind_use)这一项保证优化器优先消纳风光售氢收益price_h2 * H2_sold(t)作为负成本进入目标整体目标函数写法就是把这些项加总。碳交易成本在目标里占比未必很大但它起到的是“方向调节”作用当碳价档位升高时优化器会主动减少外购电、减少燃气消耗把更多电量导入电解槽——这正是我们想要的低碳调度行为。4. MATLAB代码骨架YALMIP下的变量、约束与阶梯碳价线性化4.1 变量定义与维度规划我习惯先把所有变量按“设备时间维”定义T取24也就是一天的调度。统一用MW作为功率单位用MWh作为能量单位碳排放用tCO2。核心变量如下T 24; % 功率变量都是1行T列的sdpvar P_chp sdpvar(1, T); % CHP发电 H_chp sdpvar(1, T); % CHP供热 P_eb sdpvar(1, T); % 电锅炉耗电 H_eb sdpvar(1, T); % 电锅炉产热 P_el sdpvar(1, T); % 电解槽耗电 P_fc sdpvar(1, T); % 氢燃料电池发电 H_fc sdpvar(1, T); % 氢燃料电池供热 Pbuy sdpvar(1, T); % 外购电 V_gas sdpvar(1, T); % 天然气输入热功率 % 储能变量T1维给初始状态留位置 S_h2 sdpvar(1, T1); % 储氢罐容量 S_tes sdpvar(1, T1); % 储热罐容量储能变量做成T1维是我的个人习惯S(1)是初始容量S(t1)是时段t结束时的容量写循环时下标不容易乱。4.2 阶梯碳价的分段线性化两个经典写法把碳交易成本写入目标函数时直接在约束里写max(E-E0, 0)是不可行的YALMIP虽然能处理epigraph形式但对MILP来说最好还是手工展开。第一种写法是基于凸性的简化E_emis sdpvar(1, 1); % 总碳排放 E_quota 500; % 免费配额tCO2 e_plus sdpvar(1, 1); % 净买入配额 e_minus sdpvar(1, 1); % 卖出配额 x1 sdpvar(1, 1); x2 sdpvar(1, 1); x3 sdpvar(1, 1); L1 100; L2 300; c1 30; c2 60; c3 100; c_sell 20; Constraints [Constraints, E_emis sum(emission_grid * Pbuy) sum(emission_gas * V_gas)]; Constraints [Constraints, E_emis - E_quota e_plus - e_minus]; Constraints [Constraints, e_plus x1 x2 x3]; Constraints [Constraints, 0 x1 L1]; Constraints [Constraints, 0 x2 L2 - L1]; Constraints [Constraints, x3 0]; Constraints [Constraints, e_plus 0, e_minus 0]; % 关键为了防止e_plus和e_minus同时非零造成虚假交易加一个互斥二元变量 z_sgn binvar(1, 1); M_big 2000; Constraints [Constraints, e_plus M_big * z_sgn]; Constraints [Constraints, e_minus M_big * (1 - z_sgn)];因为目标函数里碳价是递增的所以低价区间会先被填满不需要额外逻辑约束来保证“先买满第一档再买第二档”。如果你做的是最大化碳交易收益模型这个性质就不成立了此时才需要为每个档位加二元变量和顺序约束。4.3 约束与目标拼装接下来是能量平衡和设备约束。这部分我把循环结构直接展开成矩阵约束Constraints []; % 电功率平衡 Constraints [Constraints, Pbuy P_chp P_fc P_wind P_pv P_load P_eb P_el]; % 热功率平衡 Constraints [Constraints, H_chp H_eb H_fc (S_tes(2:T1) - S_tes(1:T)) H_load]; % CHP固定热电比 出力上下限 爬坡 Constraints [Constraints, H_chp 1.2 * P_chp]; Constraints [Constraints, P_chp 20, P_chp 150]; Constraints [Constraints, -30 P_chp(2:T) - P_chp(1:T-1) 30]; % 电锅炉 Constraints [Constraints, H_eb 0.98 * P_eb, 0 P_eb 50]; % 电解槽和燃料电池 Constraints [Constraints, 0 P_el 80]; Constraints [Constraints, 0 P_fc 60, 0 H_fc 48]; Constraints [Constraints, H_fc 0.8 * P_fc]; % 天然气消耗 Constraints [Constraints, V_gas P_chp / 0.35 H_gb_up]; % 储氢罐动态 Constraints [Constraints, S_h2(2:T1) S_h2(1:T) 0.7 * P_el - 2.1 * (P_fc H_fc)]; Constraints [Constraints, S_h2 0, S_h2 300, S_h2(1) 50, S_h2(T1) 50];上面P_wind、P_pv、P_load、H_load、H_gb_up都是外部给定的数据向量或变量。这里我写的是简化示意实际你还要包含风电消纳变量和弃风惩罚逻辑上把“可消纳风电”设为变量并约束它不超过预测出力。目标函数这样拼carbon_cost c1 * x1 c2 * x2 c3 * x3 - c_sell * e_minus; objective sum(price_buy .* Pbuy) sum(price_gas .* V_gas) ... sum(om_cost) carbon_cost ... sum(penalty_wind .* (P_wind_pred - P_wind_use)) ... - sum(price_h2 .* H2_sold);注意变量维度要对齐。price_buy我直接给一个24维向量而不是标量这样能模拟分时电价。4.4 求解器配置与BM结果提取YALMIP默认用的求解器往往不支持整数变量必须显式指定CPLEX或Gurobiops sdpsettings(solver, gurobi, gurobi.MIPGap, 0.001, verbose, 1); sol optimize(Constraints, objective, ops); if sol.problem 0 P_chp_opt value(P_chp); P_el_opt value(P_el); carbon_cost_opt value(carbon_cost); else disp(求解失败); endMIPGap设为0.001通常在工程上够用再小的话求解时间会成倍增长。T24的小系统一般几秒就出结果但如果加了跨日耦合、8760小时场景建议把聚合步长拉大或者把阶梯碳价再次松弛为连续函数。5. 算例结果怎么设计、怎么读对比实验与碳价灵敏度5.1 对比场景是论文的命根子跑通模型只是第一步结果分析才是体现价值的地方。我强烈建议至少设计三个场景场景A不考虑任何碳交易成本碳项直接从目标函数去掉。场景B固定线性碳价c 60 元/t超配额按统一单价购买。场景C阶梯式碳价按照前面c130, c260, c3100设置。三个场景下重点观察四个指标总运行成本、总碳排放量、弃风率、电解槽功率曲线。大多数时候你会看到这样的趋势场景C的总成本比场景A略高甚至可能更低因为少买了高价配额同时更积极地消纳风电总碳排放明显下降弃风率下降。这说明阶梯碳价产生的“分段压力”系统通过调整设备出力接住了而不是简单把成本转嫁到用户侧。5.2 灵敏度分析看什么做完基础对比我会继续跑一层灵敏度把L1、L2、c2各拉几个值例如把c2从50改到80看电解槽总用电量怎么变。这块数据才是真故事碳价提高初期电解槽利用率会显著上升因为系统发现与其买高价碳配额不如多制氢替代部分燃气出力但碳价进一步提高后电解槽利用率可能进入平台期因为储氢容量和燃料电池出力上限到顶了。这个饱和效应很有价值它能指引你下一步该扩容储氢罐还是增加燃料电池容量。5.3 结果图表怎么画才不浪费MATLAB画图时别直接用默认折线糊一脸。我常用的三个展示角度功率平衡堆叠图横轴24小时电负荷曲线下面用area堆叠显示CHP、风电、购电、燃料电池各自贡献一眼能看出夜间风电被电解槽消纳的效果。碳交易分段柱状图把三个档位的Q1、Q2、Q3画成柱状配合总排放折线展示阶梯碳价下系统停在哪个档位。灵敏度曲线横轴碳价纵轴弃风率或碳排放用两条曲线对比线性碳价与阶梯碳价。图表的意义不在好看而在让审稿人或验收方快速抓到结论阶梯碳价改变了调度行为电制氢消纳了原本被抛弃的风电。6. 调模型的坑和经验从数据尺度到大M取值6.1 大M取值太小截断太大病态我前面例子M_big 2000是拍脑袋的实际应该根据碳交易量的物理上界来确定。算一下最极端情况一天内全部负荷都由外购电承担购电量乘以电网折碳系数再减去免费配额就是最大可能的净买入碳配额。比如系统最大电负荷200MW一天4800MWh折碳约0.8t/MWh总排放约3840t减去配额500t净买入约3340t。这种情况下M_big取4000以上才安全。但也不能取到100000这种数量级否则整数变量对连续变量的约束在数值上会出现“假松弛”Gurobi原对偶间隙都降不下去。我通常先算一下上界再放大20%。6.2 量纲真能坑死人我早期调试出现过最诡异的错误计算结果里电解槽功率永远在零附近怎么调碳价都没反应。最后发现是电网碳排放因子用了kg/kWh而燃气排放系数用了t/MWh混在同一个表达式里碳交易成本比其他成本小了三个数量级优化器根本不care。建议开局就把所有参数列一张表统一成三个基准单位功率MW、能量MWh、碳排放tCO2。6.3 设备“白嫖”循环储能同时充放加了储能和氢系统之后如果没有互斥或运维成本优化器可能玩出“既充电又放电”的白嫖解因为效率损耗抵消不了系统套利。处理办法有两种第一种是在目标函数里给储能功率加一个很小的单位运维成本比如每MWh 0.5元破坏零成本循环第二种是用二元变量强制充放互斥。第一种实现简单且求解快我优先推荐除非审稿人特别较真。6.4 从24小时向8760小时扩展一天算例只能看趋势要做年度分析或季节对比直接扩展到8760小时会让整数变量爆炸。这时候建议把阶梯碳交易的分段决策变量从“每个时段一个”改成“全天统一一个”——因为碳交易本身是按结算周期统计的并不需要每个时段都判断档位。具体操作排放总量是全天累加结果所以碳交易变量只有一个体现在代码里就是x1,x2,x3都是1×1的标量而不是1×T。这样扩展到长周期时整数变量数量不会随T线性增加模型规模可控。最后说点实际操作中的体会阶梯碳交易加电制氢这个组合最初我以为只是两个热点的生硬拼凑算例跑透之后才明白它们的内在逻辑是互补的碳交易抬高“多发电”的隐性成本电制氢给“多余电”创造消纳通道两者在热电联产系统里形成了一组天然的低碳调节机制。做MATLAB实现时最值得花时间的地方不是设备约束本身而是把碳成本的阶梯结构用MILP形式表达清楚。这一步做对了后面套什么设备、扩到什么规模都很顺。这个模型框架我后来还改过一版加入碳捕集设备与氢负荷联动效果也很自然如果你有兴趣可以从电解槽的启停状态和储氢容量配置开始做文章——这些都是在目前这个骨架上很容易延伸的方向。