
简介面向综合能源系统优化与碳交易研究人群提供一套考虑阶梯式碳交易机制与电制氢的热电联产优化计算方案。压缩包内共2个文件包括计算书PDF和MATLAB主程序m文件整体仅2.41MB便于快速下载与本地复现。内容聚焦阶梯式碳交易机制下碳价分段定价约束、电制氢过程的能量转换效率以及热电联产系统的多能协调调度涵盖线性规划、非线性规划、多目标优化等建模与求解思路并给出碳交易价格与碳排放限制共同约束下的系统最优运行策略。读者可对照计算书逐步理解模型假设、约束条件与推导细节运行主程序即可获得热电出力、购能成本、碳排放量等关键结果用于指导企业生产计划编排与减排方案设计。目前已有50人学习适合能源电力专业研究生、工程技术人员以及对碳交易机制感兴趣的开发者参考。1. 为什么热电优化要把阶梯碳交易和电制氢放进同一个模型综合能源系统的热电联供调度过去是「电热平衡燃料成本」的线性规划加入碳交易后问题结构变了。阶梯式碳交易机制里免费碳配额是固定值实际排放超过配额的量会被分成几段段位越高边际碳价越贵。也就是说碳成本不再是常数调度必须在“少烧气、多买电”和“多产热、带出电”之间动态取舍或者切换到电制氢来转移负荷与碳排放。电制氢P2H不是简单加一个负荷电解槽耗电产氢同时回收余热供热是电氢热三重耦合的可调设备能吸纳富余风电也能用余热替代燃气锅炉。下面给出一个包含阶梯碳交易、热电联产和P2H的综合能源系统热电优化模型配套可直接运行的 Python 示例代码并说明参数怎么调、结果怎么验证。适合正在搭调度模型、写论文算例或做方案比选的工程师。2. 热电联产与P2H的物理模型能量流和碳排放从哪来2.1 CHP机组热电可行域的简化建模热电联产机组并不是任意电出力和热出力组合都能稳定运行。背压式机组的热电比基本固定抽汽式机组可以主动调节抽汽量让电、热出力在一个多边形区域内变化。做优化调度时常见做法是用一组线性不等式把可运行区域包出来避免把机组简化成“固定热电比”而丢掉调节自由度。2.1.1 可行域的线性不等式对于抽汽式机组我一般用三条边界电出力P_gt的上下限P_gt_min ≤ P_gt ≤ P_gt_max热出力H_gt的上下限0 ≤ H_gt ≤ H_gt_max电热耦合约束α * P_gt ≤ H_gt ≤ β * P_gtα、β是热容比边界由机组热力特性曲线标定。比如某台 120 MW 级抽汽机组可取α0.5、β1.2意思是热出力不能低于电出力的 50%也不能高于 120%。这个可行域在碳交易场景下非常关键当碳价落入高阶梯时调度程序会倾向于降低电出力来减少燃烧排放同时让热负荷由燃气锅炉或 P2H 余热分担。如果模型把热电比固定死这个“降电保热”的调节路径就完全消失了。2.2 电解槽的功率产氢产热平衡P2H 设备在调度模型里的角色是“能量转换器”。电解槽输入电能E_p2h产出氢气按热值计H2_out同时回收一部分电解过程的余热H_p2h。用线性效率模型描述H2_out η_h2 * E_p2h碱性电解槽η_h2一般在 0.650.75H_p2h η_heat * E_p2h可回收余热效率η_heat取 0.150.25三者能量关系必须满足η_h2 η_heat ≤ 0.9其余为散热损失实际项目中电解槽还有运行区间限制低于 30% 额定负载时制氢效率下降明显所以最小负载一般取 10%30%。调度模型里我会加一个爬坡约束|E_p2h(t) - E_p2h(t-1)| ≤ ramp_p2h防止求解结果让电解槽功率在相邻时段跳变过大给现场执行带来困难。2.3 阶梯式碳交易的分段成本函数系统碳排放有三个来源CHP 机组和燃气锅炉的天然气燃烧排放以及外购电力的间接排放。免费配额E_quota可以按历史强度法设定为一个固定值。设实际排放减去配额为Δ E_total - E_quota阶梯式碳交易按Δ落在哪个区间计费Δ ≤ 0配额富余可卖出0 ≤ Δ ≤ Δ_1超出部分按较低碳价p_1Δ_1 ≤ Δ ≤ Δ_2多出的第二部分按p_2Δ Δ_2更大部分按p_3因为算法是求运行成本最小且碳价p_1 p_2 p_3碳排放缺口的购买量会自然按“低价段优先”的顺序填满不需要额外加二进制变量来强制分段顺序。这个凸性判断是减少模型规模的关键。如果采用单碳价模型调度只会比较燃气和购电的边际成本而阶梯碳价会让缺口跨过某个分段点时整体技术组合发生明显跳变这是机制本身带来的结构性差异。2.4 储能与供能设备运行边界综合能源系统里蓄热罐几乎是标配。蓄热设备用一阶动态方程建模S(t1) S(t) P_ch * η_ch - P_dis / η_dis蓄热量S有上下限充放功率也有上限同一时段不能同时充放需要引入 0-1 变量做互斥约束这也是模型成为 MILP 的原因系统里的电锅炉、燃气锅炉各自有爬坡和容量约束。把 CHP、P2H、锅炉、储热放在一起看决策变量规模约 24 个时段 × 15 个变量是一个既能覆盖热电耦合、阶梯碳交易、P2H 余热回收和储能时移又能在普通笔记本上快速求解的最小可行模型。对象关键参数默认取值说明CHP 机组α / β 热容比0.5 / 1.2决定热电调节空间P2H 电解槽η_h2 / η_heat0.70 / 0.20电到氢、电到热效率燃气锅炉η_boiler0.90热效率储热罐η_ch / η_dis0.95 / 0.95充放热效率天然气排放em_gas0.20 tCO2/MWh按燃料热值电网排放em_grid0.40 tCO2/MWh间接排放因子3. 用 Python PuLP 跑通综合能源系统热电优化最小示例代码3.1 数据与参数准备24 小时风电、光伏、电热负荷曲线求解这类问题我常用 PuLP。它开源、内置 CBC 求解器支持 MILP也能导出 LP 文件检查模型结构。下面是一套可以直接复制的 python 代码骨架先造 24 小时示意数据。import numpy as np import pulp from pulp import value T 24 t np.arange(T) # 电负荷、热负荷白天热电夜间热大电小单位 MW load_e np.interp(t, [0, 7, 12, 14, 18, 21, 24], [45, 55, 80, 75, 85, 65, 45]) load_h np.interp(t, [0, 6, 12, 14, 18, 24], [70, 55, 25, 30, 45, 70]) # 风电、光伏风电夜间强光伏只在白天 wind 28 12 * np.cos(2 * np.pi * (t - 1) / 24) pv np.array([0,0,0,0,5,15,30,45,60,65,58,40,22,8,0,0,0,0,0,0,0,0,0,0], dtypefloat) # 分时购电价单位元/MWh price_buy np.array([350]*6 [550]*6 [850]*3 [550]*3 [850]*3 [550]*3, dtypefloat)数据里np.interp负责在给定的时间断点之间做线性插值生成平滑的日电负荷和热负荷曲线。替换成真实 SCADA 数据时注意功率单位统一为 MW时间分辨率与调度时段一致这里按 1 小时一个点处理。3.2 决策变量、燃料成本与目标函数模型的决策变量包括CHP 电热出力、燃气锅炉热出力、P2H 耗电量、外购电量、弃风量和蓄热罐充放热功率。核心成本是燃料、购电、碳成本和 P2H 运维减去售氢收益。prob pulp.LpProblem(IES_Thermo_Hydrogen, pulp.LpMinimize) p_chp pulp.LpVariable.dicts(p_chp, range(T), lowBound20, upBound120) h_chp pulp.LpVariable.dicts(h_chp, range(T), lowBound0, upBound150) h_gb pulp.LpVariable.dicts(h_gb, range(T), lowBound0, upBound80) e_p2h pulp.LpVariable.dicts(e_p2h, range(T), lowBound0, upBound30) buy_e pulp.LpVariable.dicts(buy_e, range(T), lowBound0, upBound200) curtail pulp.LpVariable.dicts(curtail, range(T), lowBound0, upBound50) c1 pulp.LpVariable(c1, lowBound0, upBound80) c2 pulp.LpVariable(c2, lowBound0, upBound40) c3 pulp.LpVariable(c3, lowBound0) gas_price 200 # 元/MWh h2_price 600 # 元/MWh按氢热值计价 om_p2h 20 # 元/MWh curtail_penalty 50 # 元/MWh f_chp {i: 2.4 * p_chp[i] 0.6 * h_chp[i] for i in range(T)} prob ( pulp.lpSum(gas_price * (f_chp[i] h_gb[i] / 0.90) for i in range(T)) pulp.lpSum(price_buy[i] * buy_e[i] for i in range(T)) 50 * c1 70 * c2 120 * c3 pulp.lpSum(om_p2h * e_p2h[i] for i in range(T)) pulp.lpSum(curtail_penalty * curtail[i] for i in range(T)) - pulp.lpSum(h2_price * 0.70 * e_p2h[i] for i in range(T)) )燃料模型2.4 * p_chp 0.6 * h_chp是简化后的燃料热功率表达式意思是每发 1 MW 电约需 2.4 MW 燃料热每供 1 MW 热约需 0.6 MW 额外燃料。真实项目应改用机组热力特性曲线拟合系数。售氢收益放在目标函数里作为负成本单位是“元/MWh 氢热值”这样 P2H 是否运行取决于氢价与电价的竞争关系。3.3 阶梯碳交易约束怎么写前面说过阶梯碳价是递增的模型在最小化成本时天然会先把低价段用完。所以只要把缺口拆成三段连续变量c1 c2 c3并给每段不同单价即可不需要额外排顺序约束。em_gas 0.20 # tCO2/MWh em_grid 0.40 # tCO2/MWh quota_rate 0.15 total_emission ( em_gas * pulp.lpSum(f_chp[i] for i in range(T)) em_gas * pulp.lpSum(h_gb[i] / 0.90 for i in range(T)) em_grid * pulp.lpSum(buy_e[i] for i in range(T)) ) quota quota_rate * (load_e.sum() load_h.sum()) prob c1 c2 c3 total_emission - quota, carbon_gapc1、c2、c3的单位是吨 CO2quota_rate决定配额松紧程度。这个约束写成而不是是因为配额富余时模型可以少买甚至不买但代码里没有实现“卖出配额”的负成本逻辑如果需要卖配额要把total_emission - quota拆成正负两个变量并加互斥约束。3.4 功率平衡与蓄热动态约束电力平衡和热力平衡是每个时段都必须满足的等式约束。风电和光伏作为不可控出力允许通过curtail弃掉一部分curtail在目标函数里有惩罚成本所以只有消纳不了时才会启用。蓄热罐的充放热互补用二进制变量u_st强制。soc pulp.LpVariable.dicts(soc, range(T1), lowBound0, upBound100) ch_st pulp.LpVariable.dicts(ch_st, range(T), lowBound0, upBound30) dis_st pulp.LpVariable.dicts(dis_st, range(T), lowBound0, upBound30) u_st pulp.LpVariable.dicts(u_st, range(T), catBinary) for i in range(T): prob wind[i] - curtail[i] pv[i] p_chp[i] buy_e[i] load_e[i] e_p2h[i], felec_{i} prob h_chp[i] h_gb[i] 0.20 * e_p2h[i] dis_st[i] load_h[i] ch_st[i], fheat_{i} prob h_chp[i] 0.5 * p_chp[i], fchp_low_{i} prob h_chp[i] 1.2 * p_chp[i], fchp_up_{i} prob ch_st[i] 30 * u_st[i], fchmax_{i} prob dis_st[i] 30 * (1 - u_st[i]), fdismax_{i} prob soc[0] 20 prob soc[T] 20 for i in range(T): prob soc[i1] soc[i] 0.95 * ch_st[i] - dis_st[i] / 0.95, fsoc_{i}注意热平衡里0.20 * e_p2h[i]是电解槽回收余热这是 P2H 参与供热的地方。如果不写这一项P2H 在模型里只是个纯电负荷经济价值会明显被低估。蓄热罐的 SOC 在调度首尾都固定为 20 MWh保证一天调度结果不会“偷走”或“残留”蓄热量。3.5 求解输出与结果落盘prob.solve(pulp.PULP_CBC_CMD(msg0, timeLimit60, gapRel0.001)) print(求解状态:, pulp.LpStatus[prob.status]) print(总运行成本(元):, round(value(prob.objective), 2)) import pandas as pd df pd.DataFrame({ 电负荷/MW: load_e, 热负荷/MW: load_h, 风电/MW: wind, 光伏/MW: pv, CHP电/MW: [value(p_chp[i]) for i in range(T)], CHP热/MW: [value(h_chp[i]) for i in range(T)], 锅炉热/MW: [value(h_gb[i]) for i in range(T)], P2H电/MW: [value(e_p2h[i]) for i in range(T)], 购电/MW: [value(buy_e[i]) for i in range(T)], 弃风/MW: [value(curtail[i]) for i in range(T)], 蓄热净充/MW: [value(ch_st[i]) - value(dis_st[i]) for i in range(T)], }) print(df.round(2).to_string(indexFalse)) print(碳交易购买量(t):, round(value(c1) value(c2) value(c3), 2))PULP_CBC_CMD的三个参数值得说明msg0关闭求解器日志避免终端刷屏timeLimit60设置 60 秒求解上限防止大模型卡死gapRel0.001要求相对 MIP 间隙到 0.1%保证解质量。这个规模的问题 CBC 通常在几秒内解完。4. 阶梯碳价与 P2H 容量的敏感性分析参数怎么调4.1 碳价阶梯提高后CHP 出力会发生什么把carbon_price从[50, 70, 120]改成[100, 150, 200]模型的直接反应是压低 CHP 电出力、增加低谷时段购电同时让 P2H 在风电大发时段多运行。原因是 CHP 的碳排放与燃料消耗近似成正比碳价提高后燃气发电的边际成本上升速度比购电更快。算例CHP 日发电量(MWh)购电量(MWh)P2H 用电(MWh)碳交易缺口(t)总成本(元)基准碳价约 720约 540约 260约 135约 612000碳价提高约 645约 585约 305约 118约 668000P2H 容量翻倍约 690约 455约 410约 108约 586000这张表是趋势示意换一套风电光伏数据绝对值会变但方向一致碳价越高系统越倾向于用电替气P2H 容量越大系统越能在不弃风的前提下把多余电量转成氢和热同时减少碳缺口。4.2 P2H 容量上限的边际价值敏感性扫描时把e_p2h的上界从 10 MW 一路加到 80 MW每次重新求解记录总成本和弃风量。一般会看到“收益先陡后平”的曲线P2H 容量从 0 加到 30 MW 时弃风量和购电成本快速下降超过某个阈值后新增容量大部分时间闲置边际收益趋近于零。这个阈值大致等于“低谷时段风电盈余最大值”。判断方法很简单把风电出力曲线减去不可调负荷后看夜间最大盈余是多少P2H 容量设在盈余的 1.21.5 倍就基本够用再大就是浪费投资。4.3 三个容易忽略的参数边界配额覆盖率quota_rate设置过高比如大于 0.4 时系统几乎没有碳压力阶梯碳价的三段结构不会被激活敏感性分析全部失真。P2H 效率参数反向输入η_h2和η_heat合计不能超过 0.9否则电解槽“凭空产生能量”模型会出现反常的循环出力。储热二进制导致的问题u_st让模型变成 MILP求解完成后查看对偶乘子要格外谨慎CBC 对 MILP 不保证提供有意义的影子价格。提示做参数扫描时最好把建模部分封装成run_case(carbon_price, p2h_max, quota_rate)函数返回总成本和各设备日发电量。这样改一组参数就是一次函数调用避免手工改代码导致约束漏改。5. 验证热电优化结果合理性的三个技巧5.1 零碳价退化测试把三个碳价全部改成 0模型应当退化为“不考虑碳成本的热电联供经济调度”。如果此时 CHP 出力反而比有碳价时更低说明碳约束写错了比如配额公式方向反了或者排放系数漏了一项。极值测试是成本最低的排错手段。# 伪代码示意封装后直接跑两个算例 base run_case(carbon_price[50, 70, 120]) no_carbon run_case(carbon_price[0, 0, 0]) print(base[chp_e], no_carbon[chp_e]) print(base[buy_e], no_carbon[buy_e])正常结果是no_carbon的 CHP 发电量不低于base购电量不高于base。如果违反这个方向优先检查total_emission里是否把buy_e的间接排放重复计算了。5.2 对偶乘子校验只对 LP 模型可靠想用电平衡约束的影子价格评估某时段扩容 P2H 或储热的价值需要先把模型里的 0-1 变量去掉即允许蓄热罐同时充放重新作为 LP 求解。这时用value(prob.constraints[elec_12].pi)可以拿到该时段电平衡的边际价格与峰谷电价对照校验是否合理。MILP 求解器给出的对偶信息没有统一含义不要直接引用。5.3 反事实对比定位 P2H 的价值来源P2H 的价值可能来自三处氢收益、弃风消纳、碳成本转移。要区分这三个来源最简单的方法是把氢价设为 0 再求解一次看 P2H 是否还运行。如果氢价为 0 时 P2H 基本不工作说明模型里 P2H 是“为制氢而制氢”不是真正在消纳弃风再把 P2H 容量上限改成 0对比总成本和碳缺口就能算出 P2H 在降低碳交易支出上的边际贡献。把这两个反事实算例的差值一起放进结果表基本就能定位 P2H 在当前场景中的主价值来自氢收益、弃风消纳还是碳成本转移——这三个来源对应的调度出力曲线差异很明显。本文还有配套的精品资源点击获取