氢能这个词,这两年但凡做能源方向的人多少都绕不开。综合能源系统里加氢能,不是赶时髦,是电、热、气这几种能量形式之外,确实需要一个能跨时长、跨介质存储的角色。我在做某园区综合能源优化调度模拟项目时就踩过不少坑,这套Matlab代码实现方案算是我反复调试后比较顺手的版本,今天把它拆开聊聊,从建模思路到代码落地,再到那些不跑一遍根本发现不了的细节问题,希望能给正在啃这块的同学一点参考。
1. 氢能在综合能源系统里的角色定位与核心价值
1.1 为什么优化调度里必须考虑氢能
做综合能源系统优化调度,首要问题是搞清楚系统里到底有哪些“能量流”。常规的园区级系统,典型配置是光伏、风电、燃气轮机、电锅炉、储能电池,外加从电网购电、从天然气管网购气。这类系统的调重点是:电负荷跟生产生活直接相关,热负荷有季节性特征,而气负荷相对刚性,来源也比较单一。
如果把氢能加进去,整个系统的能量通道就多了一条:可再生能源出力富余时,电解槽把多余电力转化为氢气,储氢罐存起来;需要时,燃料电池再把氢气转化为电和热。加上氢负荷可以直接供给氢燃料车、工业用氢,这就不再是单纯的“电力多了存电池”的问题了。
具体来说,氢能在优化调度里的价值集中在三个层面:
- 时间维度:电储能通常按小时级充放设计,氢储能则可以实现日级甚至周级的能量搬移,天然具有“跨日调峰”的属性,对连续阴天、长时间低风速场景的应对能力明显更强。
- 能量形态转换维度:电-氢-热、电-氢-气,氢是少有的能连接电、热、气三种能量网络的介质,这个属性让它成为多能互补系统的天然枢纽。
- 经济性维度:把“弃电”转化为可售卖的氢,本身就是一种收益渠道。在分时电价和阶梯气价的机制下,氢能设备的运行策略直接影响整体购能成本和碳排指标。
1.2 优化调度问题的本质:满足供需平衡的前提下让总成本最低
优化调度,无论多复杂,本质就一句话:在满足系统各项供需平衡和设备运行约束的前提下,决定每一时刻每台设备出多少力,让总成本最低。
这句话拆开有四个关键词:
- 时段:调度周期通常是24小时,最小单位是1小时,更精细的研究有15分钟分辨率的,但计算量会明显上升;
- 决策变量:各设备出力、储能充放功率、购能功率、氢流量等,其中启停状态、氢气买卖状态属于0-1整数变量;
- 约束:能量平衡、设备出力上下限、爬坡约束、储能SOC上下限与充放不能同时进行、氢罐储氢量上下限、联络线功率限制等;
- 目标函数:综合运行成本最小化,包括购电成本、购气成本、设备运行维护成本、启停成本,必要时还要加碳交易成本或弃可再生能源惩罚项。
在我做的模拟项目中,加入氢能模块后,优化问题的复杂度提升了一个明显的台阶。电和热原本可以直接对应各自的平衡方程,氢能一进来,电平衡里多了电解槽和燃料电池两个可变负荷/电源,热平衡里多了燃料电池余热,气平衡里多了电解槽原料气和燃料电池补给气,三个平衡方程全部联动,模型规模和求解难度同步上升。
2. 优化调度模型的数学构建:目标函数与约束条件
2.1 目标函数设计:成本项怎么列才合理
综合运行成本我通常按下面这个式子组织:
总成本 = 购电成本 + 购气成本 + 设备运行维护成本 + 启停成本 + 碳交易成本(可选) - 售氢收益(可选)
购电成本用分时电价乘以从电网购电功率,购气成本用天然气价格乘以购气流量,这两项基本是固定套路。设备运行维护成本需要注意单位系数的选择:燃气轮机、电锅炉这类设备的运维成本按输出功率折算,电解槽和燃料电池则按实际工作功率乘以运维系数,这是需要根据设备手册或参考论文标定参数的细节,很多初学者在这块直接拍脑袋给值,后续结果容易失真。
启停成本是加入0-1变量的直接原因。厂商的燃气轮机和电解槽每次冷态启动都有燃料损耗和寿命折损,如果不写进目标函数,求解器会让机组频繁启停来应对负荷波动,结果看起来“省了电费”,实际运行完全无法操作。
售氢收益建议单独列为一项,且要设一个约束条件:氢气产量不能超过本地下游氢需求的预测上限,否则模型会把电解槽开到最大,把一切多余电量都变成氢气,表面上看总成本很低,但氢根本卖不掉,纯属自欺欺人。
2.2 关键设备建模:电解槽、燃料电池、储氢罐怎么简化
设备建模是整个研究里最花心思的部分,简化过度会被审稿人质疑,简化不足又没法算。我最终采用的是一组兼顾精度与计算可行性的线性化模型。
电解槽(P2G):输入是电功率,输出是氢气流量,两者近似成线性关系,即:
F_H2_in(t) = η_p2g × P_el(t)
其中η_p2g是电解槽电转氢效率,范围一般在0.55~0.75之间。需要特别注意的是电解槽存在最小运行功率,低于这个功率时效率很低且设备损耗大,所以模型里加了运行下限约束,并且在目标函数里加了一个很小的启动惩罚系数,避免频繁开关。
燃料电池(FC):输入氢气,输出电和热,效率和热电比两个参数决定电气输出比例:
P_fc(t) = η_fc × F_H2_out(t) × LHV_H2 H_fc(t) = R_ht × P_fc(t)
LHV_H2是氢气低位热值,大约33.3 kWh/kg,R_ht是热电比,常见范围在0.8~1.2。燃料电池同样有最小负载率要求,我一般设置为额定功率的10%~20%。
储氢罐:建模逻辑和电池储能基本一致,核心是氢量的状态转移方程:
S_h2(t+1) = S_h2(t) + η_ch × F_H2_in(t) - F_H2_out(t) / η_dis
相比电池储能,氢罐的储放效率差异更大,电解槽和燃料电池自身的转换效率已经叠在设备模型里,罐体的充放损失相对较小,所以有的研究直接把储氢罐当作无损耗的库存来建模。我建议保留一个0.95~0.98的充放效率系数,这样更接近实际运行,运算代价也不大。
2.3 约束条件里最容易踩的坑
能量平衡约束看似简单,实际是错误的高发区。电力平衡必须覆盖园区电负荷、电解槽功耗、燃料电池出力这三方,缺一项平衡方程就会推导出荒谬结果。热平衡中燃料电池的余热要按热电比计入热源侧,却有很多代码里只算了燃气轮机的热出力,导致热源侧缺了一块。
充电和放电不能同时进行这类互斥约束,在电池储能里几乎人人都会加,但电解槽和燃料电池的“同时工作”逻辑却容易遗漏。从能量守恒角度看,两者同一时刻都工作是可能的(一边从电网买电制氢,一边燃料电池发电),但实际工程中这种状态几乎不会出现,因为它意味着能量被白白损耗掉。我习惯加一组互斥约束,让电转氢和氢转电在同一个时段内互斥,加完这组约束后,求解结果的总成本通常会降低,因为求解器不再把能量浪费在无意义的循环上了。
3. Matlab代码实现:模型落地的完整流程
3.1 代码框架与数据结构设计
Matlab代码实现,我用的是YALMIP建模框架。相比直接用求解器原生的编程接口,YALMIP的语法更接近数学公式,写优化模型几乎是一条约束一行代码,排查问题也方便。求解器用Gurobi或CPLEX都可以,免费环境下至少也要用SCS或者ECOS做一个可行性验证。
整个代码框架我拆成六个模块:
- 主脚本(main.m):定义算例场景,装配各模块,调用求解器并汇总结果
- 数据参数文件(data.m):电价、气价、负荷曲线、设备参数全部集中在这里,后续做敏感性分析时最常用
- 变量定义部分:定义所有决策变量,含连续变量和0-1整数变量
- 约束组装部分:按照设备类型和系统平衡关系逐条添加约束
- 目标函数部分:组装综合成本表达式
- 结果后处理脚本:绘制各设备出力曲线、储氢罐SOC曲线、能量平衡校验表
变量定义部分的核心代码大致长这样:
%% 时段数 T = 24; %% 连续变量 P_el = sdpvar(1, T); % 电网购电功率 P_gt = sdpvar(1, T); % 燃气轮机发电功率 P_elz = sdpvar(1, T); % 电解槽输入电功率 P_fc = sdpvar(1, T); % 燃料电池输出电功率 H_elz = sdpvar(1, T); % 电解槽产氢流量 S_h2 = sdpvar(1, T); % 储氢罐储氢量 %% 整数变量 u_gt = binvar(1, T); % 燃气轮机启停状态 u_elz = binvar(1, T); % 电解槽启停状态 u_fc = binvar(1, T); % 燃料电池启停状态3.2 约束组装的关键写法和容易出错的语法细节
电力平衡约束在YALMIP中非常直觉:
Constraints = [Constraints, P_el(t) + P_gt(t) + P_fc(t) == ... P_load(t) + P_elz(t) + P_eb(t)];这段代码的意思是大电网购电加上燃气轮机和燃料电池的出力,要等于电负荷、电解槽耗电、电锅炉耗电之和。如果设备数量多,建议把这些表达式先写成向量,再组装成约束,来做预填充。
储氢罐的SOC状态转移约束用循环逐时段写:
for t = 1:T Constraints = [Constraints, S_h2(t+1) == S_h2(t) * loss + ... H_elz(t) * eta_ch - H_fc(t) / eta_dis]; end这里面最容易出错的细节是时段的边界问题。S_h2(1)对应初始氢量,S_h2(T+1)对应调度结束后的氢量,通常要求与初始氢量相等,否则模型会把罐里的氢气全部排空来省成本——因为罐内氢气理论上已经花过成本了,不排白不排,但这在实际运行中是不可接受的。我实测过不设终端约束的情况,储氢罐会固定被清空,燃料电池会在最后几个时段疯狂发力。加上末态约束后,这个病就好解决了。
互斥约束有几种写法:
Constraints = [Constraints, u_elz(t) + u_fc(t) <= 1]; Constraints = [Constraints, P_elz(t) <= P_elz_max * u_elz(t)]; Constraints = [Constraints, P_fc(t) <= P_fc_max * u_fc(t)];这种写法比用大M法更稳定,因为YALMIP自带处理逻辑,也不容易出现数值病态。
3.3 求解配置与参数联动
求解调用部分,我通常写:
ops = sdpsettings('solver', 'gurobi', 'verbose', 2, 'showprogress', 1); % 如果机器内存有限,还可以限制求解时间 ops.gurobi.TimeLimit = 120; ops.gurobi.MIPGap = 0.01;MIPGap也就是1%的差距,对调度研究来说精度已经足够,能让求解速度提升不少。纯MILP问题用默认精度也跑得动,但如果系统规模上来了(例如节点数超过几十个),建议显式设置MIPGap,否则模型会为了证明那0.1%的差距多跑十几分钟。
求解完成后做结果检查:
optimize(Constraints, Objective, ops); %% 设备出力曲线和后处理 P_gt_opt = value(P_gt); P_elz_opt = value(P_elz); S_h2_opt = value(S_h2);需要特别注意:value函数返回的已经是求解后的数值,但在使用前要先检查problem字段是否为零。如果problem不等于0,先不要急着画图,大概率是约束不满足或求解超时,画出来的结果也没什么意义。
4. 求解器选型与参数调优经验
4.1 YALMIP + Gurobi 的配置流程及常见报错
开发环境是Matlab + YALMIP + Gurobi。第一次配置时最常遇到的问题就是把Gurobi的安装路径、许可证和Matlab的路径设置搞乱。安装时要特别注意一个点:Gurobi安装包里的matlab文件夹,要让Matlab加到搜索路径里,同时确保Gurobi求解器在YALMIP中被正常识别。
yalmiptest这条命令会把你系统里所有能用的求解器都列出来。我见过有人做完了整个模型,结果发现YALMIP里根本没有Gurobi,全靠内置求解器硬扛,效率差到无法忍受。
如果yalmiptest没有检测到Gurobi,最常见的处理方法是手动把Gurobi的matlab文件夹添加到Matlab路径中,然后重新启动Matlab,再运行yalmiptest确认。
4.2 模型规模爆炸时的降维方法
加氢能后模型规模会变大,这是必然结果,但有一些行之有效的降维方法:
- 聚合设备:同类型、同参数的设备可以聚合为一台等效设备,比如两台500kW电解槽聚合成一台1000kW的,约束数量减半,结果几乎无差异;
- 削减非关键时段分辨率:如果研究重点是调度策略,而不是设备启停时序细节,可以把调度周期从24小时改为12个2小时时段或8个3小时时段,规模成倍下降;
- 使用外逼近或基于场景的简化:对于多场景研究,可以用K-means聚类把365天聚类成3~5个典型场景来代表全年工况,这既降低了计算量,又保留了日类型之间的差异。
我在项目中就先把多台小功率电解槽聚合成了一台大电解槽,整版约束数量直接少了将近三分之一。聚合的前提是设备拓扑结构和参数一致,否则要谨慎使用。
4.3 目标函数单位一致性这个看不见的坑
电网电价通常是元/kWh,天然气价格是元/m³,设备运维成本容易顺手写成元/kW,三者单位不统一的话,目标函数里就会有人在数值上被“吃掉”。这类错误不会导致求解失败,但会让优化结果严重偏差,而且很难被发现。
我的统一做法是:
- 所有功率统一为kW,所有氢量统一为kg,所有能量统一为kWh;
- 天然气的热值按10 kWh/m³折算,氢气按33.3 kWh/kg折算;
- 所有成本统一折算到元。
比如燃气轮机的发电效率是η_gt,天然气的能量转换关系写为:
P_gt(t) = η_gt × F_ng(t) × LHV_ng
其中F_ng是天然气流量,LHV_ng是天然气低位热值。这样一气呵成,所有变量单位都是一致的,调试时逻辑清晰。
5. 实际调试中的问题排查与经验速查
5.1 四种典型异常及排查思路
模型求解后结果反直觉,或者求解器直接报错,这类问题有规律可循。我整理了一份速查表:
| 现象 | 常见原因 | 排查顺序 |
|---|---|---|
| problem返回1,无法收敛 | 约束存在冲突,尤其是互斥约束和上下限约束重叠 | 先把三分之一的约束注释掉,逐步恢复,定位到首次报错的那组约束 |
| 求解很快但结果全为0 | 目标函数单位错误或系数过小,被其他项掩盖 | 单独打印目标函数各项数值,看哪个成本在主导 |
| 储氢罐曲线异常尖峰 | 缺少充放互斥约束或末态约束 | 检查整数变量是否绑了上下限,检查SOC边界 |
| 购电量出现负数 | 忽略了与电网的功率流向限制 | 加双向功率约束,或明确购买和售电两个独立变量 |
5.2 结果合理性校验的必备清单
代码跑通不等于研究完成。我习惯在出结果后做一张“能量平衡校验表”,逐项对照每个时段的输入输出:
| 时段 | 电力供给合计 | 电力需求合计 | 偏差 |
|---|---|---|---|
| 1 | 购电+燃气轮机发电+燃料电池发电 | 电负荷+电解槽耗电+电锅炉耗电 | ≤0.01 kW |
| 2 | …… | …… | …… |
偏差大于0.01 kW基本说明约束或参数有误。热平衡用同样的方式逐时段核对,燃气轮机余热、燃料电池余热、电锅炉产热的和应当刚好覆盖热负荷。
这个工序一次都不该省。我见过有同学跑出来总成本异常低,高兴了半天,结果发现是某个平衡约束漏了一个设备项,整体方案根本不可行。
5.3 参数敏感性分析怎么做才高效
研究性文章几乎都要求敏感性分析,这一步操作起来很省事,但如果数据组织不好,会跑死人。
我的做法是写一个批量脚本,把光伏装机、氢价、电解槽效率作为三个变量,各自取三个水平(低、基准、高),形成27组算例,循环调用求解函数,最后汇总结果表格。这样一次循环下来,各种运行策略和成本的变化趋势全都清楚了。
批量跑分时需要注意的一点是每次循环都要重新初始化随机状态或数据源,否则可能出现不同算例用了相同负荷曲线而不自知。另一个细节是将结果保存为.mat文件而不是每次重新求解,这样后期画图修改格式时不需要重新跑模型。
6. 几个值得反复琢磨的进阶扩展方向
基础模型跑通之后,不妨往这几个方向延伸,能让课题的深度和价值提升一个层次:
多时间尺度协同:仅做日前24小时调度还不够,实际运行中光伏预测偏差很大会导致日内计划失效。可以考虑用“日前优化 + 日内滚动修正”两层结构,日内层以小时或分钟为步长滚动优化,偏差通过调整电解槽功率和储能出力来消除。
低碳运行与碳流追踪:双碳目标下,可以把碳交易成本或者碳捕获设备纳入模型,探讨氢能在降低系统碳排放中的作用。此时目标函数需要增加碳排放项,约束条件也要增加碳配额或碳价格的表达,整体模型会从纯MILP变成带碳约束的扩展MILP,难度适中但很有价值。
氢能与电转热多能耦合加深:电解槽本身运行时会有废热产出,这部分余热如果加以回收用于供热,可以进一步提高能源利用效率。模型里需要把电解槽热回收项加上,同时考虑热平衡约束的变化。这个方向实操性强,而且是实际工程中很关心的问题。
含氢交通负荷的互动:如果园区里有氢燃料公交或物流车,那么氢负荷就不再是固定的,而是和车辆调度相关,可以做成联动优化,探讨交通与能源两个系统之间的互动关系。这一点在逻辑上和组织架构上都有很大的讨论空间。
这些方向未必都要做全,但至少选一个点深入下去,整个课题的立足点会明显扎实很多。
最后再分享一个我在实际调试中觉得很实用的习惯:不要把模型代码和数据处理代码混在一个文件里。每个环节单独成脚本,出错时定位到某个模块,会节省大量翻代码的时间。我会先把参数文件固化,再锁定约束模块,最后主文件只保留装配和调用逻辑。这样一套结构下来,无论参数怎么改、场景怎么换,主流程几乎不需要大改,跑实验的效率会明显提升。综合能源优化调度本身就是一个系统工程,模型、数据和程序三者之间的衔接往往比单个环节里的算法设计更决定最终结果的可靠性。