最近在做一个综合能源系统的热电优化调度项目,核心是研究阶梯式碳交易机制与电制氢(Power-to-Hydrogen,P2H)设备如何协同影响系统运行。这个方向在目前的能源系统优化里确实算一个研究热点,但很多论文讲模型讲得比较理想化,落到Matlab具体实现的时候,会遇到不少实际问题。我把自己从模型建立到代码跑通的完整过程整理了一下,希望对做类似方向的同学有点帮助。
这套方案解决的核心问题是:当系统接入光伏和风电之后,如何通过合理的调度策略,在满足电、热负荷需求的前提下,综合考虑购能成本、设备运维成本和阶梯式碳交易成本,让全天运行总成本最低。适合电力系统、综合能源方向的研究生,以及从事园区级能源系统规划或调度的工程师参考。
1. 项目到底在做什么:从调度问题说起
1.1 为什么要在综合能源系统里引入阶梯式碳交易
先说背景。传统的综合能源系统优化,目标函数一般就是最小化购电购气成本和设备运行维护费用,碳约束往往用一个固定碳价折算到成本里。固定碳价的做法有个明显问题:它对碳排放的约束力度是线性的,排多排少边际成本都一样,系统的碳减排意愿完全取决于碳价定多高。碳价定低了,机组调度结果偏向经济性,碳排放量很难压下来;碳价定高了,又会过度牺牲经济性。
这时候阶梯式碳交易机制就体现出优势了。它的核心思想是:碳排放量超出免费配额越多,超出部分的碳价越高——类似阶梯电价,用分段的碳价把“多排多付”包装得更激进。这样低碳排放的机组在调度中会自动获得经济优势,系统会在经济效益激励下主动调整各设备的出力计划,而不是靠一个硬性约束去“卡”排放。
从调度模型角度来理解,阶梯式碳交易机制的引入实际上改变了目标函数的凸性构造。碳排放量作为决策变量经过分段价格函数映射后,原本线性目标被改造成了分段线性目标,而且分段点之间的决策变量是相互耦合的,必须通过0-1变量来建模,这直接决定了问题从线性规划(LP)升级为混合整数线性规划(MILP)。接下来我会重点展开这部分。
1.2 电制氢在这里扮演什么角色
电制氢设备接入综合能源系统,最直接的作用是给系统增加了“电转气”的灵活性。把低谷时段或弃风弃光时段的多余电能转化为氢能,储存起来,需要的时候再通过燃料电池发电或直接作为燃气锅炉、氢燃气轮机的燃料使用。这样一来,系统不再单纯依赖蓄电池来平抑风电光伏的波动,而是多了一条能量时移的路径。
从热电联产的角度来看,引入P2H设备最大的价值在于解耦“以热定电”的刚性约束。传统的热电联产机组为了保证供热,必须保持一定的最低电出力,导致夜间风电大发时弃风严重。有了电制氢设备之后,夜间多余的风电可以用于制氢,既消纳了可再生能源,又能在碳交易机制下通过替代部分化石燃料进一步降低系统碳排放量。
在实际工程中,P2H设备的效率模型一般是线性的:输入电功率与输出氢量之间近似成正比关系,但受设备额定容量约束。这个线性假设在项目初期足够精确,如果追求更细的建模精度,可以把效率曲线改成多项式甚至非凸的,但那样求解难度会急剧上升,这个后面说实现细节的时候再展开。
2. 核心模型拆解
2.1 阶梯式碳交易模型的数学化描述
阶梯式碳交易机制在模型层面需要处理这样几个要素:初始碳配额、实际碳排放量、碳交易区间划分、以及不同区间对应的碳交易价格。
碳配额的确定方式在文献里常见的有两种:一种是按照系统负荷大小乘以一个排放因子来分配免费配额,另一种是给定一个固定的配额总量。在这个项目里,我采用的是基准线法——根据系统的电负荷和热负荷,按照一个统一的基准排放强度来分配免费配额,公式大致是:
E_quota = δ_e * P_load + δ_h * H_load
其中δ_e和δ_h分别是电、热负荷的基准排放强度,P_load和H_load是某时段内的电、热负荷功率。这种分配方式的好处是贴合实际工程中“多出力多得配额”的原则,系统无法靠简单降低负荷来套取配额收益。
实际碳排放量E_actual则需要根据各机组的出力折算,主要包括燃煤热电联产机组、燃气锅炉燃烧排放,以及从电网购电对应的间接碳排放。具体到每台设备,碳排放量等于设备耗能量乘以燃料的碳排放因子,购电的部分就用购电量乘以电网平均排放因子来折算。
有了配额和实际排放,就可以计算碳交易量:
E_trade = E_actual - E_quota
当E_trade为正,说明系统碳排放超标,需要在碳市场购买配额;当E_trade为负,表示系统有盈余配额,可以在市场上出售获得收益。
阶梯式的差异主要体现在购买或出售配额的价格上。如果采用三阶梯结构,可以写成:
- 第一阶梯:购碳量不超过区间下限时,碳价为c1
- 第二阶梯:购碳量介于两个阈值之间时,超出部分碳价为c2 = c1 + σ
- 第三阶梯:购碳量超过第二个阈值时,再超出的部分碳价为c3 = c1 + 2σ
这种阶梯式表达对模型求解提出了新要求:目标函数里碳交易成本不再是线性的,而是分段的。在Matlab的Yalmip环境中,处理分段函数的常规思路是引入0-1整数变量和Big-M辅助变量,把每一段区间单独激活,并用大M约束保证同一时刻只能激活一个区间。我用了一个长度为T的时间序列数组来表示每小时的碳交易变量,然后引入对应的三段0-1变量,效果稳定。
2.2 电制氢与CHP等设备建模
设备建模是整个优化模型的地基,各设备的数学模型直接决定了约束条件的写法。
热电联产机组(CHP)的建模是这个系统里最关键的。我采用的可调度域模型,用两个变量分别表示电出力和热出力,然后通过一组线性不等式约束来限定其运行范围。典型的约束包括最小和最大电出力、对应的热出力上下限,以及电热出力之间的关系系数。这个可调度域的建模方式比简单的“热电比固定”的办法灵活很多,因为实际运行中CHP机组完全可以在一个范围内调节电热比。
燃气锅炉的模型相对简单,热出力乘以锅炉效率得到耗气量,需要考虑的是最大热出力和爬坡约束。电锅炉同样用一个线性模型,热出力等于耗电量乘以制热系数COP,调度时它可以快速调节功率,灵活性比燃气锅炉好,但运行成本受电价影响很大。
电制氢设备的关键参数是制氢效率和额定功率。制氢量等于输入电功率乘以制氢效率除以氢气的低位热值。在约束方面,除了功率的上下限,我还会加一个最小负载率的限制,避免设备运行在极低效率区间。一些文献会忽略这一点,但实际工程中电解槽在低负载下运行不仅效率差,还会影响设备寿命。
蓄电池储能采用带时间耦合的状态模型,电池的荷电状态(SOC)更新方程表示为上一时刻SOC加上充电功率乘以充电效率减去放电功率除以放电效率再除以电池容量。约束条件包括SOC上下限、充放电功率限值以及防止同时充放电的0-1变量约束——这个在MILP模型里很重要,否则求解器有可能会让电池同时充放电来套利。
储氢罐模型与蓄电池类似,用储氢量作为状态变量,更新方程是上一时刻储氢量加上制氢量减去供氢量,约束包括储氢容量上限和充放氢速率限制。
2.3 目标函数与约束条件
整个优化问题的目标函数是最小化一天内的总运行成本,包括四个部分:购电购气成本、设备运维成本、碳交易成本,以及为了鼓励消纳可再生能源而加入的弃风弃光惩罚成本。购电成本采用分时电价,购气成本按单位热值的天然气价格计算。设备运维成本用各设备出力乘以对应的运行维护单价。碳交易成本就是上面提到的分段函数,实际是碳交易量与对应阶梯碳价的乘积。
约束条件方面,最核心的是各时段的功率平衡约束。电功率平衡等式包含电网购电量、CHP电出力、风电出力、光伏出力之和,等于电负荷加上电锅炉耗电、P2H设备耗电和电池充电功率减去电池放电功率。热功率平衡等式包含CHP热出力、燃气锅炉热出力、电锅炉热出力之和,等于热负荷。氢平衡等式包含制氢设备的产氢量与储氢罐释放量之和,等于燃料电池消耗氢量与供氢负荷之和。
对偶性原理在这个系统里很直观:如果某个平衡约束被违反,优化一定会利用某个设备的调节能力去修复它,而成本最高的那台调节设备决定了最终平衡的边际成本。我在调试模型时会专门检查平衡约束,因为一旦出现停电或停热的“伪最优解”,八成是某个平衡方程漏写了某个变量。
设备约束除了出力上下限,还包括爬坡约束和运行状态约束。对于CHP机组,我设置了每小时的爬坡速率限制;对于储电和储氢设备,设置了充放状态互斥的0-1变量;另外还考虑了购电功率的上下限,防止模型在低电价时段无限购电。
3. Matlab实现全过程
3.1 整体代码框架与变量设计
整个Matlab程序我分成了四个模块:参数初始化模块、变量定义模块、约束构建模块和求解与结果输出模块。这种结构最大的好处是方便做参数敏感性分析——改碳价、改设备容量只需要动参数初始化模块,其他部分不用动。
参数初始化模块定义了所有基础数据,包括24小时的电负荷、热负荷、风电和光伏出力曲线、分时电价、天然气价格,以及各设备的效率、容量上下限、爬坡速率等。特别要说的是分时电价的区间划分,它直接决定了P2H设备的启停策略——如果谷电价格和平电价格差距不够大,制氢设备在谷时段运行的经济性就不明显,整个模型的结果会跟预期有较大出入。
变量定义模块使用Yalmip的sdpvar和binvar来定义决策变量。连续变量包括各设备的出力、购电购气量、碳排放量、储电储氢状态量等,0-1变量包括电池充放电状态、P2H启停状态、碳交易区间选择状态等。碳交易区间选择的状态变量是阶梯式碳交易模型能否正确求解的关键,我把每一小时要引入三组0-1变量,分别对应三个碳交易区间。
电网购电量的变量我单独定义了正变量,不采用自由变量,目的很明确——避免模型在求解过程中出现同时购电和售电的套利行为。在实际的电力市场规则下这种套利是不允许的,如果变量本身允许为负,模型会钻空子。
3.2 用Yalmip建模的关键写法
Yalmip建模最核心的一点是:尽量用向量化方式写约束,而不是写24个循环。一方面代码简洁易读,更重要的是约束数量太大时,循环写法会让Yalmip内部表达式的构建速度急剧下降。
电功率平衡约束的向量化写法大概是这样:
constraints = [constraints, P_grid + P_chp + P_wt + P_pv == ... P_load + P_eb + P_p2h + P_bat_dis - P_bat_chg];P2H设备的制氢量约束这样表达:
Q_h2 = eta_p2h * P_p2h / LHV_H2; constraints = [constraints, P_p2h >= 0.1 * P_p2h_max]; % 最小负载率约束 constraints = [constraints, P_p2h <= P_p2h_max];阶梯碳交易的分段线性化处理是整个模型最具挑战性的部分。我给每小时的碳交易量写了这样的判断逻辑:
% 三个0-1变量分别决定碳交易量落在哪个阶梯 constraints = [constraints, E_trade <= A1.*u1 + M*(1-u1)]; constraints = [constraints, E_trade >= lower_bound1.*u1 - M*(1-u1)]; % 第二、第三阶梯的约束同理,u1、u2、u3互斥等式约束M的取值要特别小心:M太小会把有效解空间误杀,M太大会导致求解时数值稳定性变差。我在实际调试中会把M设置成相应变量的上限值的10倍,通常能兼顾两方面的需求。
3.3 求解器选择与调试要点
模型构建完成后,求解器选择直接决定求解速度和稳定性。我用的是Yalmip调用Cplex,在部分对比测试中也用过Gurobi。两者的线性规划求解器差距不大,但在MILP问题上,Gurobi的节点处理效率在多数问题上略快一些。不过Cplex在约束数量非常多的情况下表现更稳健,所以我最终以Cplex作为主求解器。
求解命令很简单:
options = sdpsettings('solver','cplex','verbose',2); optimize(constraints, objective, options);但调试过程中真正让人头疼的不是求解器本身,而是模型构建的错误。我遇到过这么几种情况:约束写错导致可行域为空;目标函数里变量类型不匹配导致报错;SOC更新方程里状态变量索引错位,导致储能设备的容量在循环中被错误地放大或缩小。排查这些问题的通用方法是先把约束条件一条一条注释掉,看模型能不能恢复正常求解,然后逐段放回,定位问题约束。
另外一点经验是,模型测试时先用24小时且仅含一台CHP机组和一个P2H设备的简化系统跑通全流程,再扩展到大系统。直接用完整模型调试,一旦出错,变量太多根本无从查起。
4. 算例设置与结果分析
4.1 系统参数与场景设定
为了验证阶梯式碳交易机制和电制氢设备的效果,我设计了三组对比场景。
场景一:不含碳交易机制、不含P2H设备的基准系统,只有燃气锅炉和CHP供应热负荷。 场景二:在基准系统上引入固定碳价碳交易机制,但不含P2H设备。 场景三:完整系统,包含阶梯式碳交易机制和P2H设备。
这样做对比的意义在于:场景一和场景二对比可以得出碳交易机制对系统运行成本与碳排放量的影响;场景二和场景三对比可以得出P2H设备在碳交易环境下的“边际”价值。
系统参数方面,参考某个典型工业园区微网数据设定。CHP机组额定电功率100kW,热功率80kW;燃气锅炉额定热功率150kW;电锅炉50kW;P2H设备额定功率30kW,制氢效率按4.5kWh/Nm³氢折算;蓄电池容量50kWh;储氢罐容量10kg。风电装机100kW,光伏装机50kW。
碳交易参数的设置需要依据能源系统研究的一些常用范围来设定。免费碳配额按负荷基准排放法确定,基准排放因子电负荷取0.5kg/kWh,热负荷取0.06kg/kWh。固定碳价场景的碳价取40元/t,阶梯碳交易场景的基础碳价取25元/t,每超出一个阶梯碳价上调10元/t,整个阶梯划分为三个区间。
4.2 结果怎么对比,怎么证明方案有效
求解完成后,我会从三个维度来分析结果:运行成本与碳排放、各设备出力曲线、以及风电光伏消纳率。
运行成本方面,场景三的总成本比场景一降低了大约12%,比场景二降低了约6%。这个结果符合预期——引入碳交易后,系统有了减排激励,会主动调整运行策略;而阶梯碳价比固定碳价的惩罚机制更强,因此在控制碳排放着更有效的表现。
碳排放方面,场景三的碳排放在三个场景里最低,比场景一减少了21%。这里特别值得关注的是,减排量不仅来自碳交易机制的价格引导,还来自P2H设备的替代效应——谷电时段制氢替代了一部分燃气锅炉供热,降低了天然气消耗,碳排放自然就降下来了。
看看P2H设备在场景三里的运行曲线就会发现它的运行规律非常清晰:夜间谷电时段以满功率制氢,白天高峰期不运行,储氢量在夜间攀升,在白天供热峰时释放。这套运行策略实质上起到了“能量时间转移”的作用,和蓄热电锅炉在电价低谷蓄热是一个逻辑。不过P2H的转换链条更长——电→氢→热,中途效率损失较大,所以只有在碳价够高、谷电价够低的情况下才有经济性。
风电光伏消纳率的提升在场景三中最为显著,从场景一的不足80%提升到了接近95%。消纳率提升的直接原因是P2H设备和电锅炉能够在风电大发时段吸收多余的电力,让风电场不再因为系统调节能力不足而被迫限电。
为了让结果更有说服力,我还会做一组碳价敏感性分析,测试基础碳价从15元/t变化到55元/t时,总成本和碳排放的变化趋势。结果呈现出典型的“边际效益递减”特征:碳价从15元升到30元时碳排放下降明显,但从40元升到55元时碳排放变化趋缓,说明系统的减排潜力存在物理上限,一味靠提高碳价换减排是不经济的。
5. 常见问题和排错实录
5.1 求解不收敛或结果无法满足平衡约束
这是我在项目初期遇到最多的问题,也最让人头疼。原因通常有三个:Big-M值不合理、0-1变量过多导致求解困难、约束中存在冗余矛盾。
Big-M值方面,如果M设得太大,比如直接取1e6,求解器在分支定界过程中数值稳定性会变差,导致精度下降,甚至出现无解或溢出。如果M太小,某些本应可行的调度方案会被误杀。我在项目里对每一类需要Big-M的约束单独设置了M值,取值标准是该约束中最大可能变量值乘以一个系数(通常取10倍以上),并确保M严格大于所有可行解的实际取值。
0-1变量过多导致的求解慢问题,我的处理办法是缩减时间粒度。24小时模型每个时段都要引入多个0-1变量,如果设备种类也很多,0-1变量数量会迅速膨胀到几百个。这样会显著增加求解难度,导致求解时间变长。在做初步验证时,可以先把时间粒度改成每一时段代表2小时,一共12个时段,验证完模型逻辑后再改回1小时间隔。
约束冗余矛盾的排查方法,建议用“逐步注释法”。先只保留电功率平衡和所有设备上下限约束,求解一次;如果能解再逐段加入热平衡、氢平衡、碳交易约束;一旦加入某个约束后无解,问题一定出在那组约束里。这样定位问题非常有效。
另一个容易被忽略的地方是碳排放量的计算。CHP机组的碳排放量与它的发电和供热都有关系,如果只按电出力折算、忽略了供热部分的排放,碳交易量会明显偏低,模型结果看起来“成本很优”,实际却是算错了账。燃气锅炉的排放也要单独计入碳排放。
5.2 储氢罐状态异常或P2H设备一直停机
P2H设备在优化结果中如果始终不出力,通常要从经济性上找原因。制氢设备每生产1kg氢气,对应的电力成本,要低于它替代掉的天然气成本,设备才会启动。如果谷电价过高或天然气价格过低,经济账算不过来,设备自然闲置。这时候先检查一下谷时段的电价设定值,再看天然气单价的取值是否合理。
储氢罐SOC状态异常,比如全时段保持满容量或者空容量,一般是因为模型中缺少最小储氢量或储氢量目标值的约束。我在模型中加了一条末时段储氢量不小于初始储氢量的约束,避免优化结果把全部氢气在最后一个时段耗尽,这种“末端甩卖”策略在现实中不是稳态运行方式。
还有一种情况是储氢罐的充放速率约束没有设,导致模型在某个单时段让储氢罐瞬间充满或放空。从数学模型角度看这样没有违反任何约束,但物理上是不可能的。这类问题我统称为“模型物理一致性缺失”,在搭建约束时就要预判。
5.3 碳价参数对结果的影响有多大
碳价参数是这个模型中敏感度最高的参数之一,对优化结果的影响比对设备效率参数的影响大得多。原因是碳价直接改变了原本由经济性主导的机组出力排序,一旦碳交易成本高到足以覆盖天然气和煤的成本差异,会促使系统在调度时优先选择低碳机组,甚至改变整个系统一个小时级的出力模式。
建议做参数灵敏度分析时,不要只看单一指标变化,最好一口气输出总成本、碳排放量、弃风率、P2H设备利用率这几个关键指标的对比表格。这样能够更全面地反映出碳价变化的影响范围。比如碳价提高以后,总成本先上升,但当系统通过增加P2H出力弥补碳成本时,总成本又开始回落,这种非线性响应趋势只有通过多维指标对比才能清晰捕捉。
另外特别提醒一下,如果你的Matlab环境是R2023a以后的版本,Yalmip和Cplex/Gurobi的接口偶尔会出现兼容性问题。我踩过Cplex在R2023b下MILP求解器无法启动的坑,后来用Gurobi作为备用求解器解决的。正因为这个教训,后来在代码里加了求解器自动切换逻辑,一旦Cplex求解超时或报错,自动改调Gurobi重试。这个兜底策略在主模型调试阶段节省了大量等待时间。
6. 一些实操中的体会
这个项目从零开始到现在,我最大的感受是:数学模型写起来漂亮,落地实现才算数。阶梯式碳交易机制和电制氢设备单独拿出来都很好建模,但一旦组合进一个多能互补系统,变量之间的耦合关系会让调试工作量成倍增加。MuPad、Simulink、Yalmip这些工具箱我都试过,最终选择纯Yalmip + Cplex/Gurobi的方案,不是因为其他工具不好,而是这个组合在MILP问题上的处理效率和灵活性确实更胜一筹。
调试过程中帮了很大忙的一个习惯,是把每一小时的电功率平衡做成图形化检查。求解完成后画出各电源出力柱状堆叠图,和负荷曲线放在一起看,昼夜之间的供需匹配关系一目了然。很多时候模型报“求解成功”,但结果在工程上不合理,就是靠这种可视化检查发现的。
最后再分享一个小技巧:储能类设备的状态更新方程里,SOC的初始值不要设成0,否则前几个小时储能设备的运行会失真。我在模型里把蓄电池初始SOC设成0.5,储氢罐初始容量设为上限的一半,这样系统从仿真一开始就处于一个接近稳态的工况,比从零状态起步所得到的全天调度结果更接近实际运行情况。这个细节不处理好的话,前几个时段的优化结果会明显偏离工程直觉。