尧图网页设计 
做综合能源系统优化调度最怕的是什么不是模型复杂而是你辛辛苦苦建完模、跑完代码结果却不符合物理直觉——明明有大量弃风调度结果却在用电高峰期买高价电。我第一次在Matlab里跑含氢气氨气的综合能源系统优化调度时就栽过这个跟头后来才明白问题出在我们没有把氢和氨的时序平移特性真正建模进去。氢和氨跟电、热不一样电基本是即发即用热可以靠蓄热罐扛几个小时而氢和氨是真正能实现跨天、甚至跨周存储的化学储能介质。它们可以把风光大发时的多余电力转成化学能存起来等到负荷高峰再通过燃料电池或氨燃料机组发电相当于在物理层面给电网加了一个大容量充电宝。这套逻辑说起来简单真正落到Matlab代码里涉及设备建模、多能流平衡约束、混合整数变量处理、求解器调参等一系列问题。这篇博文就围绕含氢气氨气综合能源系统优化调度这个方向把建模思路、Matlab代码实现、调试心得和扩展方向一次性讲透。适合看这篇内容的读者正在做综合能源系统、氢能/氨能方向课题的研究生准备用Matlab做算例验证的工程师以及想了解氢氨耦合系统调度原理但对代码无从下手的初学者。我会尽量用能直接抄作业的方式写但也希望大家理解每一步背后的物理含义否则遇到求解不收敛、结果反直觉的时候排查起来会非常痛苦。1. 为什么要在综合能源系统里引入氢气氨气调度问题的本质变化1.1 氢与氨不是多了一个设备而是多了一条时间维度传统综合能源系统优化调度研究对象通常是电、热、气三种能量流。电力靠电网平衡热靠供热管网平衡天然气靠气网平衡。三者之间存在或松或紧的耦合燃气轮机同时产电产热热电联产电锅炉可以把电转成热燃气锅炉直接把气转成热。这种耦合是空间上的、同一时刻内的。氢气加入后问题立刻不一样了。电解槽把电转成氢氢气可以储存在高压储氢罐里也可以直接供给氢负荷还可以在燃料电池里重新转成电。这个电→氢→电的链路中间有存储环节所以它天然带时间属性我可以选择在风电最便宜的时候制氢在电价最高的晚高峰放氢发电。调度决策从这一时刻怎么分配能量变成了这一时刻是否应该存能量/放能量。氨作为氢的载体进一步延伸了这条链路。氨合成塔把氢和氮合成氨储氨罐比储氢罐更容易大规模存储氨裂解器又可以把氨重新分解成氢或者直接把氨送去掺烧发电。这样氢和氨就构成了一个可双向转换、可长期存储的能量走廊。从调度建模的角度看引入氢氨意味着模型中增加了储能状态变量储氢量、储氨量、反映转换效率的线性/非线性函数、以及一组跨时段耦合约束。1.2 氢氨耦合给调度带来的三个棘手问题第一个问题是多能流平衡约束的联立。电平衡、热平衡、氢平衡、氨平衡必须同时满足任何一个网络失衡其他网络都会间接受影响。比如电解槽开工不足氢产量下降不仅氢负荷受损氢燃料电池也发不出电最后可能连晚高峰的电平衡都被打破。第二个问题是储能时序耦合带来的计算复杂度上升。储氢罐和储氨罐的容量约束、进出能率约束、周期结束时的存储量约束都是跨时段约束。模型规模随调度周期线性增长24小时优化问题还比较轻松如果做成8760小时的年度生产模拟变量和约束数量会非常庞大对求解器和代码质量的要求明显提高。第三个问题是分段线性效率和启停逻辑让模型变成混合整数问题。电解槽在低负荷区间效率下降燃气轮机有最小技术出力储氢罐的充放状态可能涉及整数变量氨合成塔往往有连续运行时间约束。这导致整个优化模型通常是混合整数线性规划MILP或者混合整数非线性规划MINLP在Matlab里用Yalmip建模、再调Gurobi或CPLEX求解是最常见的做法也是我下面要详细展开的内容。2. 系统拓扑与设备建模先把能量流图变成数学关系2.1 一套典型的含氢氨综合能源系统长什么样在做Matlab实现之前必须先在纸上把系统拓扑画清楚。我在这类课题里最常用的一套结构如下供能侧风电机组、光伏阵列、上级电网购电、天然气网购气。转换侧电解槽电→氢、氢燃料电池氢→电热、燃气轮机天然气/氨/氢掺烧→电热、氨合成塔氢氮→氨、氨裂解器氨→氢、燃气锅炉气→热、电锅炉电→热。存储侧储氢罐、储氨罐、蓄电池、蓄热罐。负荷侧电负荷、热负荷、氢负荷、氨负荷比如工业用户的氨原料需求。从调度研究的角度每一个设备都是一个输入能量→输出能量的转换节点转换关系可以用效率常数或者分段线性效率曲线描述。这个系统明显比传统电热气三联供多了一层氢氨化学链但也正是这层化学链给了系统更大的调节灵活性——风光出力高了就多制氢制氨风光出力低了就放氢发电、氨裂解补氢。2.2 关键设备的稳态数学模型别一上来就写非线性对于24小时到一周的调度问题通常采用稳态模型忽略设备内部动态过程只关心单位时段内的能量转换关系。下面是我的常用建模方式风光出力直接取预测功率序列或者用场景法生成若干组典型场景。如果做鲁棒优化还会引入不确定区间。电解槽输入电功率输出氢能按热值或按质量。模型为P_H2_out η_elec * P_elec_in其中η_elec通常取0.6~0.75。有些论文会把效率和负载率拟合成分段线性函数增加求解难度但更贴近实际。氢燃料电池输入氢气输出电和热。电效率约0.45~0.55热回收率约0.3~0.4。模型为P_elec η_fc * P_H2_in产热率为H_fc η_heat_fc * P_H2_in。氨合成塔输入氢气和氮气输出氨。化学计量关系为1mol氨需要1.5mol氢和0.5mol氮。按能量计算时可以用能量转换效率0.7左右简化。氨裂解器输入氨输出氢。能量效率约0.7~0.8同样是稳态效率模型。燃气轮机如果考虑掺氨掺氢输入燃料可以是混合气输出是电和热。简化处理可以固定综合效率复杂处理则需要按燃料热值折算流量并加上掺混比例约束。储氢罐/储氨罐模型为SOC(t1) SOC(t) η_charge * P_in(t) - P_out(t)/η_discharge还要加容量上下限、单位时段充放速率上限、以及调度周期末存储量恢复约束循环调度常用。这里有一个重要提醒建模时优先用线性约束。效率常数是最简单的方式优先使用必须用变量的乘积时尽量避免两个连续变量相乘因为这会产生双线性项导致模型从MILP退化成MINLP求解难度剧增。后面我会专门讲线性化处理。2.3 能量平衡方程的书写规范综合能源系统调度里最容易写错、也最容易被审稿人追问的就是各能量网络的平衡约束。以我上面的系统为例会写成这样一组等式电平衡风电 光伏 燃料电池 燃气轮机 上级购电 蓄电池放电 电负荷 电解槽 电锅炉 蓄电池充电。每一项都要除以单位时间步长统一到MW或kW注意量纲。热平衡燃料电池余热 燃气轮机余热 燃气锅炉 蓄热罐放热 热负荷 蓄热罐蓄热。氢平衡电解槽产氢 氨裂解器产氢 储氢罐放氢 氢负荷 燃料电池耗氢 氨合成塔耗氢 储氢罐充氢。氨平衡氨合成塔产氨 储氨罐放氨 氨负荷 燃气轮机掺氨 氨裂解器耗氨 储氨罐充氨。这些等式写成Matlab代码后就是Yalmip里的Constraints [Constraints, P_wind P_pv P_fc P_gt P_buy P_bat_dis P_load P_elec P_eb P_bat_ch];这样一行行约束。重点在于等式里每个变量的符号方向必须一致一旦某个设备把方向搞反求解器给出的结果会非常诡异——比如系统明明电不够它却把电存进蓄电池。3. 优化调度模型的目标函数与约束体系决定模型聪明程度的部分3.1 目标函数经济性为主也可以叠加低碳指标纯经济调度目标函数通常是总运行成本最小化包括上级购电成本、天然气购气成本、设备运维成本按出力乘以系数、启停成本如果引入了启停变量减去售电收益如果允许系统向电网售电。在Matlab里这个目标函数写出来大概是这样Objective sum(buy_price .* P_buy) sum(gas_price .* G_gt gas_price .* G_gb) ... sum(om_elec * P_elec om_fc * P_fc om_gt * P_gt) ... sum(om_h2 * F_h2_store om_nh3 * F_nh3_store) ... sum(penalty_wind * P_wind_curtail); % 弃风惩罚项注意最后一项加入弃风惩罚系数是很有用的技巧。风光大发时如果不给一定惩罚模型可能选择大量弃风而不去制氢导致结果虽然成本最低但资源浪费严重。惩罚系数设置的量级要平衡——低于购电价高于设备运维成本这样模型才会在弃风和制氢存储之间做出合理权衡。如果研究倾向低碳化还可以把目标函数改为运行成本 碳交易成本。碳交易成本 系统碳排放量 × 碳价碳排放来源于购电和天然气消耗让模型在电/气之间做能源结构调整时自动考虑碳成本。有些论文还会做帕累托前沿分析把经济性和碳排放作为双目标用加权法或ε约束法求解这也是很有价值的扩展方向。3.2 约束体系的完整清单从工程实现角度看约束比目标函数更容易出错。下面这份清单是我在代码里一定会写的能量平衡约束上面提过的电、热、氢、氨四个平衡等式。设备出力上下限每个设备的出力变量有min和max限幅。燃气轮机还要有最小技术出力防止模型把出力压到接近零但效率极低的不合理运行点。爬坡速率约束-Ramp_down ≤ P(t1) - P(t) ≤ Ramp_up。这个约束在调度周期为1小时时对燃气轮机特别重要。储能约束容量上下限、充放功率上限、SOC递推关系、周期始末存储量一致约束。与上级网络的交互约束购电量上限联络线容量、购气量上限管网供气能力。电制氢/制氨的设备耦合约束如果电解槽同时向氢负荷和储氢罐供氢要保证供氢路径的总量关系。这通常已经体现在氢平衡等式里。启停逻辑约束可选如果燃气轮机、电解槽有启停成本需要引入二进制变量u_sta(t)并用大M法关联出力和启停状态。第7条是一个值得展开的地方。u(t)1表示设备运行那么出力满足u(t)*P_min ≤ P(t) ≤ u(t)*P_max同时把启停成本加入目标函数。对于电解槽这类频繁启停的设备有些模型还会加最小运行时间和最小停机时间约束这些都用二进制变量和大M法实现。这会显著增加求解难度但更真实。3.3 非线性项的处理MILP化是Matlab实现的关键很多用Matlab做优化的新手遇到的第一道坎就是原始模型里出现了非线性项Yalmip报错或者Gurobi拒绝求解。最常见的两类非线性项是双线性项比如氢气流量乘以氢气热值、两个连续变量相乘。分段线性函数比如设备效率随负载率变化、存储成本阶梯定价。处理双线性项的标准做法是离散化。例如要表达购买氢气的成本 购氢量 × 氢气价格如果氢价是固定常数那就没有问题如果氢价随购氢量分档就把购氢量拆成多个二进制变量对应的区间每个区间用固定单价。这样就把非线性目标替换成了混合整数线性目标。处理分段线性函数可以用Yalmip内置的pwf相关函数或者自己用大M法实现分段。但我的经验是除非论文对比需要否则在第一版模型里尽量用常数效率先把整条链路跑通。等代码确认无误再逐步把效率曲线换成分段线性观察结果变化。这样调试时不会因为非线性项引入的数值问题而难以排查。这里放一个判断标准若模型去掉所有整数变量后无论怎么缩放参数求解速度都极快秒级说明模型本身规模不大如果连续松弛模型也很难解那就要检查是不是有隐藏的双线性项或者病态约束。这个排查思路在调试阶段非常有用。4. Matlab代码实现从数据输入到结果输出的完整链路4.1 程序整体架构怎么组织代码才不让自己混乱Matlab做能源系统优化调度项目代码组织决定了后期改参数、加约束、出图分析的效率。我目前最常用的文件结构是这样main.m主脚本负责读取数据、设置时段、调用建模和求解、输出结果。data_parameters.m所有系统参数集中在这里包括设备容量、效率、价格、负荷曲线、风光出力曲线。build_system_model.m一个函数输入参数结构体输出Yalmip变量和约束集。solve_and_postprocess.m求解并整理结果输出各设备出力序列、储能SOC序列、成本明细。plot_results.m画图函数输出电/热/氢/氨平衡图、储能SOC图、源荷平衡图。这种分层的组织方式跟代码量大小没有关系哪怕只是一个24小时算例也建议从第一次写代码时就养成习惯。等你一个月后回来看代码或者导师要求换一组参数重跑时你会感谢当初这个决定。4.2 Yalmip建模核心代码解析一段可以直接改用的骨架下面是核心建模代码的简化版本我按照一个包含风电、光伏、电解槽、氢燃料电池、燃气轮机、氨合成塔、储氢罐的小系统来写。每个变量的注释都标了物理含义方便大家对照自己的系统修改。%% 参数设置 T 24; % 调度周期单位小时 dt 1; % 时间步长小时 % 读入数据假设已有结构体 paras % paras.P_wind_max, paras.P_pv_max: 各时段最大可再生出力 % paras.P_load, paras.H_load: 电负荷、热负荷 % paras.elec_price, paras.gas_price: 分时电价、气价 %% 定义变量 P_wind sdpvar(1, T); % 风电实际出力 P_pv sdpvar(1, T); % 光伏实际出力 P_elec sdpvar(1, T); % 电解槽输入电功率 P_gt sdpvar(1, T); % 燃气轮机输出电功率 H_gt sdpvar(1, T); % 燃气轮机输出热功率 F_h2 sdpvar(1, T); % 电解槽产氢功率kW按氢的高位热值折算 S_h2 sdpvar(1, T1); % 储氢罐存储量 P_buy sdpvar(1, T); % 上级购电功率 P_fc sdpvar(1, T); % 燃料电池发电功率 F_fc_in sdpvar(1, T); % 燃料电池耗氢功率 F_nh3 sdpvar(1, T); % 氨合成塔耗氢功率即合成氨消耗的氢 P_curt sdpvar(1, T); % 弃风弃光量 %% 目标函数运行成本最小 cost_buy sum(paras.elec_price .* P_buy); cost_gas sum(paras.gas_price .* (F_gt2gas * P_gt F_gb_gas)); % 燃气轮机耗气燃气锅炉耗气 cost_om sum(paras.om_elec * P_elec paras.om_fc * P_fc paras.om_gt * P_gt); cost_curt sum(paras.penalty_wind * P_curt); Objective cost_buy cost_gas cost_om cost_curt; %% 约束构建 Constraints []; % 电平衡 Constraints [Constraints, P_wind P_pv P_fc P_gt P_buy ... paras.P_load P_elec P_curt]; % 热平衡简化燃气轮机余热 燃气锅炉 Constraints [Constraints, H_gt H_gb paras.H_load]; % 氢平衡 Constraints [Constraints, F_h2 F_nh3_split paras.H2_load F_fc_in F_nh3]; % 风光出力约束 Constraints [Constraints, 0 P_wind paras.P_wind_max]; Constraints [Constraints, 0 P_pv paras.P_pv_max]; Constraints [Constraints, 0 P_curt paras.P_wind_max paras.P_pv_max]; % 电解槽约束 Constraints [Constraints, 0 P_elec paras.P_elec_max]; Constraints [Constraints, F_h2 paras.eta_elec * P_elec]; % 储氢罐SOC递推 Constraints [Constraints, S_h2(2:T1) S_h2(1:T) ... paras.eta_store * F_h2 - F_fc_in / paras.eta_fc]; Constraints [Constraints, 0 S_h2 paras.S_h2_max]; Constraints [Constraints, S_h2(1) paras.S_h2_initial]; Constraints [Constraints, S_h2(T1) paras.S_h2_initial]; % 周期性约束 % 燃气轮机 Constraints [Constraints, 0 P_gt paras.P_gt_max]; Constraints [Constraints, 0 H_gt paras.H_gt_max]; Constraints [Constraints, H_gt paras.heat_power_ratio * P_gt];这段代码是高度简化的骨架只保留了核心平衡关系和主要设备。实际项目中设备类别更多、约束更细但思路完全一样。我可以明确地说我这个项目里超过90%的bug都出在某个变量忘了加进平衡等式或者某个变量在等式里方向写反而不是出在求解器本身。4.3 求解器配置与求解策略YalmipGurobi最常见的黄金组合Yalmip是一个建模层本身不求解需要调用底层求解器。我的建议是MILP用Gurobi找不到Gurobi用CPLEX都没有就用MATLAB内置的intlinprog。Gurobi在教学用户里有免费许可证30天试用也很容易申请。求解代码很简单ops sdpsettings(solver, gurobi, verbose, 2, ... gurobi.TimeLimit, 300, ... gurobi.MIPGap, 0.005); optimize(Constraints, Objective, ops);solver指定求解器verbose2可以输出详细求解日志TimeLimit是限制求解时间MIPGap是终止最优性间隙。这几个参数在实际研究里非常重要如果模型规模很大不给时间限制Gurobi可能跑几个小时后才返回一个最优解这在做灵敏度分析时要反复求解几十次的时候是完全不能接受的。求解成功后结果提取方式P_wind_opt value(P_wind); P_gt_opt value(P_gt); S_h2_opt value(S_h2);注意如果optimize返回的sol里的sol.problem不等于0说明求解失败。最常见的失败类型是Infeasible不可行或Numerical issues数值问题。这时候不要急着调参数先用Yalmip的诊断命令check(Constraints)查看哪条约束不满足再用sol.solverinput.info看求解器层面的报错信息。4.4 结果可视化与指标评估让数据自己说话求解没问题后出图是给导师/审稿人/自己看懂结果的重要手段。我的出图清单一般包括电功率平衡堆叠图横轴是时间纵轴是功率把风电、光伏、购电、燃气轮机、燃料电池堆在一起叠加电负荷曲线可以直接看到哪些时段缺电、哪些时段有弃电。氢平衡图电解槽产氢、燃料电池耗氢、储氢罐SOC变化重点关注储氢罐的充放时段。储能SOC曲线图储氢罐、蓄电池的SOC变化直观反映系统的时序转移能力。成本构成饼图购电成本、购气成本、运维成本各占比例。画图用Matlab原生plot、stairs、area就够不需要额外工具箱。值得一提的小技巧是用stairs画阶梯图比plot更符合调度结果的分段常数特性视觉上更专业。另外在论文里出图时字体统一用Arial或Times New Roman坐标轴字号10.5pt或11pt线宽2左右导出为PDF矢量图比用截图清晰很多。5. 调试心得那些容易让模型跑不起来的坑5.1 单位不一致导致的量纲错误是最隐蔽的坑我在跑氢气氨气系统调度时遭遇过最让人崩溃的一次bug求解结果里储氢罐的SOC变化曲线特别奇怪明明是晚上电价低谷充了很多氢但第二天早上的储量却比预期少了一大截。排查了很久才发现是单位混用。电功率单位是MW但氢气那边我一开始用kg/h后来改成按能量算kW时漏掉了氢的高位热值系数结果储氢罐的进氢量比实际小了一个数量级。这类问题在单一能量系统里几乎不会出现因为大家都用MW但含氢氨的多能流系统里氢既有质量流量kg/h又有能量流量kW氨既有摩尔流量kmol/h又有能量流量混用几乎是必然事件。我的建议是在代码开头写一个unit_conversion函数把所有单位统一到一套基准上并且做一次手算校核——比如电解槽输入10MW电按效率0.7每小时的产氢能量应为7MW对应氢气质量约为0.21吨/小时。把这个数带入代码看到结果是否符合预期再继续下一步。5.2 不可行解定位不要大海捞针用逐条检查如果optimize返回不可行最常见的原因是约束之间存在矛盾。最典型的矛盾某个时段电负荷非常高但所有电源的出力上限加在一起都不够满足负荷。这在普通电力系统里很少见但在含氢氨系统里很容易因为电解槽抢电而出现——电解槽在低谷时段想多制氢但它把低谷时段本就不多的电都吃掉了导致系统必须从电网买高价电甚至出现电缺额。定位不可行约束的方法是把约束分成几组暂时去掉一组观察是否变可行。比如先只保留电平衡约束求解成功后再加入热平衡再成功就加入氢平衡。这样逐组加入很快能找到矛盾源头。Yalmip的check(Constraints)命令也能给出每一条约束的残差配合使用效果更好。5.3 求解器参数调优与性能问题MILP求解器在一个模型上卡死大多数时候不是求解器不行而是模型本身存在优化空间。我从实践中总结的性能调优优先级是去掉多余的整数变量。很多设备其实不需要二进制变量只要连续出力变量就能跑。比如电解槽能不能用连续变量表示启停如果能接受极低负荷效率差的误差就不要用0/1变量。收紧变量的上下界。设备出力上限是500MW但你只会在0~200MW区间运行那就直接把上限写成200MW。上下界越紧分支定界的效率越高。减少大M取值。大M法里的M值不要盲目取1e6太大容易引起数值问题取该约束可能出现的最大量级×1.2就够了。设置合理的MIPGap。做灵敏度分析时0.1%的最优性间隙已经足够精确不必追求0.01%以下的严格最优。问题规模实在太大时考虑滚动时域优化。把24小时分成四个6小时的窗口前一个窗口的终端储能状态作为下一个窗口的初始状态牺牲最优性换来几乎线性的求解时间节省。6. 从确定性调度到不确定性优化后续可以这样扩展6.1 源荷不确定性场景法与鲁棒优化怎么融入现有代码目前讲的全部是确定性调度——前提是风光出力、负荷都是已知精确值。但真实系统中预测误差不可避免。扩展方向主要有两个随机优化场景法对风光出力生成若干个典型场景模型在期望值意义上做决策。代码结构上只需要把原来的确定性变量扩展成场景索引维度例如P_gt(s, t)目标函数变为各场景成本乘概率再求和。约束里又分成两类第一类是可调变量电解槽、燃气轮机出力第二类是表达约束平衡等式表达约束需要在每个场景下满足。鲁棒优化不确定集合用区间表示比如风光出力在每个时段都在预测值的±15%范围内波动。模型目标是最坏情况下的成本最小化。Matlab里可以用Yalmip直接建模但需要比较熟悉对偶转化的数学推导。我个人觉得对一个刚把确定性代码跑通的研究者来说优先做场景法因为它的代码改动逻辑更直观审稿人也更容易理解鲁棒优化适合后续深入研究因为它对数学功底和求解器配置都有更高要求。6.2 多时间尺度协调日前调度日内调整氢和氨的引入天然带出了多时间尺度调度的需求。日前调度决定各设备的基础运行计划日内调整根据最新预测修正机组出力。实现上可以把两种时间尺度的决策变量分为两层日前变量是整数变量和基础连续变量日内变量是在日前计划基础上的修正量。这是一个比较高级的扩展建议前面的代码彻底跑通之后再做。6.3 电-氢-氨-碳耦合的完整闭环如果往更宏观的方向想含氢氨系统的调度还可以跟碳捕集、碳利用相结合——比如把天然气重整制氢产生的二氧化碳捕集下来用于合成甲醇或尿素氨合成的氮气也可以来自空分装置。整个系统就变成了一个电-氢-氨-碳多能流耦合的复杂网络优化调度问题会进一步复杂但研究价值也更高。Matlab代码实现的思路依然不变把每个设备看成能量转换节点把每个网络看成平衡约束用Yalmip建模用MILP求解用灵敏度分析验证结果稳定性。就我自己的切身体验来说做这类研究最难的不是数学推导也不是某个约束不会写而是保持对物理过程的直觉。每次求解完先别着急看成本数字逐个时段看一遍关键设备出力和储能SOC曲线问自己一个问题这个结果在物理上说得通吗如果风光大发时电解槽没开工、电价峰值时燃气轮机没满发、储氢罐到周期末莫名清空那肯定有地方建模出错了。代码是人写的bug在所难免但沿着物理直觉去排查一定比对着报错信息瞎猜要快得多。最后分享一个小操作经验在你终于把模型调通、结果也合理之后把边界条件比如储氢罐容量、电解槽效率往两个极端各改一次再跑一遍。如果结果变化符合物理预期说明模型稳健如果结果突然乱跳说明那组参数附近可能存在数值敏感区写论文时最好避开或者补充说明。这个习惯能帮你省下跟审稿人来回扯皮的大量时间。 含氢气氨气综合能源系统优化调度：Matlab建模与求解实践
做综合能源系统优化调度最怕的是什么不是模型复杂而是你辛辛苦苦建完模、跑完代码结果却不符合物理直觉——明明有大量弃风调度结果却在用电高峰期买高价电。我第一次在Matlab里跑含氢气氨气的综合能源系统优化调度时就栽过这个跟头后来才明白问题出在我们没有把氢和氨的时序平移特性真正建模进去。氢和氨跟电、热不一样电基本是即发即用热可以靠蓄热罐扛几个小时而氢和氨是真正能实现跨天、甚至跨周存储的化学储能介质。它们可以把风光大发时的多余电力转成化学能存起来等到负荷高峰再通过燃料电池或氨燃料机组发电相当于在物理层面给电网加了一个大容量充电宝。这套逻辑说起来简单真正落到Matlab代码里涉及设备建模、多能流平衡约束、混合整数变量处理、求解器调参等一系列问题。这篇博文就围绕含氢气氨气综合能源系统优化调度这个方向把建模思路、Matlab代码实现、调试心得和扩展方向一次性讲透。适合看这篇内容的读者正在做综合能源系统、氢能/氨能方向课题的研究生准备用Matlab做算例验证的工程师以及想了解氢氨耦合系统调度原理但对代码无从下手的初学者。我会尽量用能直接抄作业的方式写但也希望大家理解每一步背后的物理含义否则遇到求解不收敛、结果反直觉的时候排查起来会非常痛苦。1. 为什么要在综合能源系统里引入氢气氨气调度问题的本质变化1.1 氢与氨不是多了一个设备而是多了一条时间维度传统综合能源系统优化调度研究对象通常是电、热、气三种能量流。电力靠电网平衡热靠供热管网平衡天然气靠气网平衡。三者之间存在或松或紧的耦合燃气轮机同时产电产热热电联产电锅炉可以把电转成热燃气锅炉直接把气转成热。这种耦合是空间上的、同一时刻内的。氢气加入后问题立刻不一样了。电解槽把电转成氢氢气可以储存在高压储氢罐里也可以直接供给氢负荷还可以在燃料电池里重新转成电。这个电→氢→电的链路中间有存储环节所以它天然带时间属性我可以选择在风电最便宜的时候制氢在电价最高的晚高峰放氢发电。调度决策从这一时刻怎么分配能量变成了这一时刻是否应该存能量/放能量。氨作为氢的载体进一步延伸了这条链路。氨合成塔把氢和氮合成氨储氨罐比储氢罐更容易大规模存储氨裂解器又可以把氨重新分解成氢或者直接把氨送去掺烧发电。这样氢和氨就构成了一个可双向转换、可长期存储的能量走廊。从调度建模的角度看引入氢氨意味着模型中增加了储能状态变量储氢量、储氨量、反映转换效率的线性/非线性函数、以及一组跨时段耦合约束。1.2 氢氨耦合给调度带来的三个棘手问题第一个问题是多能流平衡约束的联立。电平衡、热平衡、氢平衡、氨平衡必须同时满足任何一个网络失衡其他网络都会间接受影响。比如电解槽开工不足氢产量下降不仅氢负荷受损氢燃料电池也发不出电最后可能连晚高峰的电平衡都被打破。第二个问题是储能时序耦合带来的计算复杂度上升。储氢罐和储氨罐的容量约束、进出能率约束、周期结束时的存储量约束都是跨时段约束。模型规模随调度周期线性增长24小时优化问题还比较轻松如果做成8760小时的年度生产模拟变量和约束数量会非常庞大对求解器和代码质量的要求明显提高。第三个问题是分段线性效率和启停逻辑让模型变成混合整数问题。电解槽在低负荷区间效率下降燃气轮机有最小技术出力储氢罐的充放状态可能涉及整数变量氨合成塔往往有连续运行时间约束。这导致整个优化模型通常是混合整数线性规划MILP或者混合整数非线性规划MINLP在Matlab里用Yalmip建模、再调Gurobi或CPLEX求解是最常见的做法也是我下面要详细展开的内容。2. 系统拓扑与设备建模先把能量流图变成数学关系2.1 一套典型的含氢氨综合能源系统长什么样在做Matlab实现之前必须先在纸上把系统拓扑画清楚。我在这类课题里最常用的一套结构如下供能侧风电机组、光伏阵列、上级电网购电、天然气网购气。转换侧电解槽电→氢、氢燃料电池氢→电热、燃气轮机天然气/氨/氢掺烧→电热、氨合成塔氢氮→氨、氨裂解器氨→氢、燃气锅炉气→热、电锅炉电→热。存储侧储氢罐、储氨罐、蓄电池、蓄热罐。负荷侧电负荷、热负荷、氢负荷、氨负荷比如工业用户的氨原料需求。从调度研究的角度每一个设备都是一个输入能量→输出能量的转换节点转换关系可以用效率常数或者分段线性效率曲线描述。这个系统明显比传统电热气三联供多了一层氢氨化学链但也正是这层化学链给了系统更大的调节灵活性——风光出力高了就多制氢制氨风光出力低了就放氢发电、氨裂解补氢。2.2 关键设备的稳态数学模型别一上来就写非线性对于24小时到一周的调度问题通常采用稳态模型忽略设备内部动态过程只关心单位时段内的能量转换关系。下面是我的常用建模方式风光出力直接取预测功率序列或者用场景法生成若干组典型场景。如果做鲁棒优化还会引入不确定区间。电解槽输入电功率输出氢能按热值或按质量。模型为P_H2_out η_elec * P_elec_in其中η_elec通常取0.6~0.75。有些论文会把效率和负载率拟合成分段线性函数增加求解难度但更贴近实际。氢燃料电池输入氢气输出电和热。电效率约0.45~0.55热回收率约0.3~0.4。模型为P_elec η_fc * P_H2_in产热率为H_fc η_heat_fc * P_H2_in。氨合成塔输入氢气和氮气输出氨。化学计量关系为1mol氨需要1.5mol氢和0.5mol氮。按能量计算时可以用能量转换效率0.7左右简化。氨裂解器输入氨输出氢。能量效率约0.7~0.8同样是稳态效率模型。燃气轮机如果考虑掺氨掺氢输入燃料可以是混合气输出是电和热。简化处理可以固定综合效率复杂处理则需要按燃料热值折算流量并加上掺混比例约束。储氢罐/储氨罐模型为SOC(t1) SOC(t) η_charge * P_in(t) - P_out(t)/η_discharge还要加容量上下限、单位时段充放速率上限、以及调度周期末存储量恢复约束循环调度常用。这里有一个重要提醒建模时优先用线性约束。效率常数是最简单的方式优先使用必须用变量的乘积时尽量避免两个连续变量相乘因为这会产生双线性项导致模型从MILP退化成MINLP求解难度剧增。后面我会专门讲线性化处理。2.3 能量平衡方程的书写规范综合能源系统调度里最容易写错、也最容易被审稿人追问的就是各能量网络的平衡约束。以我上面的系统为例会写成这样一组等式电平衡风电 光伏 燃料电池 燃气轮机 上级购电 蓄电池放电 电负荷 电解槽 电锅炉 蓄电池充电。每一项都要除以单位时间步长统一到MW或kW注意量纲。热平衡燃料电池余热 燃气轮机余热 燃气锅炉 蓄热罐放热 热负荷 蓄热罐蓄热。氢平衡电解槽产氢 氨裂解器产氢 储氢罐放氢 氢负荷 燃料电池耗氢 氨合成塔耗氢 储氢罐充氢。氨平衡氨合成塔产氨 储氨罐放氨 氨负荷 燃气轮机掺氨 氨裂解器耗氨 储氨罐充氨。这些等式写成Matlab代码后就是Yalmip里的Constraints [Constraints, P_wind P_pv P_fc P_gt P_buy P_bat_dis P_load P_elec P_eb P_bat_ch];这样一行行约束。重点在于等式里每个变量的符号方向必须一致一旦某个设备把方向搞反求解器给出的结果会非常诡异——比如系统明明电不够它却把电存进蓄电池。3. 优化调度模型的目标函数与约束体系决定模型聪明程度的部分3.1 目标函数经济性为主也可以叠加低碳指标纯经济调度目标函数通常是总运行成本最小化包括上级购电成本、天然气购气成本、设备运维成本按出力乘以系数、启停成本如果引入了启停变量减去售电收益如果允许系统向电网售电。在Matlab里这个目标函数写出来大概是这样Objective sum(buy_price .* P_buy) sum(gas_price .* G_gt gas_price .* G_gb) ... sum(om_elec * P_elec om_fc * P_fc om_gt * P_gt) ... sum(om_h2 * F_h2_store om_nh3 * F_nh3_store) ... sum(penalty_wind * P_wind_curtail); % 弃风惩罚项注意最后一项加入弃风惩罚系数是很有用的技巧。风光大发时如果不给一定惩罚模型可能选择大量弃风而不去制氢导致结果虽然成本最低但资源浪费严重。惩罚系数设置的量级要平衡——低于购电价高于设备运维成本这样模型才会在弃风和制氢存储之间做出合理权衡。如果研究倾向低碳化还可以把目标函数改为运行成本 碳交易成本。碳交易成本 系统碳排放量 × 碳价碳排放来源于购电和天然气消耗让模型在电/气之间做能源结构调整时自动考虑碳成本。有些论文还会做帕累托前沿分析把经济性和碳排放作为双目标用加权法或ε约束法求解这也是很有价值的扩展方向。3.2 约束体系的完整清单从工程实现角度看约束比目标函数更容易出错。下面这份清单是我在代码里一定会写的能量平衡约束上面提过的电、热、氢、氨四个平衡等式。设备出力上下限每个设备的出力变量有min和max限幅。燃气轮机还要有最小技术出力防止模型把出力压到接近零但效率极低的不合理运行点。爬坡速率约束-Ramp_down ≤ P(t1) - P(t) ≤ Ramp_up。这个约束在调度周期为1小时时对燃气轮机特别重要。储能约束容量上下限、充放功率上限、SOC递推关系、周期始末存储量一致约束。与上级网络的交互约束购电量上限联络线容量、购气量上限管网供气能力。电制氢/制氨的设备耦合约束如果电解槽同时向氢负荷和储氢罐供氢要保证供氢路径的总量关系。这通常已经体现在氢平衡等式里。启停逻辑约束可选如果燃气轮机、电解槽有启停成本需要引入二进制变量u_sta(t)并用大M法关联出力和启停状态。第7条是一个值得展开的地方。u(t)1表示设备运行那么出力满足u(t)*P_min ≤ P(t) ≤ u(t)*P_max同时把启停成本加入目标函数。对于电解槽这类频繁启停的设备有些模型还会加最小运行时间和最小停机时间约束这些都用二进制变量和大M法实现。这会显著增加求解难度但更真实。3.3 非线性项的处理MILP化是Matlab实现的关键很多用Matlab做优化的新手遇到的第一道坎就是原始模型里出现了非线性项Yalmip报错或者Gurobi拒绝求解。最常见的两类非线性项是双线性项比如氢气流量乘以氢气热值、两个连续变量相乘。分段线性函数比如设备效率随负载率变化、存储成本阶梯定价。处理双线性项的标准做法是离散化。例如要表达购买氢气的成本 购氢量 × 氢气价格如果氢价是固定常数那就没有问题如果氢价随购氢量分档就把购氢量拆成多个二进制变量对应的区间每个区间用固定单价。这样就把非线性目标替换成了混合整数线性目标。处理分段线性函数可以用Yalmip内置的pwf相关函数或者自己用大M法实现分段。但我的经验是除非论文对比需要否则在第一版模型里尽量用常数效率先把整条链路跑通。等代码确认无误再逐步把效率曲线换成分段线性观察结果变化。这样调试时不会因为非线性项引入的数值问题而难以排查。这里放一个判断标准若模型去掉所有整数变量后无论怎么缩放参数求解速度都极快秒级说明模型本身规模不大如果连续松弛模型也很难解那就要检查是不是有隐藏的双线性项或者病态约束。这个排查思路在调试阶段非常有用。4. Matlab代码实现从数据输入到结果输出的完整链路4.1 程序整体架构怎么组织代码才不让自己混乱Matlab做能源系统优化调度项目代码组织决定了后期改参数、加约束、出图分析的效率。我目前最常用的文件结构是这样main.m主脚本负责读取数据、设置时段、调用建模和求解、输出结果。data_parameters.m所有系统参数集中在这里包括设备容量、效率、价格、负荷曲线、风光出力曲线。build_system_model.m一个函数输入参数结构体输出Yalmip变量和约束集。solve_and_postprocess.m求解并整理结果输出各设备出力序列、储能SOC序列、成本明细。plot_results.m画图函数输出电/热/氢/氨平衡图、储能SOC图、源荷平衡图。这种分层的组织方式跟代码量大小没有关系哪怕只是一个24小时算例也建议从第一次写代码时就养成习惯。等你一个月后回来看代码或者导师要求换一组参数重跑时你会感谢当初这个决定。4.2 Yalmip建模核心代码解析一段可以直接改用的骨架下面是核心建模代码的简化版本我按照一个包含风电、光伏、电解槽、氢燃料电池、燃气轮机、氨合成塔、储氢罐的小系统来写。每个变量的注释都标了物理含义方便大家对照自己的系统修改。%% 参数设置 T 24; % 调度周期单位小时 dt 1; % 时间步长小时 % 读入数据假设已有结构体 paras % paras.P_wind_max, paras.P_pv_max: 各时段最大可再生出力 % paras.P_load, paras.H_load: 电负荷、热负荷 % paras.elec_price, paras.gas_price: 分时电价、气价 %% 定义变量 P_wind sdpvar(1, T); % 风电实际出力 P_pv sdpvar(1, T); % 光伏实际出力 P_elec sdpvar(1, T); % 电解槽输入电功率 P_gt sdpvar(1, T); % 燃气轮机输出电功率 H_gt sdpvar(1, T); % 燃气轮机输出热功率 F_h2 sdpvar(1, T); % 电解槽产氢功率kW按氢的高位热值折算 S_h2 sdpvar(1, T1); % 储氢罐存储量 P_buy sdpvar(1, T); % 上级购电功率 P_fc sdpvar(1, T); % 燃料电池发电功率 F_fc_in sdpvar(1, T); % 燃料电池耗氢功率 F_nh3 sdpvar(1, T); % 氨合成塔耗氢功率即合成氨消耗的氢 P_curt sdpvar(1, T); % 弃风弃光量 %% 目标函数运行成本最小 cost_buy sum(paras.elec_price .* P_buy); cost_gas sum(paras.gas_price .* (F_gt2gas * P_gt F_gb_gas)); % 燃气轮机耗气燃气锅炉耗气 cost_om sum(paras.om_elec * P_elec paras.om_fc * P_fc paras.om_gt * P_gt); cost_curt sum(paras.penalty_wind * P_curt); Objective cost_buy cost_gas cost_om cost_curt; %% 约束构建 Constraints []; % 电平衡 Constraints [Constraints, P_wind P_pv P_fc P_gt P_buy ... paras.P_load P_elec P_curt]; % 热平衡简化燃气轮机余热 燃气锅炉 Constraints [Constraints, H_gt H_gb paras.H_load]; % 氢平衡 Constraints [Constraints, F_h2 F_nh3_split paras.H2_load F_fc_in F_nh3]; % 风光出力约束 Constraints [Constraints, 0 P_wind paras.P_wind_max]; Constraints [Constraints, 0 P_pv paras.P_pv_max]; Constraints [Constraints, 0 P_curt paras.P_wind_max paras.P_pv_max]; % 电解槽约束 Constraints [Constraints, 0 P_elec paras.P_elec_max]; Constraints [Constraints, F_h2 paras.eta_elec * P_elec]; % 储氢罐SOC递推 Constraints [Constraints, S_h2(2:T1) S_h2(1:T) ... paras.eta_store * F_h2 - F_fc_in / paras.eta_fc]; Constraints [Constraints, 0 S_h2 paras.S_h2_max]; Constraints [Constraints, S_h2(1) paras.S_h2_initial]; Constraints [Constraints, S_h2(T1) paras.S_h2_initial]; % 周期性约束 % 燃气轮机 Constraints [Constraints, 0 P_gt paras.P_gt_max]; Constraints [Constraints, 0 H_gt paras.H_gt_max]; Constraints [Constraints, H_gt paras.heat_power_ratio * P_gt];这段代码是高度简化的骨架只保留了核心平衡关系和主要设备。实际项目中设备类别更多、约束更细但思路完全一样。我可以明确地说我这个项目里超过90%的bug都出在某个变量忘了加进平衡等式或者某个变量在等式里方向写反而不是出在求解器本身。4.3 求解器配置与求解策略YalmipGurobi最常见的黄金组合Yalmip是一个建模层本身不求解需要调用底层求解器。我的建议是MILP用Gurobi找不到Gurobi用CPLEX都没有就用MATLAB内置的intlinprog。Gurobi在教学用户里有免费许可证30天试用也很容易申请。求解代码很简单ops sdpsettings(solver, gurobi, verbose, 2, ... gurobi.TimeLimit, 300, ... gurobi.MIPGap, 0.005); optimize(Constraints, Objective, ops);solver指定求解器verbose2可以输出详细求解日志TimeLimit是限制求解时间MIPGap是终止最优性间隙。这几个参数在实际研究里非常重要如果模型规模很大不给时间限制Gurobi可能跑几个小时后才返回一个最优解这在做灵敏度分析时要反复求解几十次的时候是完全不能接受的。求解成功后结果提取方式P_wind_opt value(P_wind); P_gt_opt value(P_gt); S_h2_opt value(S_h2);注意如果optimize返回的sol里的sol.problem不等于0说明求解失败。最常见的失败类型是Infeasible不可行或Numerical issues数值问题。这时候不要急着调参数先用Yalmip的诊断命令check(Constraints)查看哪条约束不满足再用sol.solverinput.info看求解器层面的报错信息。4.4 结果可视化与指标评估让数据自己说话求解没问题后出图是给导师/审稿人/自己看懂结果的重要手段。我的出图清单一般包括电功率平衡堆叠图横轴是时间纵轴是功率把风电、光伏、购电、燃气轮机、燃料电池堆在一起叠加电负荷曲线可以直接看到哪些时段缺电、哪些时段有弃电。氢平衡图电解槽产氢、燃料电池耗氢、储氢罐SOC变化重点关注储氢罐的充放时段。储能SOC曲线图储氢罐、蓄电池的SOC变化直观反映系统的时序转移能力。成本构成饼图购电成本、购气成本、运维成本各占比例。画图用Matlab原生plot、stairs、area就够不需要额外工具箱。值得一提的小技巧是用stairs画阶梯图比plot更符合调度结果的分段常数特性视觉上更专业。另外在论文里出图时字体统一用Arial或Times New Roman坐标轴字号10.5pt或11pt线宽2左右导出为PDF矢量图比用截图清晰很多。5. 调试心得那些容易让模型跑不起来的坑5.1 单位不一致导致的量纲错误是最隐蔽的坑我在跑氢气氨气系统调度时遭遇过最让人崩溃的一次bug求解结果里储氢罐的SOC变化曲线特别奇怪明明是晚上电价低谷充了很多氢但第二天早上的储量却比预期少了一大截。排查了很久才发现是单位混用。电功率单位是MW但氢气那边我一开始用kg/h后来改成按能量算kW时漏掉了氢的高位热值系数结果储氢罐的进氢量比实际小了一个数量级。这类问题在单一能量系统里几乎不会出现因为大家都用MW但含氢氨的多能流系统里氢既有质量流量kg/h又有能量流量kW氨既有摩尔流量kmol/h又有能量流量混用几乎是必然事件。我的建议是在代码开头写一个unit_conversion函数把所有单位统一到一套基准上并且做一次手算校核——比如电解槽输入10MW电按效率0.7每小时的产氢能量应为7MW对应氢气质量约为0.21吨/小时。把这个数带入代码看到结果是否符合预期再继续下一步。5.2 不可行解定位不要大海捞针用逐条检查如果optimize返回不可行最常见的原因是约束之间存在矛盾。最典型的矛盾某个时段电负荷非常高但所有电源的出力上限加在一起都不够满足负荷。这在普通电力系统里很少见但在含氢氨系统里很容易因为电解槽抢电而出现——电解槽在低谷时段想多制氢但它把低谷时段本就不多的电都吃掉了导致系统必须从电网买高价电甚至出现电缺额。定位不可行约束的方法是把约束分成几组暂时去掉一组观察是否变可行。比如先只保留电平衡约束求解成功后再加入热平衡再成功就加入氢平衡。这样逐组加入很快能找到矛盾源头。Yalmip的check(Constraints)命令也能给出每一条约束的残差配合使用效果更好。5.3 求解器参数调优与性能问题MILP求解器在一个模型上卡死大多数时候不是求解器不行而是模型本身存在优化空间。我从实践中总结的性能调优优先级是去掉多余的整数变量。很多设备其实不需要二进制变量只要连续出力变量就能跑。比如电解槽能不能用连续变量表示启停如果能接受极低负荷效率差的误差就不要用0/1变量。收紧变量的上下界。设备出力上限是500MW但你只会在0~200MW区间运行那就直接把上限写成200MW。上下界越紧分支定界的效率越高。减少大M取值。大M法里的M值不要盲目取1e6太大容易引起数值问题取该约束可能出现的最大量级×1.2就够了。设置合理的MIPGap。做灵敏度分析时0.1%的最优性间隙已经足够精确不必追求0.01%以下的严格最优。问题规模实在太大时考虑滚动时域优化。把24小时分成四个6小时的窗口前一个窗口的终端储能状态作为下一个窗口的初始状态牺牲最优性换来几乎线性的求解时间节省。6. 从确定性调度到不确定性优化后续可以这样扩展6.1 源荷不确定性场景法与鲁棒优化怎么融入现有代码目前讲的全部是确定性调度——前提是风光出力、负荷都是已知精确值。但真实系统中预测误差不可避免。扩展方向主要有两个随机优化场景法对风光出力生成若干个典型场景模型在期望值意义上做决策。代码结构上只需要把原来的确定性变量扩展成场景索引维度例如P_gt(s, t)目标函数变为各场景成本乘概率再求和。约束里又分成两类第一类是可调变量电解槽、燃气轮机出力第二类是表达约束平衡等式表达约束需要在每个场景下满足。鲁棒优化不确定集合用区间表示比如风光出力在每个时段都在预测值的±15%范围内波动。模型目标是最坏情况下的成本最小化。Matlab里可以用Yalmip直接建模但需要比较熟悉对偶转化的数学推导。我个人觉得对一个刚把确定性代码跑通的研究者来说优先做场景法因为它的代码改动逻辑更直观审稿人也更容易理解鲁棒优化适合后续深入研究因为它对数学功底和求解器配置都有更高要求。6.2 多时间尺度协调日前调度日内调整氢和氨的引入天然带出了多时间尺度调度的需求。日前调度决定各设备的基础运行计划日内调整根据最新预测修正机组出力。实现上可以把两种时间尺度的决策变量分为两层日前变量是整数变量和基础连续变量日内变量是在日前计划基础上的修正量。这是一个比较高级的扩展建议前面的代码彻底跑通之后再做。6.3 电-氢-氨-碳耦合的完整闭环如果往更宏观的方向想含氢氨系统的调度还可以跟碳捕集、碳利用相结合——比如把天然气重整制氢产生的二氧化碳捕集下来用于合成甲醇或尿素氨合成的氮气也可以来自空分装置。整个系统就变成了一个电-氢-氨-碳多能流耦合的复杂网络优化调度问题会进一步复杂但研究价值也更高。Matlab代码实现的思路依然不变把每个设备看成能量转换节点把每个网络看成平衡约束用Yalmip建模用MILP求解用灵敏度分析验证结果稳定性。就我自己的切身体验来说做这类研究最难的不是数学推导也不是某个约束不会写而是保持对物理过程的直觉。每次求解完先别着急看成本数字逐个时段看一遍关键设备出力和储能SOC曲线问自己一个问题这个结果在物理上说得通吗如果风光大发时电解槽没开工、电价峰值时燃气轮机没满发、储氢罐到周期末莫名清空那肯定有地方建模出错了。代码是人写的bug在所难免但沿着物理直觉去排查一定比对着报错信息瞎猜要快得多。最后分享一个小操作经验在你终于把模型调通、结果也合理之后把边界条件比如储氢罐容量、电解槽效率往两个极端各改一次再跑一遍。如果结果变化符合物理预期说明模型稳健如果结果突然乱跳说明那组参数附近可能存在数值敏感区写论文时最好避开或者补充说明。这个习惯能帮你省下跟审稿人来回扯皮的大量时间。