前段时间有个做综合能源规划的同行拿着风电和光伏的出力曲线来找我,说园区想加储能,准备直接上锂电池。我让他先别急着下单——如果调度目标只是把24小时内的用电成本压到最低,锂电池足够;可一旦把视角拉长到周、月甚至跨季节,氢和氨就绕不开了。这也是我最近一直在做的“含氢气氨气综合能源系统优化调度研究”最核心的出发点。这篇文章把整个项目的建模思路、约束推导、Matlab代码骨架以及调试过程中踩过的坑完整梳理一遍,算是一份能直接照着改的实操笔记。
这套系统到底能做什么,适合谁看?如果你在做微电网、综合能源系统、电氢耦合或者储能容量配置方向的研究,又不想只停留在概念图层面,想用Matlab跑出一套真正可复现的调度结果,那这篇内容应该能帮你节省不少时间。我会先讲清楚为什么氢氨联动值得做,再给出从物理拓扑到数学模型的完整转换过程,最后落到代码实现和算例分析上。全程不绕弯子。
1. 氢氨联动:为什么“电—氢—氨”会进入调度模型
1.1 传统储能在长周期调度面前显得吃力
先聊一个基础问题:为什么不能只靠电池?
电池储能的特点是响应速度快、循环效率高(一般在85%到95%),但它的问题是能量密度和自放电率。按现在的技术水平,锂电池电站的持续放电时长设计大多在2到4小时,少数项目做到6到8小时。应对日内的削峰填谷完全够用,可一旦碰上连续阴雨、连续静稳天气,风光的出力缺口会持续两三天甚至更久,电池就得把容量放大到几乎不经济的程度。
氢储能不一样。电解槽把富余电力变成氢气,储存在高压气态罐或者地下盐穴里,需要的时候再通过燃料电池发电。这个链路虽然“电到氢再到电”的综合效率不高,大概只有30%到40%,但储能的时长上限远高于电池。氢气罐本质上就是一个可以按天、按周、按季节存放的能量仓库,时间维度一拉长,它就体现出独特的战略价值。
1.2 氨不是终点,是高效的储氢载体
储氢解决了“长时间”问题,但氢气本身的储运又成了新麻烦。
标准状况下氢气的密度只有0.0899千克每立方米,别说运输,就是想在园区里找个位置放个大储罐,都要为容积发愁。高压气态储氢做到35兆帕或者70兆帕,体积密度还是有限;液态储氢需要把温度压到零下253摄氏度,液化过程本身就要消耗掉氢气能量的30%以上,代价很高。
氨在这里出现的意义就体现出来了。氨分子式是NH3,含氢量按质量计算约17.6%,在常压下零下33摄氏度就能液化,而且液氨的储运基础设施在化工行业里已经跑了很多年,非常成熟。更关键的是,氨的体积能量密度比液氢还要高,在相同体积下能携带更多有效氢。
所以“电—氢—氨”这条链路,本质上是用合成氨把“大规模储氢”问题转化成“大规模储氨”问题。氢不适合长途运输,氨可以;氢不适合在常温常压下长期保存,氨可以。氨在这个系统里就是一个高效的储氢载体。
1.3 调度层面多出的自由度
从优化调度的角度来看,加入氢和氨之后最直接的变化是多出了好几个可控环节。
电能不能直接变成氨,中间要经过电解槽制氢、氢气进入合成氨装置这两步;反过来,氨要变成电也不是一步到位,需要先经过氨分解或重整制氢,再进入燃料电池。每一个转换环节都对应若干台设备,每台设备都有启停状态、出力上下限、爬坡速率和运行成本。这等于在传统“电源—负荷”的刚性平衡之间,塞进了一组非常灵活的时间转移模块。
举个简单例子:凌晨风电大发、电价很低甚至为负,系统完全可以加大电解槽功率,把多余的电变成氢,再把氢合成氨存起来;等到傍晚用电高峰,氨分解制氢,氢气进燃料电池发电。这一进一出,相当于把凌晨的廉价电搬到了傍晚再用。调度模型要做的,就是决定每一台设备在每个时段该以什么功率运行,才能让整体成本最低。
从这个意义上说,含氢气氨气的综合能源系统优化调度,不再是简单的机组组合和负荷分配问题,而是一个多能互补、多时间尺度耦合的复杂优化问题。
2. 系统拓扑与建模准备:先画能流,再谈数学
2.1 系统的物理构成和能流关系
拿到一个实际项目,我不会急着建数学模型,而是先把系统的物理拓扑画清楚。这张图画明白了,后面所有约束和变量就都有了出处。
一套典型的含氢气氨气综合能源系统,大致可以分成四个板块:
- 供能侧:风电、光伏,加上外部电网作为备用电源。
- 转换侧:电解槽(电转氢)、储氢罐、合成氨装置(氢转氨)、氨储罐、氨分解装置(氨转氢)、燃料电池(氢转电)。
- 负荷侧:常规电负荷、氢负荷(比如加氢站需求)、氨负荷(比如化工用户或船用燃料需求)。
- 辅助系统:考虑热电联产时还有余热回收模块,不过第一版模型可以先忽略。
能流方向是:风电和光伏优先供给电负荷,富余电力进电解槽制氢;氢气一部分直接供给氢负荷,一部分进入储氢罐,还有一部分送到合成氨装置变成氨;需要发电的时候,储氢罐里的氢气或者储氨罐里的氨通过重整/分解重新变成氢气,进入燃料电池。整个系统在电力、氢气、氨三种能量载体之间来回切换,耦合关系非常强。
2.2 为什么最终收敛成一套混合整数线性规划
把这套系统写成优化模型,第一步要决定采用什么数学结构。我最开始想用纯线性规划,结果发现行不通,原因在设备启停。
电解槽、合成氨装置、燃料电池这些设备都有一个特点:不能像水龙头一样随意从0功率直接跳到额定功率。实际工程里,电解槽一般有最小负载率限制,比如额定功率的20%,低于这个值设备无法稳定运行;同时启动和停机还有额外成本。这些要求天然需要0到1的整数变量来描述“开/关”状态,模型自然就变成了混合整数线性规划,简称MILP。
MILP的好处在于:只要目标函数和约束都是线性的,求解器就能给出全局最优解,而不是像非线性规划那样可能陷入局部最优。对于调度这类需要可信结论的问题,全局最优性非常重要,这也是我在项目里坚持用MILP框架的根本原因。
2.3 时间尺度和分段粒度怎么定
调度模型的时间粒度也很重要,直接影响求解规模。
我第一版用的是全天24个时段、每时段1小时。这个粒度足够展示氢氨系统的日运行规律,比如凌晨制氢、傍晚发电,模型规模也不大,适合反复调参验证。如果要做周调度或者季度调度,可以保持1小时粒度扩展到168时段或者更长时间,代价是求解时间指数级上升。这时候通常会把电解槽、储氢罐这类慢动态设备单独拿出来做长时间尺度的聚合建模,把燃料电池等快响应设备留在短时间尺度里。两步走,才能在精度和求解效率之间找到平衡。
3. 目标函数与约束体系:调度模型的数学骨架
3.1 目标函数:成本最小化如何拆项
优化调度要回答的核心问题是:满足负荷需求的前提下,怎么样运行最省钱。所以目标函数我定义为系统总运行成本最小化。
总运行成本包括五个部分:
一是购电成本,系统从外部电网买电的费用,按分时电价计算。二是设备运维成本,电解槽、燃料电池、合成氨装置、氨分解装置运行都会产生损耗,按运维系数乘以运行功率估算。三是启停成本,设备每启动或停机一次都有人工检查和损耗成本,用0到1变量的差值来描述。四是弃风弃光惩罚,为了让模型尽量消纳可再生能源,在目标函数里加一个较大的惩罚系数,让弃风光在成本上“不划算”。五是碳排放成本或环境惩罚,这部分在碳交易试点地区可以直接用碳价折算。
对应到Matlab代码里,目标函数可以写成类似下面的样子:
objective = sum(sum(Price_electricity .* P_buy)) ... % 购电成本 + sum(C_om_electrolyzer * P_el + C_om_fc * P_fc) ... % 运维成本 + sum(C_start_el * z_start_el + C_stop_el * z_stop_el) ... % 启停成本 + sum(rho_curtail * (P_curtail_wind + P_curtail_pv)) ... % 弃风弃光惩罚 + sum(C_co2 * E_co2); % 碳排放成本实际项目中可以根据研究重点增减项。如果只关心经济性,前三项就够;如果做高比例可再生能源消纳分析,弃风弃光惩罚必须加,否则模型会“故意”弃掉多余的风光来省运维成本。
3.2 电氢氨全链路约束
约束条件是整个模型的灵魂,也是最容易出错的地方。我把约束分成四组来说。
第一组是电力平衡约束。任意一个时段,风力发电功率、光伏发电功率、外部购电功率、燃料电池发电功率之和,要等于电负荷、电解槽耗电功率之和,再加上弃风电量和弃光电量。弃风弃光在约束里作为松弛变量存在,保证模型在极端场景下仍然有可行解。
第二组是氢气平衡约束。电解槽产氢量加上储氢罐放出量,要满足氢负荷、燃料电池耗氢量以及合成氨装置用氢量之和。储氢罐的充放逻辑不是简单的同时进出,而是通过储氢罐的荷电状态来衔接前后时段,相当于一个时间耦合约束。
第三组是氨平衡约束。合成氨装置的产氨量加上氨储罐的库存变化,要满足氨负荷和氨分解装置用氨量。氨储罐的模型和储氢罐类似,也有容量上下限和当日初始/末尾状态要固定的约束,从而保持日调度的周期性。
第四组是外部电网交互约束。购电功率不能超过联络线容量,一般来说我不允许系统向电网反送电,因为涉及上网电价和交易机制,第一版模型没必要引入这个复杂度。如果想做双向互动,只需要加一个售电变量,结构上不复杂。
3.3 非线性效率曲线分段线性化
设备效率不是常数,电解槽在20%负载率和80%负载率的时候效率能差好几个百分点。如果直接用非线性效率曲线,模型会变成混合整数非线性规划,求解难度大幅上升。
我的做法是分段线性化。具体来说,把设备的输入功率范围切成几段,每段用一个线性效率表示。只要分段的段数足够多(一般取3到5段),拟合精度对最终调度结果的影响可以控制在很小范围内,而模型结构依然保持MILP。
举一个电解槽的例子。假设额定功率是5兆瓦,我将0到5兆瓦按最小负载率切成三个区间:0到1兆瓦区间设备不运行,1到3兆瓦区间效率较低,3到5兆瓦区间效率较高。引入分段变量后,每段的功率不能独立取值,需要满足相邻段的次序约束——不能跳过第二段直接用第三段。这个逻辑在YALMIP里可以用implies和binvar实现,但要注意避免引入过多二进制变量导致求解变慢。
4. Matlab建模与求解:YALMIP环境下搭一套可运行调度程序
4.1 工具链怎么选
我用的是Matlab加YALMIP工具箱再加一个MILP求解器。YALMIP不是求解器,它是一个建模层,让用户用比较接近数学公式的语法来描述优化问题,然后把问题转成求解器能懂的标准格式。求解器我推荐三个选择,根据授权情况自行安排:
| 求解器 | 类型 | 适合场景 | 备注 |
|---|---|---|---|
| Gurobi | 商业 | 大规模MILP,求解速度快 | 学术免费,性能强 |
| CPLEX | 商业 | 大规模MILP,兼容性好 | 老牌求解器,稳健 |
| CBC | 开源 | 中小规模MILP | 免费,速度一般 |
| MATLAB内置intlinprog | 内置 | 中小规模MILP | 无需额外安装 |
如果只是做课程设计或者入门练习,直接用Matlab自带的intlinprog就够;如果论文算例规模大,或者要跑周调度,建议用Gurobi,速度差别非常明显。
ops = sdpsettings('solver', 'gurobi', 'verbose', 2); ops.gurobi.MIPGap = 0.001; % 设置1%以内的gap,平衡速度与精度4.2 从参数定义到决策变量
建模之前先把数据整理成固定的数据结构。我的习惯是用MATLAB的struct,把负荷、可再生能源出力、设备参数、能源价格分开存放。这样后面对照修改某个参数非常方便。
决策变量分两类。连续变量包括各时段的购电功率、风电光伏实际出力、电解槽功率、燃料电池功率、储氢罐充放功率、氨储罐充放功率、弃风弃光量等。二进制变量包括每台设备的启停状态、启动动作和停机动作。
启动动作和停机动作是从启停状态衍生出来的。如果电解槽在t时段运行、t-1时段不运行,那么t时段有一个启动事件;反过来就是停机事件。这个关系需要一组线性约束来保证,最常用的写法是:
% z_start >= u(t) - u(t-1) % z_stop >= u(t-1) - u(t) % z_start >= 0, z_stop >= 0YALMIP里可以直接用diff配合循环来写。注意第一时段没有前一时段状态,要单独给定初始启停状态,否则约束会越界。
4.3 约束装配与求解的核心代码
下面给出一段可以运行的主干代码框架,它不是一个完整的工程代码,但包含了从变量定义到求解的全部核心步骤。实际使用时把参数表替换成自己的数据即可。
%% 参数准备 T = 24; % 时段数,单位小时 P_wind_forecast = ...; % 风电预测出力,1xT P_pv_forecast = ...; % 光伏预测出力,1xT P_load = ...; % 电负荷,1xT H_load = ...; % 氢负荷,1xT N_load = ...; % 氨负荷,1xT %% 定义决策变量 P_buy = sdpvar(1, T); % 购电功率 P_el = sdpvar(1, T); % 电解槽输入功率 P_fc = sdpvar(1, T); % 燃料电池输出功率 P_wind_use = sdpvar(1, T); % 风电实际使用功率 P_pv_use = sdpvar(1, T); % 光伏实际使用功率 P_curtail = sdpvar(1, T); % 弃风弃光总功率 u_el = binvar(1, T); % 电解槽启停状态,1为运行 u_fc = binvar(1, T); % 燃料电池启停状态 SOC_h2 = sdpvar(1, T+1); % 储氢罐荷电状态 SOC_nh3 = sdpvar(1, T+1); % 氨储罐荷电状态 %% 约束条件 C = []; % 电力平衡 C = [C, P_wind_use + P_pv_use + P_buy + P_fc == P_load + P_el + P_curtail]; C = [C, P_wind_use <= P_wind_forecast, P_pv_use <= P_pv_forecast]; C = [C, P_wind_use >= 0, P_pv_use >= 0, P_curtail >= 0]; % 电解槽上下限与启停关联 P_el_min = 0.2 * P_el_rated; % 最小负载率20% C = [C, P_el_min * u_el <= P_el <= P_el_rated * u_el]; % 燃料电池上下限与启停关联 C = [C, P_fc_min * u_fc <= P_fc <= P_fc_rated * u_fc]; % 储氢罐状态更新 C = [C, SOC_h2(2:T+1) == SOC_h2(1:T) + (eta_h2_in * H2_in - H2_out / eta_h2_out)]; C = [C, SOC_h2(1) == SOC_h2(T+1)]; % 日始日终一致 C = [C, SOC_h2_min <= SOC_h2 <= SOC_h2_max]; % 燃料电池耗氢量计算 H2_fc = P_fc * HeatRate_fc; % 热电关系,按实际参数转换 H2_to_NH3 = ...; % 进入合成氨的氢量 H2_in = eta_el * P_el; % 电解槽产氢量 H2_out = H2_fc + H_load; % 储氢罐放出量 %% 目标函数 objective = sum(Price_elec .* P_buy) + ...; %% 求解 ops = sdpsettings('solver', 'gurobi', 'verbose', 2); optimize(C, objective, ops); %% 结果提取 P_buy_opt = value(P_buy); P_el_opt = value(P_el); SOC_h2_opt = value(SOC_h2);这段代码只是骨架,实际运行前还要处理不少细节。比如储氢罐的动态方程里,H2_in、H2_out是衍生变量,需要写清楚它们和电解槽功率、燃料电池功率之间的关系;合成氨装置的用氢量要从氨负荷反推出来。这些量之间环环相扣,建议在写出全部约束后用YALMIP的check函数一行一行验证约束是否满足,而不是直接去看目标函数值。
4.4 结果输出与可视化
求解完成后,除了直接看目标函数值,我几乎必做三张图:第一张是电力平衡图,把风电、光伏、购电、燃料电池发电、电解槽耗电五个分量堆叠在同一张面积图里,一眼就能看出每个时段电从哪里来、到哪里去;第二张是储氢罐和氨储罐的SOC曲线,观察它是否满储满放;第三张是各设备启停状态的时间轴,用横向柱状图表示,方便发现一些不合理的高频启停。
figure; area([P_wind_use_opt; P_pv_use_opt; P_buy_opt; P_fc_opt]'); legend('风电', '光伏', '购电', '燃料电池'); xlabel('时段/h'); ylabel('功率/MW'); title('系统电力平衡图');如果发现某台设备出现频繁启停的“抖振”现象,一般不是求解器的问题,而是目标函数里缺少启停成本或者启停成本系数太小。找到原因后不要盲目调大惩罚系数,先分析目标函数各项的量级再调整,会更有效。
5. 算例验证:一天24小时的调度结果分析
5.1 测试数据与设备参数
下面这道算例是基于一个假设园区设计的,数据经过脱敏处理,但量级和趋势都贴近实际。风电和光伏都装了一定容量,某典型日的出力曲线呈现“夜间风电大、白天光伏大”的组合特征;电负荷呈现白天双峰特点,晚高峰出现在19点前后。
设备参数如下:
| 设备 | 额定功率 | 最小负载率 | 效率 | 运维成本系数 |
|---|---|---|---|---|
| 电解槽 | 8 MW | 20% | 75%(含电源变换) | 0.02 元/kWh |
| 储氢罐 | 4000 kg | — | 充放效率90% | 0.005 元/kg |
| 合成氨装置 | 5 t/h | 40% | 80% | 0.03 元/kg |
| 氨储罐 | 100 t | — | 充放效率95% | 0.01 元/kg |
| 燃料电池 | 3 MW | 30% | 50%(电效率) | 0.06 元/kWh |
| 氨分解装置 | 2 t/h | 30% | 85% | 0.04 元/kg |
电价采用分时电价:谷段(0点到7点、22点到24点)0.25元/kWh,平段(8点到10点、16点到17点)0.55元/kWh,峰段(11点到15点、18点到21点)0.95元/kWh。弃风弃光惩罚系数设为2元/kWh,确保模型优先消纳。
5.2 三种方案的对比结果
我做了三组对比:方案A是含完整氢氨链路的综合能源系统,方案B把氢氨系统去掉、只保留锂电池储能(容量按等效能量4兆瓦时),方案C完全不配置储能。
主要指标如下:
| 方案 | 总运行成本(万元) | 弃风弃光率 | 购电量(MWh) |
|---|---|---|---|
| A 氢氨综合系统 | 4.82 | 1.2% | 38.5 |
| B 锂电池储能 | 5.36 | 6.8% | 44.2 |
| C 无储能 | 6.10 | 19.5% | 52.3 |
氢氨综合系统的总运行成本比无储能状态低约21%,比锂电池方案低约10%。这组数据说明,在这个场景下氢氨系统不仅把弃风弃光率从19.5%压到1.2%,还因为多出了跨时段能量转移能力,明显减少了峰时段高价购电量。
5.3 从结果看氢氨系统的运行规律
再看优化求解给出的运行策略,规律性很强。
凌晨时段风电大发、电价很低,电解槽进入满负荷运行,产出的氢气一部分直接供给氢负荷,另一部分进入合成氨装置生成液氨存进氨储罐。这一阶段燃料电池是停机状态,因为此时发电不划算,不如把氢先存起来。
白天光伏出力上来后,电解槽的功率会根据实时电价和光伏预测动态调整,储氢罐在午间光伏大发时继续充入氢气。傍晚电价进入峰段,燃料电池启动,储氢罐开始放氢发电,缓解晚高峰的电负荷压力。如果储氢罐容量不够,氨储罐里的氨会经过氨分解装置重新变成氢气,再进燃料电池,相当于给系统加了第二层“跨时段保险”。
这个“凌晨制氢—白天储氢—傍晚放氢发电”的模式,是氢氨系统在日前调度结果里最典型的运行规律。你如果跑出来的结果不是这样,大概率是某个参数量纲或者约束方向写错了。
6. 实际开发中的坑与调参心得
6.1 电解槽最小负载率带来的“不可运行区间”
这是我第一次建模时踩得最深的坑。
电解槽的功率本来是一个连续变量,按上下限约束可以直接写0到额定功率。但实际上很多电解槽在额定功率的0到20%区间根本无法稳定运行,一开机就是“要么不开,要么至少20%负载”。如果模型允许功率落在0和20%之间的“灰色地带”,求解器很可能会给出一个看起来成本很低、实际物理上完全不可行的方案。
解决办法就是引入启停变量,把功率变成“0”和“20%到100%之间”的互斥状态:
C = [C, P_el_min * u_el <= P_el <= P_el_rated * u_el];这条约束在自然语言里的意思是:如果u_el = 0,P_el被夹逼为0;如果u_el = 1,P_el的可行范围是[最小负载率,额定功率]。所有具有类似特性的设备,比如合成氨装置和燃料电池,都要用这个思路处理。
6.2 求解器性能与MIP gap的取舍
氢氨系统相比传统储能的调度模型,最大的变化是二进制变量数量增多。每台设备在每个时段都有启停状态,24个时段、6台设备,就是144个二进制变量,这还只是一个“小算例”。
二进制变量一多,MILP的求解时间会明显上升。有一次我调大系统规模到168时段的周调度,Gurobi跑了将近40分钟还没出结果。后来在求解器里设置合理的最优间隙MIP gap,比如0.1%或者0.5%,问题求解时间从几十分钟降到两三分钟,而目标函数值和理论最优值之间只差不到0.1%。
这个思路在科研和工程里都很实用。调度问题本身就不需要追求绝对的千万分之一精度,一个0.1%的最优间隙在工程上完全可以接受。不同的问题是“花30分钟把成本从4.8200压到4.8195”,还是在3分钟内拿到4.8215的结果,后者往往更符合实际决策需要。
6.3 索引、量纲和数据预处理
Matlab的索引从1开始,但很多论文公式里的时段从t=1到t=T,储氢罐的SOC状态却要涉及t=0到t=T。这让代码里的下标很容易错位。我的经验是,在代码里把SOC定义成长度为T+1的向量,SOC(1)表示初始状态,SOC(k)对应第k-1个时段结束后的状态。这样约束写起来虽然长一点,但不容易出现时序错位。
量纲问题是另一个容易翻车的地方。电解槽功率单位是兆瓦,产氢量单位是千克每小时,储能罐容量单位是吨,把“单位”混在一起写约束,数据会差出好几个数量级。我的习惯是所有能量相关量统一用标准单位,氢和氨的质量流直接用千克,再在参数表里通过热值或当量系数把功率和流量关联起来。转换系数也建议显式写成变量,不要用魔法数字硬编码,后面调参时能少掉很多头发。
6.4 容易被忽视的热量平衡约束
大部分初学者建氢氨模型时只关注电、氢、氨三种能量载体,容易忽略热负荷。
电解槽工作时会产生热量,燃料电池更明显,合成氨反应本身是强放热反应,氨分解则是强吸热反应。如果园区有热负荷需求,这些废热完全可以回收利用,但如果系统没有热负荷,多余的热量如何处理就成了一个工程问题。模型层面最简单的处理是加一个“可抛弃余热”变量,并用一个约束限制换热设备的容量上限。如果研究重点不是热电联产,可以在一开始假设余热免费可弃,后面再扩展。
我自己却建议不要完全不管它。校园、园区型综合能源系统往往有热水或蒸汽需求,把“电热氢氨”四网放在一起建模,不仅更贴近工程实际,而且审稿人和工程评审会认为你把系统耦合关系考虑得更完整。扩展也不难,只需在目标函数里增加热负荷供应成本,再加一组热平衡约束即可。
6.5 从简化系统起步,逐步加设备
最后一条经验最实在:不要一上来就把系统建完整。
我建议先从“风机+光伏+电解槽+储氢罐+燃料电池+负荷”的简化系统起步,把代码跑通、把约束验证无误,再逐步加入合成氨装置、氨储罐、氨分解装置。每加一个设备,就把那一个设备的目标函数贡献和关键约束单独测试一遍。这样做的好处是,一旦结果出现异常,你能很快定位是新引入设备的问题还是原有模块被破坏。
这套程序我前后改了两周才稳定下来,最深的体会就是“给调度问题做减法比做加法重要”。优化模型不是越大越好,也不是约束越多越接近实际,而是要在可求解性和物理准确性之间找到平衡点。你先用简化版本把核心链路跑通,再一层一层把细节叠加上去,整个项目的推进过程就会顺畅很多。