这两年综合能源系统的调度优化,基本上已经到了“不算碳就不好意思打招呼”的阶段。我最近在复现一个园区级热电联供系统的优化算例,核心就落在标题里这四件事上:MATLAB、阶梯式碳交易、电制氢、热电优化。这四个词拆开看都不新鲜,但把它们捏在一个模型里,求解难度和工程细节立刻就不一样了。
这套东西解决什么问题?一句话就是:在一个同时供应电和热的能源系统里,加入阶梯式碳交易机制,再把电制氢这一条“电→氢→热/电”的柔性环节引进去,用MATLAB搭一个混合整数线性规划模型,算出未来24小时每个设备的出力计划,让系统总运行成本最低,同时碳排放成本可控、新能源消纳更好。
这个方向很适合两类人:一类是正在做综合能源系统优化相关毕业设计或课题的学生,另一类是刚接触园区能源调度、想知道碳约束到底怎么落到数学模型里的工程师。我会把从设备建模、阶梯碳价线性化、YALMIP代码到求解器调参、踩坑排查的完整经验都撸一遍,不整虚的,直接给能跑的东西。
1. 系统边界和主目标:先想清楚再动手建模
很多初学者拿到这类题目,第一反应就是堆设备:燃气轮机、余热锅炉、电锅炉、电解槽、氢燃料电池、光伏、风电、储能、蓄热罐,全塞进一个模型里。结果二进制变量爆炸,求解器跑十几个小时还不收敛,最后只能把算例删到只剩两个设备,失去研究意义。
这类项目的正确打开方式,是先画清楚系统边界,再明确优化目标,最后才是列方程。
1.1 综合能源系统热电优化的“最小可行框架”
我习惯把系统简化成能源枢纽(Energy Hub)视角:左侧是输入,中间是转换设备,右侧是负荷。输入侧就三类——外网购电、天然气、风光可再生能源;转换设备按用途分三层:
- 供电层:燃气轮机、氢燃料电池、光伏和风电
- 供热层:余热锅炉(回收燃气轮机排烟余热)、电锅炉
- 耦合层:电解槽(电→氢)、储氢罐、氢燃料电池(氢→电和热)、蓄热罐
负荷侧就两类:电负荷和热负荷。注意这里“热负荷”是广义的,既包括供暖、工艺蒸汽,也包括生活热水。把它们拆开后,整个系统变成了一个非常清晰的“电—热—氢”三网耦合图。
选择这个框架的原因很朴素:它抓住了综合能源系统最核心的两个耦合特性。第一,热电联产机组强迫性地把电出力和热出力绑在一起,改变电出力就会改变热出力;第二,电制氢环节引入了时间转移能力——多余的风光电变成氢存下来,晚上再通过燃料电池发回电或供热,相当于一个化学形态的储能,比电化学储能的能量密度高得多。
1.2 为什么碳交易偏偏选“阶梯式”
传统碳交易在优化模型里通常被写成线性碳价:每吨碳排放成本固定,比如单价300元/吨,目标函数里加一项(碳价 \times 排放量)就完事。线性碳价方便,但不符合实际政策导向——现实中碳配额交易价格一般随缺口增大而上升,排得越多,买配额越贵,这就是阶梯式(累进式)碳价的含义。
在数学模型上,阶梯式碳价天然是一个分段线性函数。价格递增使得这个分段函数是凸的,而凸的分段线性函数在MILP框架内处理非常优雅:不需要引入非线性的乘积项,只需要增加几个连续变量和二进制变量,把“超排配额量”分配到不同的价格区间里。
更关键的是,阶梯式碳价会给优化问题注入一个非常强的信号:当碳排放缺口可能跨入更高价格档位时,模型会自发觉醒——减少燃气轮机出力、启动电解槽储氢、或增加外购绿电,因为阶梯碳价放大了减排的边际收益。这就是为什么这类型题目近几年在学术论文里特别火,伦理基础是“碳约束市场化”,数学基础是“分段线性化”,两者一结合,模型和故事都成立了。
1.3 电制氢环节的真实作用
“电制氢”(Power-to-Hydrogen,P2H)听上去高大上,本质就是一个电解水的过程:电力输入电解槽,把水分解成氢气和氧气。在综合能源系统里,它承担的角色远不只是“产氢”,更像一个灵活的电力负荷和能量时移器。
举个实际例子,白天光伏大发或夜间风电过剩时,电价低甚至出现弃风弃光,电解槽加大功率把多余的电能转化成氢气存进储氢罐。等到电价高的峰时,氢燃料电池再把氢气转化为电力和热力供给负荷。一来二去,新能源消纳率上去了,系统总运行成本反而下降了。
不过P2H在模型里有个容易踩的坑:效率。电解槽的电-氢转换效率一般只有60%到75%(按氢气低热值计),氢燃料电池的氢-电转换效率大约40%到55%。两个环节串起来,从电到电的往返效率约30%到40%。这意味着氢循环不可能是零成本的“白嫖”,它一定是通过“低电价储氢、高电价放电”的价差来赚钱的。模型的参数设置必须忠实反映这个效率损失,否则你会得出“全天开着电制氢就能省钱”的荒谬结论。
2. 设备模型逐个击破:燃气轮机、电解槽、储能一个都别含糊
建模是整个项目的地基。我见过太多人花80%的时间在调求解器上,结果模型本身有个地方少了爬坡约束、换热效率写反,怎么调都是白费功夫。这里把每个关键设备的数学模型拆开讲,尽量说清楚“为什么这么写”。
2.1 燃气轮机+余热锅炉:电热耦合怎么表达
燃气轮机是系统的发电主力,也是碳排放的主要来源。它的模型至少有四个约束:出力上下限、爬坡限制、电量与热量的耦合关系、燃料成本线性化。
出力上下限比较好写:某台燃气轮机的电出力 (P_{gt}) 必须在最小稳定技术出力 (P_{gt,min}) 和额定容量 (P_{gt,max}) 之间:
[ P_{gt,min} \le P_{gt}(t) \le P_{gt,max} ]
爬坡约束反映的是机组响应速度,保证上一时段到下一时段的出力变化不超过某个阈值 (\Delta P_{gt}):
[ -\Delta P_{gt} \le P_{gt}(t) - P_{gt}(t-1) \le \Delta P_{gt} ]
热电耦合是这类系统的灵魂。如果是背压式供热机组,热出力与电出力基本成正比,可以写成一个简单的线性关系 (H_{gt}(t) = k_{th} \cdot P_{gt}(t))。如果是抽凝式机组,可行域是一个多边形,常见表达是:
[ H_{gt}(t) \le \alpha_1 P_{gt}(t) + \beta_1, \quad H_{gt}(t) \ge \alpha_2 P_{gt}(t) + \beta_2 ]
也就是说,热出力上限随着电出力变化——电出力大,排烟余热多,能抽出的热也多;但电出力低到一定程度,热出力最小受稳定燃烧限制,不能无限低。
余热锅炉的模型更简单,它把燃气轮机排烟中的余热拾取出来变成供热。设余热回收效率为(\eta_{hrb}),则:
[ H_{gt,out}(t) = \eta_{hrb} \cdot H_{gt}(t) ]
这里我强烈建议在代码里把燃气轮机排烟热作为一个单独的中间变量存下来,不要直接写成 (H_{gt,out} = \eta_{hrb} \cdot k_{th} \cdot P_{gt})。中间变量虽然多了一行代码,但调试时能清楚看到“排烟热是多少、回收热是多少、散失热是多少”,排查热不平衡问题会轻松很多。
天然气燃料成本的计算也要注意。如果按立方米气价折算到发电功率上,公式是:
[ C_{fuel}(t) = \frac{c_{gas}}{\eta_{gt} \cdot LHV_{gas}} \cdot P_{gt}(t) ]
其中 (c_{gas}) 是天然气到厂价(元/Nm³),(LHV_{gas}) 是天然气低热值(kWh/Nm³),(\eta_{gt}) 是燃气轮机发电效率。习惯了这两个单位换算的坑之后,你就会明白为什么很多论文里干脆把燃料成本系数当成一个常数 (\lambda_{fuel})(单位:元/MWh)来写,因为这样可以省掉不少单位换算烦恼。
2.2 电解槽和氢能回电:两个效率参数搞定模型
电解槽在系统里是一个纯电负荷,输入电功率 (P_{el}(t)),输出氢气的等效热值功率 (P_{h2,out}(t))。为了避免在“kg氢气”和“MWh”之间来回换算出错,我推荐直接用等效功率表达:
[ P_{h2,out}(t) = \eta_{el} \cdot P_{el}(t) ]
举个例子,电解槽效率 (\eta_{el}=0.65),输入10 MW电,产出的氢气若全部燃烧放热相当于6.5 MW的热功率。这样储氢罐的SOC(荷氢状态)更新公式就很自然了。设储氢量 (S_{h2}) 以MWh为计量单位,充放氢效率分别为 (\eta_{ch2}) 和 (\eta_{dis,fc})(放电时按燃料电池效率折算):
[ S_{h2}(t+1) = S_{h2}(t) + \left( \eta_{ch} \cdot P_{h2,out}(t) - \frac{P_{fc}(t)}{\eta_{fc}} \right) \cdot \Delta t ]
注意最后一项是反着除的:燃料电池发电功率 (P_{fc}) 是输出,要用氢气量去填补这个输出,分母上要除一个效率。
这里有一个非常隐蔽的问题:符号方向。如果 (\Delta t) 的单位是小时,功率单位是MW,那么能量的单位就是MWh,没问题。但如果你哪天把功率单位设成了kW,而时间单位还是小时,能量单位就成了kWh,和储氢罐容量单位不一致,求解器的数值尺度立刻就会病态。我建议整个项目统一用MW和MWh,这能省下大量排查时间。
氢燃料电池放热的部分也需要建模。燃料电池发电和放热可以近似写成联产关系:(H_{fc}(t) = \mu_{fc} \cdot P_{fc}(t)),其中 (\mu_{fc}) 是热电比,一般在0.8到1.1之间。这意味着氢燃料电池在夜间顶峰放电时通常还能贡献可观的热出力,在“热电双缺”的时段特别值钱。
2.3 电锅炉、储能与热负荷侧的约束设计
电锅炉是弥补热出力不足的简单设备:热出力等于电输入乘以电热效率 (\eta_{eb}),典型效率在0.95到0.98区间。它的模型只有上下限和效率,没什么花头。但要注意,电锅炉的存在会让“电-热”耦合进一步强化:供电紧张时,电锅炉反而成了负担;供热量紧张时,它又能用便宜电补热。因此它的调度结果往往和电价、热负荷峰值密切相关。
蓄热罐的模型和蓄电池几乎一样,只是存储的对象从电能变成热能。设蓄热罐的储存量为 (S_{tes}(t)),充放热功率分别 (H_{tes,ch}(t))、(H_{tes,dis}(t)),考虑散热损失系数 (\eta_{tes,loss}):
[ S_{tes}(t+1) = \eta_{tes,loss} \cdot S_{tes}(t) + \left( H_{tes,ch}(t) \cdot \eta_{tes,ch} - \frac{H_{tes,dis}(t)}{\eta_{tes,dis}} \right) \cdot \Delta t ]
蓄热罐的加入让热侧有了时间转移能力,可以在电价低谷时段用电锅炉蓄热,在电价高峰时段放热,从而削峰填谷。注意蓄热罐的充放不能同时进行,通常需要引入二进制变量,但这会显著增加模型复杂度。如果只是日前调度且时间不长,可以用 (H_{tes,ch}(t) + H_{tes,dis}(t) \le H_{tes,max}) 这样线性化的方式,避免同时充放的不合理结果,同时不用二进制变量。
3. 阶梯式碳交易建模:分段线性化写进MILP的两种做法
碳交易建模是整个项目里最“学术”的部分,也是查阅论文时最容易被绕晕的一环。我先从底层逻辑讲明白。
3.1 碳配额核算:基准法怎么落到代码里
碳排放交易机制里,系统每时段有一个免费碳配额额度。这个额度的计算方式常用基准法:根据系统的发电量乘以一个碳排放基准值确定,比如:
[ E_{quota}(t) = \alpha_{quota} \cdot P_{gt}(t) ]
(\alpha_{quota}) 的单位是 tCO₂/MWh,取值一般在0.4到0.8之间,具体看你参考哪个行业基准。如果系统还有其他化石能源设备,配额也要相应叠加。
系统实际碳排放主要来自燃气轮机的天然气燃烧。简化计算可以用排放系数乘以燃料消耗量:
[ E_{actual}(t) = \beta_{CO2} \cdot \frac{P_{gt}(t)}{\eta_{gt}} ]
(\beta_{CO2}) 是燃气碳排放因子,单位 tCO₂/MWh(按燃气低热值计算),常用值约0.2 tCO₂/MWh。当然,这个值在不同文献里有差异,原因主要是天然气的组分不同。代码里写成一个可调参数就对了。
配额缺口就是:
[ \Delta E(t) = \max(0, E_{actual}(t) - E_{quota}(t)) ]
3.2 把阶梯碳价写进MILP的两种做法
最经典的方法是“分段变量法”。假设阶梯碳价设置了三个区间,缺口量分别为 (x_1, x_2, x_3),各区间长度上限 (L_1, L_2, L_3),对应碳价 (\lambda_1 < \lambda_2 < \lambda_3)。那么总碳成本就是:
[ C_{carbon} = \lambda_1 x_1 + \lambda_2 x_2 + \lambda_3 x_3 ]
关键是这几个 (x_i) 不能乱拆。比如缺口量是600吨,必须让 (x_1=500),(x_2=100),而不能让(x_1=100),(x_2=500)——因为后一种分配虽然总缺口一样,但碳成本算得不对(实际上要按区间顺序累加)。
要让模型自动按顺序填满区间,一种做法是引入二进制变量 (z_1, z_2),其中 (z_1=1) 表示“已经进入第二区间”,(z_2=1) 表示“已经进入第三区间”。约束写出来就是:
[ 0 \le x_1 \le L_1, \quad 0 \le x_2 \le L_2 \cdot z_1, \quad 0 \le x_3 \le L_3 \cdot z_2 ]
以及顺序约束:
[ z_2 \le z_1, \quad x_1 \ge L_1 \cdot (z_1 - z_2), \quad x_2 \ge L_2 \cdot z_2 ]
第一组约束保证了只有进入某区间,该区间的变量才能取正值;第二组顺序约束保证了阶梯必须从低到高依次填满。这样即使是价格递增的凸分段函数,也能用线性约束严格表达,不会出现区间倒挂。
另一种做法是用SOS2(Special Ordered Set of Type 2)变量,直接把分段线性函数包进去。不过SOS2在一些求解器里处理效率不高,而且很多人不熟悉,我建议先用二进制变量法,稳定可靠,调试门槛低。
3.3 目标函数组装:五项成本一起放进求解器
到这里,系统的目标函数可以总结为五项成本的总和最小化:
- 外网购电成本(购电功率×分时电价)
- 购气成本(燃气轮机燃料消耗×气价)
- 碳排放成本(阶梯碳价带来的那部分)
- 设备运行维护成本(各设备出力×单位运维成本系数)
- 弃风弃光惩罚(如果新能源出力没有全部消纳,给一个惩罚项)
这五项全部是决策变量的线性表达式,加上上一层列出的所有线性约束,最终模型就是标准的混合整数线性规划(MILP)。这也是为什么我不推荐在这个框架里用非线性求解器的原因——所有核心逻辑都可以线性化表达,用商用MIP求解器不仅能保证全局最优解,求解速度还快得多。
4. MATLAB实操:YALMIP建模主流程与三组场景对照
理论部分说完了,现在上实操。我用MATLAB加YALMIP工具箱搭模型。如果你用的是R2020a之后版本,直接命令行install_YALMIP就能装。求解器我习惯用Gurobi,速度最快;没有License也可以用开源的SCIP或者MATLAB自带的intlinprog,算例规模不大时它们也能跑。
4.1 典型日数据准备与单位统一
任何优化模型都要先有输入数据。一个标准算例需要准备:
- 未来24小时电负荷曲线(MW)
- 未来24小时热负荷曲线(MW)
- 未来24小时光伏出力预测曲线(MW)
- 未来24小时风电出力预测曲线(MW)
- 分时电价(元/MWh)
- 天然气到厂价(元/Nm³)
- 各设备的容量、效率、爬坡率、运维成本
我建议把所有输入数据写成一个结构体数组,比如params,最后统一传入优化函数。这样改参数只需要改一处。
单位方面,强烈建议全部都写清楚。如果电负荷是MW,电价是元/MWh,那么购电成本就是元,目标函数数值在(10^4)量级,碳成本在(10^3)量级,数值尺度会比较均匀。如果混用了kW和MWh,答案会出现(10^3)的数量级差,求解器在数值容差范围内很可能会优化出离谱结果。
4.2 YALMIP建模核心代码与求解器配置
下面给一个简化但可运行的主干代码,完整版建议按设备种类拆成多个build_*.m函数,方便排查。
%% 决策变量定义 P_gt = sdpvar(1, 24); % 燃气轮机电出力 H_gt = sdpvar(1, 24); % 燃气轮机排烟热(回收前) P_el = sdpvar(1, 24); % 电解槽输入电功率 P_fc = sdpvar(1, 24); % 氢燃料电池输出电功率 P_eb = sdpvar(1, 24); % 电锅炉输入电功率 P_buy = sdpvar(1, 24); % 从外网购电功率 P_sell = sdpvar(1, 24); % 向外部售电功率 S_h2 = sdpvar(1, 25); % 储氢量(MWh) C_seq = sdpvar(3, 24); % 碳排放缺口的三阶梯分解量 z_seq = binvar(2, 24); % 阶梯区间选择二进制变量 %% 约束组装 Constraints = []; % 电功率平衡:发电+购电 = 负荷+电解槽+电锅炉+售电 Constraints = [Constraints, P_gt + P_pv + P_wind + P_fc + P_buy == ... P_load + P_el + P_eb + P_sell]; % 热功率平衡:排烟余热回收 + 电锅炉供热 = 热负荷 Constraints = [Constraints, eta_hrb * H_gt + eta_eb * P_eb == H_load]; % 燃气轮机上下限与爬坡 Constraints = [Constraints, P_gt_min <= P_gt <= P_gt_max]; Constraints = [Constraints, -ramp_gt <= diff(P_gt) <= ramp_gt]; % 电解槽与储氢约束 Constraints = [Constraints, 0 <= P_el <= P_el_max]; Constraints = [Constraints, 0 <= P_fc <= P_fc_max]; Constraints = [Constraints, S_h2(2:end) == S_h2(1:end-1) + ... (eta_el * P_el - P_fc / eta_fc)]; Constraints = [Constraints, 0 <= S_h2 <= S_h2_max]; Constraints = [Constraints, S_h2(1) == S_h2_init, S_h2(25) == S_h2_init]; % 日始日末持平 % 阶梯碳交易约束 E_actual = beta_co2 / eta_gt * P_gt; % 实际碳排放 E_quota = alpha_quota * P_gt; % 配额碳排 Delta_E = E_actual - E_quota; % 缺口量 Constraints = [Constraints, Delta_E == sum(C_seq, 1)]; Constraints = [Constraints, 0 <= C_seq(1,:) <= L1]; Constraints = [Constraints, 0 <= C_seq(2,:) <= L2 .* z_seq(1,:)]; Constraints = [Constraints, 0 <= C_seq(3,:) <= L3 .* z_seq(2,:)]; Constraints = [Constraints, z_seq(2,:) <= z_seq(1,:)]; % 顺序填充约束(保证先填满低档) Constraints = [Constraints, C_seq(1,:) >= L1 .* (z_seq(1,:) - z_seq(2,:))]; Constraints = [Constraints, C_seq(2,:) >= L2 .* z_seq(2,:)]; %% 目标函数 carbon_cost = lambda1 * sum(C_seq(1,:)) + lambda2 * sum(C_seq(2,:)) + lambda3 * sum(C_seq(3,:)); fuel_cost = (c_gas / (eta_gt * LHV_gas)) * sum(P_gt); buy_cost = sum(c_buy .* P_buy); sell_income = sum(c_sell .* P_sell); om_cost = c_om_gt * sum(P_gt) + c_om_el * sum(P_el) + c_om_fc * sum(P_fc); Objective = buy_cost - sell_income + fuel_cost + carbon_cost + om_cost; %% 求解 ops = sdpsettings('solver', 'gurobi', 'gurobi.MIPGap', 0.01); ops.verbose = 1; optimize(Constraints, Objective, ops);上面这段代码不是完整的,但给出了建模的骨架。注意几个细节:
diff(P_gt)得到23个差值,爬坡约束默认需要首尾时段也要管,可以另外用P_gt(1) - P_gt_init处理第一时段。- 储氢约束里 (S_{h2}(25) == S_{h2}(1)) 是很多论文里都会加的“日始日末持平”约束,目的是让调度结果可以周期复现,避免第一天净放空储氢、最后一天净充满这样的作弊方案。
- 阶梯碳约束里的 (L1, L2, L3) 在代码里要确保是标量或与
C_seq维度匹配的向量;如果每时段单独计算配额,(L) 可以设成足够大的常数,但最好先算一下实际峰值,否则大M取值不当会引发数值问题。
4.3 三组场景对照实验的设计与结果解读
我建议不要只跑一个算例就交差,做三组对照,最能体现阶梯碳价和电制氢各自的价值:
- Case 1:无碳交易。目标函数只包含购电、购气、运维成本。
- Case 2:引入阶梯式碳交易,但不配置电解槽和燃料电池。
- Case 3:阶梯式碳交易 + 电制氢 + 储氢罐 + 氢燃料电池。
跑完三组,对比总成本和碳排放量,通常能看到很有意思的结论。Case 2 对比 Case 1,总成本可能上升(因为碳成本计入),但碳排放显著下降;Case 3 对比 Case 2,总成本可能不升反降,因为电解槽在低谷电价时段制氢、高峰时段燃料电池放电,套利收益抵消了一部分设备运维成本,同时进一步压低了碳排放。
我也遇到过Case 3反而更贵的情况,原因是电解槽和燃料电池的运维成本系数设得偏高、氢燃料电池效率设得偏低,经济性被抹掉了。这时候不用怀疑模型错了,调低运维成本系数或提高效率,再观察结果。这正是对照实验的作用:它能帮你暴露参数设置是否合理。
4.4 结果后处理:用图把“热电碳氢”四条线说清楚
优化结果出来之后,光看数值总成本远远不够。我一般至少画四张图:
- 电功率平衡堆叠图:燃气轮机、光伏、风电、购电、燃料电池、电解槽、电锅炉各自的正负贡献。
- 热功率平衡堆叠图:余热回收、电锅炉供热、蓄热罐充放热、热负荷曲线。
- 储氢罐SOC曲线:看储氢量一天内的波动规律。
- 碳排放缺口与阶梯区间分布图:看每个时段的碳排放缺口落在哪个价格档。
画图用MATLAB的area函数最合适。比如:
figure; h = area(t, [P_gt; P_pv; P_wind; P_fc; P_buy]'); % 各发电源堆叠 legend({'GT','PV','Wind','FC','Grid'});这些图能直观看出模型的调度逻辑是否合理。比如正常情况下,电解槽应该在凌晨谷电时段开启,燃料电池在傍晚高峰时段放电;如果结果出现白天电解槽满发、晚上燃料电池闲置,那大概率是电价曲线或设备效率设反了。
5. 踩坑记录:从不可行解到数值病态的五类排查经验
做完前面几步,模型能不能一次跑通,全看这里了。我总结了一下自己复现这类模型时反复踩过的坑,列成五类。
5.1 不可行解:先查平衡方程和上限值
YALMIP报错Infeasible problem是最常见的开头。90%的情况是功率平衡约束写反了,比如把+P_fc写成了-P_fc;还有10%是某台设备出力上限设置得太小,在负荷峰值时段无论如何都凑不够电。
排查技巧:把约束拆开,逐个用assign和check验证。YALMIP的check(Constraints)函数能告诉你哪些约束残差很大。或者干脆先跑一个不包含碳交易约束和储氢约束的版本,如果它可行,再一步一步把约束加回去,找到出问题的那一步。
5.2 数值尺度:p.u.归一化比你想的重要
当你的系统容量是几十MW,而碳配额缺口只有几十吨,目标函数里会出现(10^6)和(10^2)的跨度,这会让求解器陷入数值病态。解决办法就是做标幺化(p.u.):所有功率量除以基准功率 (P_{base}=100) MW,成本量除以基准成本 (10^6) 元,这样模型里所有数值都在0到10的范围内。
很多成熟的课题组直接就把所有模型写成p.u.格式,目的就是省掉后续调求解器数值容差的麻烦。
5.3 二进制变量不收敛:大M别贪大
阶梯碳交易里的大M值如果设成1e6,二进制变量配合这个量级,很容易导致MIP求解器花费大量时间在割平面和分支定界上,而且可能因为数值误差导致错误剪枝。
经验法则是:大M只比理论最大值大5%。比如某时段碳排放缺口最大不会超过50吨,那么(L2)写成55、(L3)写成60就够了。开始写代码时可以保守一点,但跑通后要记得回头收紧。
5.4 阶梯碳价分段验证:先把解固定再画曲线
如果你怀疑碳交易这一段线性化有问题,有一个很有效的验证手段:先把其他决策变量固定为某组已知可行值,只留碳排放缺口相关的变量,然后极小化总碳成本,看模型是否自动按顺序填满阶梯。如果它把缺口直接分配到高价段而没有填低价段,说明顺序约束写错了。
这种“固定变量后单独验证一段逻辑”的调试方法,对MILP模型特别有效,比盯着一大堆约束看半天强得多。
5.5 求解器参数:MIPGap和求解时长怎么权衡
商用求解器默认的MIP Gap通常是1e-4,对24时段、几十个二进制变量的模型来说其实有点严。实际算下来,最优解和次优解的差距可能只有0.1%,但却要多花几倍时间。
我一般把MIP Gap设置成1e-2(1%),对“运行成本优化”这类工程问题,1%的次优性完全可接受。如果你非要证明自己是全局最优,可以调到1e-3,但要有心理准备,某些时段可能会跑几个钟头。
另外求解器启动参数也值得调。Gurobi的MIPFocus如果设成1(最优性优先)或2(可行性优先),在不可行解和次优解之间平衡时会有明显差异。先用默认参数跑一轮看结果,再根据情况调整。
做完整套流程,我的体会是:这类型项目的难点从来不是某个设备的模型有多复杂,而是不同环节之间的耦合太容易出各种隐蔽错误。一个碳配额基准值的改动、一个大M取值的变化、一个功率平衡符号的翻转,都能让结果面目全非。所以每一次跑通,都建议把所有参数的算例结果存档,形成自己的参数基线。等以后做敏感性分析、多场景对比时,基线数据能帮你快速定位到底是参数变了还是代码改坏了。
如果你正要复现类似系统,我还有一个非常具体的建议:先别看整篇论文的复杂模型,先搭一个只有燃气轮机、电网购电、电负荷、热负荷、简单碳价的“最小模型”,跑通之后再加电制氢,再加阶梯碳价,加蓄热罐。每一步都在原有结果上做差分验证,别一上来就摊大饼。这样下来,大概率一天之内就能把整套模型扛下来。