氨储能最近在综合能源系统里讨论度很高我这次把一整套“电转氨风光火双层优化调度”的Matlab代码完整跑了一遍。这套模型适合做园区级或区域级能源规划的同学参考同样的逻辑换成氢储能、锂电池储能只需要替换储能单元约束就能复用。文章会从系统结构、模型约束、YalmipCplex的双层实现、以及我实际调试中踩过的坑讲起力求让你拿到思路后能直接照着搭。先说清楚这套代码解决的典型场景风电场和光伏电站出力波动大火电机组负责兜底但弃风弃光问题严重尤其是夜间风电大发时段电网消纳不掉。引入电转氨Power-to-Ammonia环节后富余电力通过电解水制氢再与空分氮气合成氨存入储氨罐等负荷高峰或者风光出力不足时再用氨燃气轮机发电或者走火电掺氨燃烧通道。本质上是用“氨”这个化学载体把电变成可长期储存、可跨季节调用的燃料氨储能同时具备能量密度高、储运便捷的特点这正是它这几年在行业里热度上升的原因。1. 整体设计思路为什么偏偏是氨又为什么双层1.1 系统拓扑里多了哪些环节先把我用的系统结构交代清楚。模型里不是只有“风机光伏火电”三台机组并联而是多了一套完整的“电转氨-储氨-氨发电”回路。整个网络可以分为四块供能侧、转换侧、储能侧和负荷侧。供能侧包括风电场、光伏电站、火电机组还可以接外网购电通道转换侧是电解槽电→氢、合成氨装置氢氮→氨、氨燃气轮机氨→电储能侧就是储氨罐做的是长时间尺度的氨量平衡负荷侧是常规电负荷代码里我按典型日96个时段每15分钟一个点处理。功率流向上有两条主线。第一条是传统电通道风光火直接送电给负荷不够再向外部电网购电第二条是化学储能通道负荷低谷或新能源大发时电功率进入电解槽和合成氨装置变成液氨存入储氨罐负荷高峰时储氨罐放氨给氨燃气轮机重新变成电。这套结构的关键在于火电机组和氨燃气轮机都能作为调节电源参与调度清洁能源、常规火电、化学储能三种角色同时在场双层优化调度才有意义。1.2 为什么氨储能比氢储能更适合综合能源我在分析模型之前做过一轮对比氢储能和氨储能经常被放在一起讨论但工程落脚点差别很大。氢的优点是能量转换链路短电解水出来就能直接用于燃料电池发电缺点也很明显氢气体积能量密度太低需要高压35MPa以上或者低温-253℃储运泄漏风险和安全成本高长距离运输非常麻烦。氨在常温下加压即可液化约1MPa、温度20℃左右就能保持液态储运条件比氢温和得多体积能量密度大约是压缩氢的1.5到2倍。从产业链成熟度看液氨的存储、输送、加注在化工行业已经很成熟港口和储罐设施可以改造复用。更实际的一点是氨本身就是绿氢的大规模载体电解水制氢后氢和氮气合成氨再把氨运到远端用户侧分解或者直接燃烧比直接运氢可行得多。当然氨储能不是没有代价。整个“电→氢→氨→电”的能量链条长往返效率测算下来只有25%到35%左右远低于锂电池储能的80%以上。所以氨储能不适合做短时高频调节它更适合做长时间尺度、跨季节的容量型储能。在这个模型里我把它定位成“备用容量新能源消纳工具”而不是替代火电去干AGC调频的活。1.3 双层优化把容量配置和运行调度拆开这套代码最大的特征是“双层优化调度”不是单层模型里同时决策容量和运行功率。我需要解释清楚为什么要这样拆。如果做成单层决策变量既包含电解槽额定功率、储氨罐容量这类“规划变量”又包含每个时段机组出力、氨流量等“运行变量”模型规模会迅速膨胀。96个时段乘以十几个变量再叠加0-1启停变量完全塞给求解器硬解混合整数非线性规划基本会让Cplex直接“罢工”而且上层规划的年度成本和下层运行的经济性目标混在一起量纲不一致权重系数调起来非常痛苦。双层的好处我个人理解主要是两点。第一是时间尺度解耦上层管“年”甚至“十年”的容量投资下层管“日内96点”的运行调度两者的决策频率和数学模型天然不同。第二是求解器分工明确上层用启发式算法我用的粒子群PSO搜索容量方案下层用YalmipCplex求解报给的容量参数下的最优运行成本内层返回的成本值构成外层适应度函数迭代前进。打个比方这就像做家庭预算先决定买多大的冰箱、装不装光伏板上层规划再根据每天天气和电价决定开关空调、什么时候洗衣服下层运行。两层决策分开算比一次性把所有问题都揉在一起要清楚得多。2. 核心模型与约束怎么建风-光-火-氨的数学表达2.1 风光出力与场景生成别让91个场景拖垮求解器风光出力的不确定性是这个模型绕不开的坎。我处理的方法是先产生大量随机场景再做场景削减。风速概率分布一般用Weibull分布光照强度用Beta分布通过逆变换采样可以生成出力的随机序列。蒙特卡洛方法跑500个场景下来每个场景都对应一组全天96点风电、光伏出力曲线。如果把这500个场景全部放进下层调度模型每轮双层迭代就要调用下层求解500次单次下层即使Cplex只花0.5秒一轮也要4分多钟50轮迭代下来就是三个多小时这个代价基本不能接受。我的做法是做K-means聚类削减。把500个场景聚成12类每一类取聚类中心作为典型场景权重设为该类的样本占比。这样下层调度只需处理12个场景计算量下降两个数量级精度损失却不大因为聚类中心保留了原始场景的主要波动特征和峰谷形状。代码里关键参数设置如下num_scenarios_initial 500; % 初始随机场景数 num_scenarios_reduced 12; % 削减后的典型场景数 K_opt 12; % 聚类数 [centers, weights] kmeans_reduce(initial_scenarios, K_opt);权重会进入下层目标函数每个典型场景的日运行成本乘该场景的权重累加起来就是全年运行成本的近似。这一步直接影响最终容量配置方案如果场景削减做得太粗比如只留3个场景峰谷特征丢了储氨罐容量会被明显低估。2.2 电转氨与氨储能的建模细节这部分是整个模型的核心需要把能量流拆开写清楚。电解槽模型我按线性效率处理输入电功率 (P_{PTH}(t))输出氢功率按低热值计为 (\eta_{PTH} \cdot P_{PTH}(t))(\eta_{PTH}) 取0.7。合成氨环节输入氢功率 (P_{H2}(t))输出氨功率 (P_{NH3_syn}(t)\eta_{PtA} \cdot P_{H2}(t))这个环节的 (\eta_{PtA}) 取0.9左右包含了氢气压缩、氮气制备、反应器辅助能耗的折算。氨燃气轮机的发电效率 (\eta_{PA}) 取0.4。我做过一个特别直观的测算把能量流换算成实际物理量你感受一下整个链路的效率。假设1MWh电进入电解槽效率0.7得到0.7MWh氢化学能换算成氢气质量大约是21kg这些氢气与氮气合成氨理论上能得到约119kg氨氨的低热值约为5.17kWh/kg所以119kg氨对应约0.62MWh氨化学能再经过40%效率的氨发电最终只能回送约0.25MWh电。整个“电-氢-氨-电”往返效率约25%这就是我前面说的氨储能适合长时储能、不适合短时高效回电的根本原因。储氨罐的动态模型是下层调度里最重要的约束之一[ E_{NH3}(t1) E_{NH3}(t) \Delta t \cdot (\eta_{fill} \cdot P_{NH3_syn}(t) - P_{NH3_dis}(t)) ]其中 (\eta_{fill}) 是注入效率(P_{NH3_dis}(t)) 是氨燃气轮机的氨消耗功率。储氨罐容量约束设置为[ E_{min} \leq E_{NH3}(t) \leq E_{max} ]注意一个细节这里的能量单位是MWh按氨的低热值折算不是立方米或吨。实际工程中有储罐的物理容积限制代码里 (E_{max}) 是通过氨的热值转换为等效能量存储量。我在初始版本里直接用了质量单位后来发现目标函数里能量和功率量纲不统一导致权重系数不断出问题最后统一改成MWh才理顺。这个坑后面会再细讲。2.3 火电机组与氨燃气轮机的运行约束火电机组在这个系统里承担的是基荷加调节的角色。它的建模不需要太复杂但要包含几个关键约束。出力上下限约束 (P_{fire}^{min} \le P_{fire}(t) \le P_{fire}^{max})爬坡约束 (|P_{fire}(t1)-P_{fire}(t)| \le R_{fire})以及最小启停时间约束。最小启停时间是0-1变量会引入混合整数性质。我习惯把火电出力离散化处理得太细因为每个时段都加0-1变量会让求解速度显著变慢。实操中我的做法是如果只研究经济调度而不研究机组组合就只保留出力上下限和爬坡约束把机组始终设为在线状态如果研究场景确实需要启停那就用Yalmip的binvar声明启停变量下层模型升级成MILP问题。氨燃气轮机的约束类似但多了氨燃料消耗与电出力之间的转换关系[ P_{PA}(t) \eta_{PA} \cdot P_{NH3_dis}(t) ]这里 (P_{NH3_dis}(t)) 就是储氨罐释放的氨功率(P_{PA}(t)) 就是氨燃气轮机的发电功率。氨燃气轮机最大电出力在代码里作为上层规划的决策变量之一也就是容量配置的结果。如果氨储能容量配置过大、回收成本过高上层粒子群会自动倾向于把氨轮机容量压到很低反之如果风光渗透率高、弃电惩罚重氨轮机容量就有存在的价值。这个博弈过程正是双层优化的核心魅力。3. Matlab代码实现从Yalmip建模到双层迭代3.1 工具箱选型与代码整体架构Matlab建模求解方面我推荐的组合是Yalmip建模Cplex求解器Yalmip由Johan Löfberg维护上层用PSO时不需要额外装优化工具箱底层求解只需要Cplex或Gurobi。Yalmip的好处是你能用接近数学表达式的语法描述约束代码可读性强、调试方便坏处是有时求解器识别不准确需要在sdpsettings里明确指定。我建议的代码文件结构如下main_upper.m % 主程序上层PSO inner_scheduler.m % 下层调度模型 create_scenarios.m % 场景生成与削减 data_load_profile.mat % 负荷数据 data_wind_pv.mat % 风光出力数据 params_initialize.m % 系统参数初始化3.2 上层粒子群优化的实现细节上层PSO的核心任务是在容量可行域里搜索一组使“年度综合成本最小”的方案。粒子维度对应四个决策变量电解槽额定功率、合成氨装置容量、储氨罐容量、氨燃气轮机容量。粒子群的位置和速度更新用标准公式v w * v c1 * rand * (pbest - x) c2 * rand * (gbest - x); x x v;权重w我还是习惯从0.9线性降到0.4前期重视全局搜索后期重视局部收敛c1、c2都取2.0粒子群规模设20个迭代40次。粒子变量的上下界从火电容量占系统峰值负荷比例来推算比如电解槽最大容量取系统峰值负荷的0.5倍储氨罐最大能量容量取单日负荷电量的0.3倍。每个粒子的适应度计算是先调用下层调度模块返回该容量方案下的全年运行成本再叠加容量投资折算到每年的成本。投资成本里我考虑了电解槽、合成氨装置、储氨罐按体积或能量容量、氨燃气轮机四类设备的单位投资系数和寿命折现。为了避免下层无解导致粒子适应度变成NaN我在内层设置了一个大罚值返回1e8。3.3 下层调度模型的Yalmip与Cplex实现下层调度模型是整个代码的引擎每个粒子都要调用一次。我把模型目标设成“典型日运行成本最小”包含购电成本、火电燃料成本、弃风弃光惩罚、碳交易成本和氨储能设备运行维护费用。Yalmip中声明决策变量T 96; P_wt sdpvar(1, T); P_pv sdpvar(1, T); P_fire sdpvar(1, T); P_PTH sdpvar(1, T); % 电解槽输入电功率 P_syn_NH3 sdpvar(1, T); % 合成氨输出氨功率 P_NH3_gen sdpvar(1, T); % 氨燃气轮机输出电功率 E_NH3 sdpvar(1, T1); % 储氨罐能量状态 P_grid sdpvar(1, T); % 购电功率正为购电负为售电功率平衡约束是整个模型的基石[ P_{wt}(t) P_{pv}(t) P_{fire}(t) P_{NH3_gen}(t) P_{grid}(t) P_{load}(t) P_{PTH}(t) ]每条约束在Yalmip里写起来非常直观Constraints []; for t 1:T Constraints [Constraints, P_wt(t) P_pv(t) P_fire(t) P_NH3_gen(t) P_grid(t) P_load(t) P_PTH(t)]; end目标函数可以直接用权重向量乘以变量数组启动Cplex求解Objective sum(C_grid .* P_grid) sum(C_fire .* P_fire) ... C_curtail * sum((P_wt_avail - P_wt) (P_pv_avail - P_pv)) ... C_carbon * emission_factor * sum(P_fire); ops sdpsettings(solver, cplex, verbose, 0); optimize(Constraints, Objective, ops);这里有一个我调试时发现的重要细节企业实际在写论文或报告时会把弃风弃光作为惩罚项放进目标但如果惩罚系数设得太小优化器会更倾向切掉新能源而不是调用氨储能导致氨储能设备闲置。我的处理是把弃电惩罚系数设为火电燃料成本的1.5倍左右让“弃电”变得足够贵这样优化器才有动力去充电转氨。3.4 双层迭代的数据传递与结果回传双层迭代是这套代码的主循环上层给出容量方案下层求解后返回运行成本上层再更新粒子位置直到达到最大迭代次数。我在实际跑的时候发现有人会把上层PSO的内层调用写成串行循环跑完40次迭代就要等上很久。建议用parfor做并行粒子计算但要注意内层调用的Cplex如果是共享许可证并行时每个worker都会占用一个license需要先确认你的Cplex许可证支持并行。如果没有并行许可证就把粒子群规模调小到16个迭代次数限制到30次总体耗时也可控。收敛判断也不能只靠最大迭代次数。我加了两个额外条件一是全局最优适应度在连续5代内的相对变化小于1e-4二是gbest位置在连续3代内不再变化。如果同时触发就提前退出迭代省下不必要的计算时间。最后的输出建议保存四个文件容量配置结果、典型日各设备出力曲线、储氨罐能量状态曲线、以及收敛曲线。4. 调试实录我踩过的坑和排查技巧4.1 Yalmip报错与求解器配置问题新手最常见的问题是Yalmip报“No suitable solver found”。这不是Yalmip本身的问题而是没有装求解器或者求解器路径没配好。我建议的检查顺序是yalmiptest运行这个命令能看到当前Yalmip能找到哪些求解器。如果Cplex被识别那大概率是许可证问题如果根本没出现Cplex先在Matlab里addpath到Cplex安装目录下的matlab子文件夹再重启Matlab。还有一个隐蔽问题是新版Cplex的Matlab接口文件夹路径是类似于C:\Program Files\IBM\ILOG\CPLEX_Studio2210\cplex\matlab\x64_win64的要addpath到这个层级而不是只加到cplex\matlab。如果遇到“MILP incompatible”之类的报错先检查你有没有使用binvar声明0-1变量然后检查Yalmip版本老版本对二进制变量和连续变量的混合约束支持偶尔会出现bug升级到GitHub上的最新版Yalmip能解决大部分问题。4.2 下层求解无解时的处理与罚值设置双层算法跑到中间内层Cplex返回无解是常有的事。原因通常是容量参数太小比如电解槽容量给到0.1倍峰值负荷但夜间风电大发时段可消纳功率超过这个容量功率平衡约束找不到可行解。我的处理方式是把内层优化从“硬约束”改成“软约束罚函数”负荷平衡约束加一个正负松弛变量目标函数里加一个极大罚系数。松弛变量哪怕被惩罚触发上层PSO也能得到可行解而不是返回NaN导致粒子群直接崩掉。罚系数建议取10的6次方以上这样正常可行解根本不会激活松弛变量但无解时会有一个有限数值反馈回上层粒子群仍有梯度信息可循。4.3 氨储能效率低导致优化器不用氨这是我在真实项目里遇到过的最要命的问题。运行结果出来后发现氨燃气轮机几乎不发电优化器宁可让火电满发甚至高价购电也不用储氨罐里的氨。原因非常简单氨储能“电→氢→氨→电”的往返效率只有25%而火电直接燃烧发电的效率有40%以上燃料成本上氨电根本打不过火电。解决方案有三个层次。第一层是在目标函数里加入碳交易成本让火电的碳排放成本显性化如果碳价足够高氨电的相对经济性就提升了。第二层是设一个“新能源消纳率约束”强制风光出力占比不低于某个阈值氨储能成为消纳的必需手段第三层是给氨燃气轮机加最小运行时长约束保证研究场景下氨储能至少被调用一定次数。如果你只是复现代码逻辑我建议直接用第二层约束既简单又贴合系统安全的实际需求。4.4 模型求解速度的压榨技巧整套代码最耗时的就是下层调度被反复调用Cplex每次都要重新构建模型并求解混合整数规划。我分享三个亲测有效的提速手段。第一个是控制场景数12个典型日场景已经足够支撑调度精度再多基本就是浪费。第二个是松弛求解精度在sdpsettings里给Cplex设置cplex.epgap, 0.01意思是Gap不超过1%就停运行时间能够缩短一半以上精度损失很小。第三个是避免在循环体内重复创建Yalmip变量先在外层一次性把所有sdpvar变量定义好内层只改系数矩阵和数据构建完Constraints和Objective后立即求解不要反复做变量声明。可选最后分享一个关于这套代码的个人体会我在实际跑通整个双层流程后发现容量配置结果对储氨罐初始容量特别敏感。初始 (E_{NH3}(1)) 设成满罐的一半和设成空罐得到的最终配置能差出30%以上。因为氨储能的价值就在于调节新能源的峰谷差初始库存不同长期平衡约束的平均水位也不同。建议你在代码里把初始氨量设成储罐容量的50%再把储罐末状态约束 (E_{NH3}(T1) E_{NH3}(1)) 加上让调度周期内氨量保持“周而复始”的稳态。这样才能避免“一次性耗尽存氨”的失真结果容量配置才会可信。这个细节不算复杂但直接影响模型结论的工程意义希望你能在代码里保留。
