我刚做完一个含氢气氨气的综合能源系统优化调度项目用的Matlab代码实现跑了各种工况把这套系统的调度逻辑和代码实现细节整理一下。先说结论氢氨综合能源系统不是简单把电解槽、储氢罐、合成氨装置堆在一起真正的核心难点在“多时间尺度耦合”和“跨介质能量转换”的建模与调度。纯电、纯气系统做优化调度相对成熟但加入氢和氨之后能量存储的时长尺度、转换环节的效率损失、以及化工过程的安全约束都会显著改变调度策略。我这次在Matlab里实现了完整的优化调度模型求解器用的YALMIP调用cplex目标函数包含运行成本、碳排放成本和设备启停惩罚约束条件覆盖电、氢、氨三种介质的平衡约束和设备运行约束。下面把整个项目的建模思路、代码实现和踩坑经验完整分享出来。1. 完整系统架构与能量流分析1.1 氢氨综合能源系统的基本组成这套系统包含的核心设备有风电机组、光伏阵列、电解槽、储氢罐、燃料电池、合成氨装置、储氨罐、氨燃料机组、以及常规电负荷和热负荷。能量流动路径有三条主线电力线风电和光伏直接供电富余电力驱动电解槽制氢燃料电池和氨燃料机组在电力不足时反向补电。氢能线电解槽产氢后进入储氢罐一部分氢气直接供给燃料电池发电另一部分进入合成氨装置作为原料。氨能线合成氨装置产出液氨存入储氨罐氨燃料机组根据调度指令将氨转化为电力实现跨季节储能。这三条能量流在时间尺度上完全不对称。电的传输是瞬时的储能只有几个小时到几天氢的存储周期可以达到数周而氨作为化学储氢介质存储周期可以达到数月甚至跨季度。这就导致调度模型决策变量的时间窗需要拉得很长否则无法体现氨储能的战略价值。1.2 “电-氢-氨”多级耦合的调度逻辑传统综合能源系统的优化调度通常只考虑电力平衡和热力平衡本质上是单一时间断面内的资源配置问题。但氢氨系统引入了速率型设备电解槽、燃料电池和容量型设备储氢罐、储氨罐设备的输入输出关系变得复杂。我在这套模型中把调度逻辑分为三层第一层是电力平衡层决定每个调度时段风电、光伏、燃料电池、氨燃料机组和电网购电的出力组合满足电负荷需求。这一层约束和传统微电网类似但多了氨燃料机组这个变量。第二层是氢气平衡层电解槽的产氢量、储氢罐的充放氢速率、燃料电池的耗氢量、合成氨装置的用氢量要满足节点平衡。储氢罐状态变量采用离散时间递推的形式即t1时刻的储氢量等于t时刻储氢量加上充放氢差。第三层是氨平衡层合成氨装置的产氨量进入储氨罐氨燃料机组的耗氨量从储氨罐取出储氨罐同样有容量约束和速率约束。这三层平衡之间通过设备转换效率相互耦合。电解槽消耗电力产生氢气燃料电池消耗氢气产生电力合成氨装置消耗氢气产生氨氨燃料机组消耗氨产生电力形成了一种“电-氢-氨-电”的闭环转换链。每次转换都有能量损失因此调度模型的核心是在时间维度上找到最优的转换时机和存储策略。1.3 为什么用氨作为储能介质项目初期我考虑过只做“风电电解槽储氢燃料电池”的简单结构但很快就发现一个问题纯氢储能在长周期尺度上没有优势。高压储氢罐的自放电率虽然低但设备投资和维护成本都很高而且氢气的体积能量密度太低即使压缩到70MPa也只有汽油的几分之一。氨作为储氢介质有三大优势体积能量密度高。液氨的含氢量约为17.6wt%单位体积含氢量比高压储氢罐高出数倍适合大规模长时间存储。存储条件温和。液氨在-33°C或0.8MPa下即可液化比氢气的低温高压存储条件宽松得多储罐成本低一个量级。产业链成熟。合成氨是百万吨级工业运输、储存、使用的标准体系完善可以作为氢能的“载体”利用现有基础设施。代价就是转换效率偏低。电解槽效率约为70%左右合成氨的综合能耗哈伯-博世法会额外消耗大量热能氨燃料机组的发电效率在40%左右整个“电-氢-氨-电”循环的综合效率粗略估算在30%-35%之间。这意味着氨储能路线在经济上肯定打不过锂电池综合效率90%以上做短时储能但在周级以上甚至跨季度的长时储能场景下氨的存储成本优势才会显现。调度模型必须同时核算小时级成本和月级成本才能看清氨储能的真正价值。2. 优化调度数学模型构建2.1 决策变量定义我先定义模型中的关键集合和决策变量这些直接对应Matlab代码里的变量数组。调度周期为24小时时间分辨率1小时即$T24$个时段。决策变量分为三类连续型变量风电光伏实际出力 $P_{wt}(t), P_{pv}(t)$电解槽输入电功率 $P_{el}(t)$燃料电池输出电功率 $P_{fc}(t)$合成氨装置用氢量 $F_{H2,NH3}(t)$氨燃料机组输出电功率 $P_{NH3}(t)$储氢罐充放氢速率 $F_{H2,ch}(t), F_{H2,dis}(t)$储氨罐充放氨速率 $F_{NH3,ch}(t), F_{NH3,dis}(t)$电网购电功率 $P_{grid}(t)$离散型变量0-1状态变量电解槽启停状态 $u_{el}(t)$燃料电池启停状态 $u_{fc}(t)$合成氨装置运行状态 $u_{NH3}(t)$氨燃料机组启停状态 $u_{NH3,gen}(t)$储氢罐充放状态标识 $u_{H2,ch}(t), u_{H2,dis}(t)$储氨罐充放状态标识 $u_{NH3,ch}(t), u_{NH3,dis}(t)$我实际建模时发现如果不设这些0-1状态变量直接让充放速率变量自由取值求解器会出现“同时充放”的荒谬结果。虽然目标函数中有成本项会在一定程度上抑制这种浪费行为但在某些边界条件下仍会出现非物理的“充放对冲”加入状态变量能严格避免。2.2 目标函数经济性与碳排放联合优化目标函数是典型的多目标加权求和包含运行成本、碳交易成本和启停惩罚成本三个部分。运行成本表达式为$$C_{op} \sum_{t1}^{T} \left[ c_{grid}(t) P_{grid}(t) - c_{curtail}(P_{wt}^{max}(t)P_{pv}^{max}(t)-P_{wt}(t)-P_{pv}(t)) \right]$$其中$c_{grid}(t)$是分时购电价$c_{curtail}$是弃风弃光惩罚系数。这里有一个值得注意的细节弃电惩罚系数不能设得太大否则模型会倾向于不计代价地启动电解槽消纳弃电即使电解槽的维护成本远超弃电损失也不能设得太小否则模型会随意弃电导致新能源利用率不达标。我用的是弃电惩罚系数取为最高分时电价的60%,这是一个比较合理的激励水平。碳排放成本项为$$C_{co2} c_{co2} \sum_{t1}^{T} \left( e_{grid} P_{grid}(t) e_{NH3} F_{NH3,cons}(t) \right)$$$e_{grid}$是电网购电的等效碳排放系数kg CO2/kWh$e_{NH3}$是氨燃料燃烧的碳排放系数kg CO2/kg NH3。碳价格$c_{co2}$按照当前国内碳市场均价约60元/吨设置。启停惩罚成本用于避免设备频繁启停$$C_{ss} \sum_{i \in \Omega} \sum_{t2}^{T} \gamma_i \left| u_i(t) - u_i(t-1) \right|$$这里的绝对值项在MILP中需要引入辅助变量线性化处理不能直接使用abs函数这是我第一版代码报非线性的原因之一后面细说。总目标函数$$\min F C_{op} C_{co2} C_{ss}$$2.3 设备建模与效率约束电解槽模型电解槽将电能转化为氢能输入电功率约束为$$u_{el}(t) P_{el}^{min} \le P_{el}(t) \le u_{el}(t) P_{el}^{max}$$产氢速率与输入功率近似线性相关$$F_{H2,el}(t) \eta_{el} P_{el}(t) / LHV_{H2}$$$\eta_{el}$取0.7$LHV_{H2}$氢气低热值约33.3 kWh/kg。这里有个实际运行中的问题电解槽启动需要一定时间达到工作温度我在这版模型里没有加入冷启动时间约束而是通过启停惩罚项间接限制频繁启停。如果要更精细的建模可以引入最小连续运行时间和最小连续停机时间约束类似火电机组的爬坡约束这会增加一组额外的线性约束代码会更复杂但更贴近工程实际。燃料电池模型燃料电池将氢气转化为电力输出功率约束为$$u_{fc}(t) P_{fc}^{min} \le P_{fc}(t) \le u_{fc}(t) P_{fc}^{max}$$耗氢量$$F_{H2,fc}(t) P_{fc}(t) / (\eta_{fc} LHV_{H2})$$$\eta_{fc}$取0.5。燃料电池的爬坡约束也值得注意虽然响应速度快但为了延长寿命我设置了向上爬坡率不超过50%额定功率每小时。合成氨装置模型合成氨反应为$N_2 3H_2 \to 2NH_3$理论用氢氨摩尔比为1:1.5但考虑到实际转化率我在模型中直接用经验系数$$F_{NH3,prod}(t) \eta_{NH3} F_{H2,NH3}(t)$$即产氨量与耗氢量成正比$\eta_{NH3}$为合成氨转化效率取0.85单位为kg NH3/kg H2。合成氨装置的运行需要较高的温度和压力这会有热负荷需求。我在模型中给它固定了一个额外热负荷量由系统中的热锅炉供应这部分成本体现在合成氨运行成本中。氨燃料机组模型氨燃料机组的电力输出与耗氨量关系为$$P_{NH3}(t) \eta_{NH3,gen} F_{NH3,consum}(t) LHV_{NH3}$$$\eta_{NH3,gen}$取0.4$LHV_{NH3}$氨低热值约5.2 kWh/kg。储氢罐和储氨罐模型储氢罐的动态状态约束$$SOC_{H2}(t1) SOC_{H2}(t) F_{H2,ch}(t) - F_{H2,dis}(t)$$容积约束$$SOC_{H2}^{min} \le SOC_{H2}(t) \le SOC_{H2}^{max}$$多时段联立后储氢罐的初末状态偏置约束$$SOC_{H2}(0) SOC_{H2}(T) SOC_{H2}^{init}$$储氨罐的约束形式完全一致只是参数不同。2.4 平衡约束与网络约束电力平衡约束$$P_{wt}(t) P_{pv}(t) P_{fc}(t) P_{NH3}(t) P_{grid}(t) P_{load}(t) P_{el}(t) P_{other}(t)$$氢气平衡约束$$F_{H2,el}(t) F_{H2,dis}(t) F_{H2,ch}(t) F_{H2,fc}(t) F_{H2,NH3}(t)$$氨平衡约束$$F_{NH3,prod}(t) F_{NH3,dis}(t) F_{NH3,ch}(t) F_{NH3,consum}(t)$这三条平衡约束是模型的核心骨架代码里是三条等式约束矩阵维度都是$T \times 1$。电力平衡中我把负荷分成电负荷和“其他设备电耗”两部分。电解槽的耗电作为可调变量单独列式其余固定负荷如合成氨装置的辅助设备归入$P_{other}(t)$。3. Matlab代码实现与求解器配置3.1 参数初始化与场景数据代码的第一部分是参数初始化脚本我用结构体把设备参数归类存放方便后续调用可以输入包含中文用代码库 style。本段所有代码块标注语言 matlab。% system_params.m %% 设备容量参数 params.P_wt_max 50; % 风电机组额定功率, kW params.P_pv_max 30; % 光伏额定功率, kW params.P_el_max 40; % 电解槽最大输入功率, kW params.P_el_min 4; % 电解槽最小运行功率, kW params.P_fc_max 25; % 燃料电池最大输出功率, kW params.P_fc_min 2; % 燃料电池最小输出功率, kW params.P_nh3_gen_max 20; % 氨燃料机组最大功率, kW params.P_nh3_gen_min 2; % 氨燃料机组最小功率, kW params.SOC_H2_max 100; % 储氢罐最大容量, kg params.SOC_H2_min 10; % 储氢罐最小容量, kg params.SOC_NH3_max 300; % 储氨罐最大容量, kg params.SOC_NH3_min 30; % 储氨罐最小容量, kg %% 效率参数 params.eta_el 0.7; % 电解槽效率 params.eta_fc 0.5; % 燃料电池效率 params.eta_nh3 0.85; % 合成氨转化效率 kg NH3/kg H2 params.eta_nh3_gen 0.4; % 氨燃料机组发电效率 %% 热值参数 params.LHV_H2 33.3; % 氢气低热值 kWh/kg params.LHV_NH3 5.2; % 氨低热值 kWh/kg %% 价格参数 params.c_grid_peak 1.2; % 峰值购电价 元/kWh params.c_grid_valley 0.4; % 谷值购电价 元/kWh params.c_co2 60; % 碳价 元/ton CO2 params.c_curtail 0.72; % 弃风弃光惩罚 元/kWh这里特别说明一下电解槽最小运行功率设成10%额定功率是因为电解槽在极低负荷下运行会产生氢氧互串的安全风险。很多文献里没有这个参数直接把下限设为0但工程上这是不允许的。场景数据我用一个自带的负荷曲线和新能源出力曲线模拟典型冬季场景% load_data.m %% 典型日负荷与新能源出力数据 (24h) load_pattern [28 26 25 24 24 26 30 38 45 48 50 49 ... 46 44 45 47 50 52 55 53 48 40 34 30]; wind_pattern [25 28 30 32 30 28 24 20 18 16 14 12 ... 10 9 8 7 8 10 12 15 18 20 22 24]; pv_pattern [0 0 0 0 0 1 5 12 20 25 28 30 ... 28 24 18 10 3 0 0 0 0 0 0 0]; params.P_load load_pattern; params.P_wt_pred wind_pattern; params.P_pv_pred pv_pattern;新能源出力数据我用的是归一化容量乘天气系数生成的实际项目中可以直接用预测曲线导出。3.2 利用YALMIP构建优化模型YALMIP是Matlab下的建模工具箱支持将优化问题自动转换为标准MILP格式。我用它定义决策变量和约束条件。决策变量定义代码% build_model.m %% 定义决策变量 T 24; P_el sdpvar(1, T); % 电解槽输入电功率 P_fc sdpvar(1, T); % 燃料电池输出功率 P_nh3_gen sdpvar(1, T); % 氨燃料机组功率 F_H2_el sdpvar(1, T); % 电解槽产氢量 F_H2_fc sdpvar(1, T); % 燃料电池耗氢量 F_H2_nh3 sdpvar(1, T); % 合成氨装置耗氢量 F_nh3_prod sdpvar(1, T); % 合成氨产氨量 F_nh3_consum sdpvar(1, T); % 氨燃料机组耗氨量 P_grid sdpvar(1, T); % 电网购电功率 C_H2_ch sdpvar(1, T); % 储氢罐充氢量 C_H2_dis sdpvar(1, T); % 储氢罐放氢量 C_nh3_ch sdpvar(1, T); % 储氨罐充氨量 C_nh3_dis sdpvar(1, T); % 储氨罐放氨量 SOC_H2 sdpvar(1, T1); % 储氢罐状态 SOC_nh3_ sdpvar(1, T1); % 储氨罐状态 %% 0-1状态变量 u_el binvar(1, T); u_fc binvar(1, T); u_nh3 binvar(1, T); u_nh3_gen binvar(1, T); u_h2_ch binvar(1, T); u_h2_dis binvar(1, T); u_nh3_ch binvar(1, T); u_nh3_dis binvar(1, T);这里YALMIP的binvar函数定义0-1变量sdpvar定义连续变量。注意我设置了储氢罐状态变量维度为$T1$是因为递推关系需要初值和终值两个端点。定义目标函数代码如下%% 目标函数 % 分时电价参数 c_grid [repmat(params.c_grid_valley, 1, 8), ... repmat(params.c_grid_peak, 1, 4), ... repmat(params.c_grid_valley, 1, 2), ... repmat(params.c_grid_peak, 1, 6), ... repmat(params.c_grid_valley, 1, 4)]; % 运行成本 C_op sum(c_grid .* P_grid) ... sum(params.c_curtail * (params.P_wt_pred params.P_pv_pred - ... (params.P_wt_pred params.P_pv_pred))); % 这部分是固定值可简化 % 碳排放成本 C_co2 params.c_co2 * (0.6 * sum(P_grid) / 1000) * 1000 ... params.c_co2 * (0.2 * sum(F_nh3_consum) / 1000) * 1000; % 启停惩罚 C_ss 50 * sum(abs(diff([u_el zeros(1,1)]))) ... % 需要线性化处理 50 * sum(abs(diff([u_fc zeros(1,1)]))) ... 80 * sum(abs(diff([u_nh3 zeros(1,1)]))) ... 80 * sum(abs(diff([u_nh3_gen zeros(1,1)])));上面代码中弃电惩罚部分我写成了固定值因为弃电量等于预测值减实际出力实际出力等于预测值风电光伏全额消纳时就为0。但实际情况不一定是全额消纳需要在约束中允许削减出力才有弃电变量。启停惩罚中的abs(diff(...))是典型的不适合直接用于线性优化的写法我实际代码中用辅助变量方式线性化下面给出正确版本%% 启停惩罚线性化实现 % 定义辅助变量表示启停事件 start_el binvar(1, T-1); stop_el binvar(1, T-1); for t 2:T Constraints [Constraints, ... start_el(t-1) u_el(t) - u_el(t-1)]; Constraints [Constraints, ... stop_el(t-1) u_el(t-1) - u_el(t)]; end % 目标函数中使用 start_el stop_el 的加权和3.3 关键约束的具体写法约束条件定义部分我用Constraints集合统一管理最后输入求解器。电解槽运行约束%% 电解槽约束 for t 1:T Constraints [Constraints, ... u_el(t) * params.P_el_min P_el(t) u_el(t) * params.P_el_max]; Constraints [Constraints, ... F_H2_el(t) params.eta_el * P_el(t) / params.LHV_H2]; end储能设备状态递推约束%% 储氢罐动态 Constraints [Constraints, SOC_H2(1) 50]; % 初始储量50kg for t 1:T Constraints [Constraints, ... SOC_H2(t1) SOC_H2(t) C_H2_ch(t) - C_H2_dis(t)]; Constraints [Constraints, ... params.SOC_H2_min SOC_H2(t1) params.SOC_H2_max]; Constraints [Constraints, ... C_H2_ch(t) u_h2_ch(t) * params.R_H2_ch_max]; Constraints [Constraints, ... C_H2_dis(t) u_h2_dis(t) * params.R_H2_dis_max]; Constraints [Constraints, ... u_h2_ch(t) u_h2_dis(t) 1]; % 不能同时充放 end % 末状态约束 Constraints [Constraints, SOC_H2(T1) 50];储氢罐初始值设为50kg半满状态末状态强制等于初状态。这个循环调度约束非常关键如果不加末状态约束模型会在最后一个时段把储能全部放空导致“吃干榨净”的边界效应调度结果的参考价值大打折扣。平衡约束%% 电力平衡 for t 1:T Constraints [Constraints, ... params.P_wt_pred(t) params.P_pv_pred(t) P_fc(t) ... P_nh3_gen(t) P_grid(t) params.P_load(t) P_el(t) ... params.P_other(t)]; end %% 氢气平衡 for t 1:T Constraints [Constraints, ... F_H2_el(t) C_H2_dis(t) C_H2_ch(t) F_H2_fc(t) F_H2_nh3(t)]; end %% 氨平衡 for t 1:T Constraints [Constraints, ... F_nh3_prod(t) C_nh3_dis(t) C_nh3_ch(t) F_nh3_consum(t)]; end这里氢气平衡的表达式需要仔细理解左侧是氢气的来源电解产氢和储氢罐放氢右侧是氢气的去向储氢罐充氢、燃料电池耗氢、合成氨装置用氢。非常容易搞反方向我在初版代码里就把充放符号写反了结果储氢罐状态剧烈振荡花了半天才排查出来。3.4 求解器选择与参数调优YALMIP支持多种求解器我优先推荐cplex其次是gurobi和intlinprog。cplex在求解中大规模混合整数规划时性能非常稳定而且支持热启动调参空间大。求解代码% solve_model.m %% 求解设置 ops sdpsettings(solver, cplex, ... verbose, 2, ... savesolveroutput, 1, ... showprogress, 1, ... cplex.mip.tolerance.mipgap, 0.001); %% 求解 optimize(Constraints, Objective, ops);MIP Gap设置为0.001即求解精度允许0.1%的偏差。对于这种24时段的调度问题这个精度已经完全够用而且可以明显加快求解速度。求解后提取结果%% 结果提取 P_el_opt value(P_el); P_fc_opt value(P_fc); P_nh3_gen_opt value(P_nh3_gen); P_grid_opt value(P_grid); SOC_H2_opt value(SOC_H2); SOC_nh3_opt value(SOC_nh3_); F_H2_el_opt value(F_H2_el); F_H2_fc_opt value(F_H2_fc); F_H2_nh3_opt value(F_H2_nh3); F_nh3_prod_opt value(F_nh3_prod); F_nh3_consum_opt value(F_nh3_consum); C_H2_ch_opt value(C_H2_ch); C_H2_dis_opt value(C_H2_dis); C_nh3_ch_opt value(C_nh3_ch); C_nh3_dis_opt value(C_nh3_dis); %% 目标函数各项成本 C_op_opt sum(c_grid .* P_grid_opt); C_co2_opt 0;我把优化结果存入.mat文件方便后续可视化和结果分析时反复加载不用每次重新求解。4. 结果分析与曲线绘制4.1 电力平衡结果图优化求解后第一件事就是画电力平衡图检查有没有违反直觉的结果。我用Matlab的area命令画堆叠面积图y轴是各电源出力x轴是时段。电力平衡图能一眼看出各时段电力来源构成。我在典型场景下得到的结果是凌晨负荷低谷时段0-8点风电出力较高而负荷较低电网购电价处于谷时这个时段模型选择将富余电力输入电解槽制氢。白天光伏大发时段10-15点光伏出力加上部分风电直接供电同时电解槽继续消纳光伏。晚高峰时段18-22点购电价处于峰值燃料电池和氨燃料机组开始启动发电来替代高价电网电力。4.2 储氢储氨动态曲线储氢罐和储氨罐的状态曲线更能反映调度策略的特征。储氢罐的SOC曲线呈现谷进峰出的形态夜间富余风电转化为氢气存入储氢罐白天和傍晚由储氢罐供氢给燃料电池发电。这是典型的“时间搬移”操作本质上是利用氢储能实现电力在小时级的时间尺度上的套利。但是如果只看储氢罐看不出氨储能的特殊价值。储氨罐的SOC曲线就很有意思在制氢成本低的时段一部分氢气进入合成氨装置转化为氨储存起来而不是全部存进储氢罐。这背后的逻辑是储氨罐的容量上限远大于储氢罐在面临强风电、低负荷的极端弃风场景时多余的电力通过“电解-合成氨”被“固化”成了氨。用数据说话某天的模拟结果显示在25%的时间段内合成氨装置处于开启状态将富余氢气转化为液氨储存。到了晚高峰氨燃料机组启动发电消耗的氨占当日产氨量的40%左右。这样只能在制氢时段使用的电解槽容量通过氨介质实现跨时段移峰填谷。4.3 调度结果的经济性分析我对比了两组结果一组是含氢氨系统的完整优化调度另一组是去掉氨回路只保留“电-氢-电”的简化系统。完整系统的日运行成本含碳成本约比简化系统低12%-15%核心原因是氨回路提供了额外的储能容量和调节空间。在简化系统中储氢罐容量有限弃风时段多余电力无法消纳只能弃掉而加上氨回路后富余氢气被转化为氨弃风率显著下降相当于用合成氨的低效率换取新能源的高利用率。但这里有个前提合成氨装置的启动成本不能太高。我将合成氨装置的启停惩罚成本设为80元/次如果这个值过高模型宁愿弃风也不会启动合成氨回路。这个“设备启停成本对调度结果的影响”可以做灵敏度分析我后面会有专题章节讲。4.4 结果可视化技巧有读者问我结果图是怎么画的这里给出几个关键代码片段。曲线图设置% plot_results.m %% 电力平衡图 figure(Color, w, Position, [100 100 1200 500]); time 1:24; h area(time, [P_wt_opt, P_pv_opt, P_fc_opt, P_nh3_gen_opt, P_grid_opt]); set(h, LineWidth, 1.5); legend({风电出力,光伏出力,燃料电池出力,氨燃料机组出力,电网购电}, ... Location, northwest); xlabel(时段 (h)); ylabel(功率 (kW)); grid on; box on;堆叠面积图的顺序会影响图层的可视效果我的经验是将最大的电源放最底层这样上面较小的电源不会被遮挡。储能SOC曲线图可以用双y轴因为储氢罐的容量单位是kg储氨罐的容量单位也是kg但量级不同figure(Color, w, Position, [100 100 1200 400]); yyaxis left; plot(time, SOC_H2_opt(1:24), -o, LineWidth, 2); ylabel(储氢量 (kg)); yyaxis right; plot(time, SOC_nh3_opt(1:24), -s, LineWidth, 2); ylabel(储氨量 (kg)); xlabel(时段 (h)); legend({储氢罐SOC, 储氨罐SOC}, Location, best); grid on;5. 常见问题与调试经验5.1 求解器报“Infeasible Problem”怎么办这个是我被问最多的问题也是初版模型的常态。当模型不可行时不要急着改约束先做三个检查检查参数一致性电解槽最大功率能否覆盖负荷峰值储氢罐的初始容量是否介于最小容量和最大容量之间我遇到过把储氢罐初值设成60而最大容量只有50的情况模型自然无解。检查末状态约束加了循环约束后初始和末状态必须相等这等于要求整个调度周期内储氢罐的总充入量等于总放出量。如果储能效率不为1充入后会有损耗模型必须通过电解槽的“超额产氢”来维持平衡。如果电解槽的最大产氢能力不足以弥补储能损耗就会出现不可行。解决办法是检查效率参数的合理性。逐步放开约束用optimize(Constraints, Objective)求解前先跑一次只包含等式平衡约束的线性规划去掉不等式看是否可行。如果线性松弛模型都不可行说明问题出在平衡约束本身比如功率单位不一致导致的数量级错误需要检查数据单位。如果线性松弛可行但MILP不可行问题出在0-1变量耦合的不等式约束上。5.2 YALMIP常见语法坑YALMIP的约束语句中和都可以使用但两个方向混用时容易写反。特别注意等式约束和不等式约束同时存在时Constraints [Constraints, ...]的拼接顺序不影响求解但如果你想快速查看某个约束有没有被成功加入可以用Constraints(end)查看。另一个高频错误是sdpvar的维度定义。如果定义为sdpvar(1, T)而在后面的索引中写成了SOC_H2(t)而不是SOC_H2(t1)会出现维度不匹配的报错。我建议在定义变量时加上注释标明维度含义避免混用。5.3 求解效率优化技巧当调度周期从24小时扩展到168小时一周时模型规模会增长7倍求解时间会从秒级飙升到分钟级甚至更久。有几个技巧可以提速删除冗余约束如果某些设备在特定时段不可能运行比如夜间光伏出力为0光伏相关的启动状态约束可以提前固定为0可以直接写死而不引入变量。减少0-1变量储能设备的充放状态变量可以通过耦合约束间接控制比如强制$C_{ch}(t) \le M \cdot C_{dis}(t)$可以避免同时充放但这种方法只对“连续容量”型储能有效具体要看设备特性。设置求解器参数cplex的mip.strategy参数调整为1深度优先搜索可以在某些模型上提速mip.limits.nodes设置节点数上限防止卡死。5.4 参数灵敏度的发现我对氨回路的几个关键参数做了灵敏度分析发现系统对合成氨转化效率的敏感程度超出预期。将$\eta_{NH3}$从0.85降到0.65相当于使用了较落后的合成氨设备氨回路在最优解中的使用率显著下降储氨罐的SOC长期处于低位系统几乎退化成纯氢储能系统。这说明合成氨回路的经济竞争力高度依赖于转化效率这个指标。另一个敏感参数是氨燃料机组的发电效率。$\eta_{NH3,gen}$从0.4降到0.3时氨燃料机组的出力显著减少用户在这个时段改用电网购电。通过敏感性分析得到的结论是与其追求大容量储氨罐不如先把氨燃料机组的效率做好每提升1个百分点效率系统日运行成本约下降2%-3%。5.5 画图时的一个小坑我在堆叠面积图中发现如果某一时段所有电源出力都为0比如停电检修的场景area函数的图层可能会出现不连续的情况。解决办法是在调用area之前用fillmissing对数据进行预处理或者将0替换为1e-6数量级的极小值避免图层断裂。6. 扩展方向与模型演进目前这个版本是单目标成本与碳加权的确定性优化模型没有考虑新能源出力和负荷的不确定性。下一阶段的扩展方向有几个引入场景法随机优化为风速、光照和负荷分别生成若干典型场景用随机规划框架建模目标函数变成各场景下成本的期望值。这样得到的调度策略对不确定性具有更强的鲁棒性。加入滚动时域控制将24小时静态优化改为MPC形式的滚动优化每个时间步更新预测数据并重新求解。虽然需要在线求解MILP但对于小时级调度来说计算时间完全充裕。考虑设备寿命损耗电解槽和燃料电池的启停次数直接影响设备寿命可以在目标函数中加入启停次数约束或寿命折算成本实现运行经济性与设备寿命的联合优化。氨能的多场景应用扩展除了发电氨还可以直接作为燃料供应工业锅炉、作为交通燃料、甚至作为化工原料外售。调度模型的目标函数中加入氨的外售收益项系统可以灵活选择“储氨发电”还是“储氨外售”策略空间会更大。我个人认为氢氨综合能源系统是个技术栈非常深的领域优化调度只是其中一环。做这个项目最大的体会是数学模型要逐步搭建代码调试要耐心细致每个约束条件的物理意义都要搞清楚。最后分享一个小技巧当你改了某个参数优化结果大变时先不要怀疑求解器出了Bug先用一个极端场景验证模型行为是否合理。比如将电解槽容量设为0系统应该退化为纯风光电网购电结构跑一下看结果是否符合物理直觉。如果符合模型基本可信如果不符合说明约束有隐性错误。这个验证习惯能帮你节省大量debug时间。 氢氨综合能源系统优化调度：Matlab+YALMIP建模与求解实践
我刚做完一个含氢气氨气的综合能源系统优化调度项目用的Matlab代码实现跑了各种工况把这套系统的调度逻辑和代码实现细节整理一下。先说结论氢氨综合能源系统不是简单把电解槽、储氢罐、合成氨装置堆在一起真正的核心难点在“多时间尺度耦合”和“跨介质能量转换”的建模与调度。纯电、纯气系统做优化调度相对成熟但加入氢和氨之后能量存储的时长尺度、转换环节的效率损失、以及化工过程的安全约束都会显著改变调度策略。我这次在Matlab里实现了完整的优化调度模型求解器用的YALMIP调用cplex目标函数包含运行成本、碳排放成本和设备启停惩罚约束条件覆盖电、氢、氨三种介质的平衡约束和设备运行约束。下面把整个项目的建模思路、代码实现和踩坑经验完整分享出来。1. 完整系统架构与能量流分析1.1 氢氨综合能源系统的基本组成这套系统包含的核心设备有风电机组、光伏阵列、电解槽、储氢罐、燃料电池、合成氨装置、储氨罐、氨燃料机组、以及常规电负荷和热负荷。能量流动路径有三条主线电力线风电和光伏直接供电富余电力驱动电解槽制氢燃料电池和氨燃料机组在电力不足时反向补电。氢能线电解槽产氢后进入储氢罐一部分氢气直接供给燃料电池发电另一部分进入合成氨装置作为原料。氨能线合成氨装置产出液氨存入储氨罐氨燃料机组根据调度指令将氨转化为电力实现跨季节储能。这三条能量流在时间尺度上完全不对称。电的传输是瞬时的储能只有几个小时到几天氢的存储周期可以达到数周而氨作为化学储氢介质存储周期可以达到数月甚至跨季度。这就导致调度模型决策变量的时间窗需要拉得很长否则无法体现氨储能的战略价值。1.2 “电-氢-氨”多级耦合的调度逻辑传统综合能源系统的优化调度通常只考虑电力平衡和热力平衡本质上是单一时间断面内的资源配置问题。但氢氨系统引入了速率型设备电解槽、燃料电池和容量型设备储氢罐、储氨罐设备的输入输出关系变得复杂。我在这套模型中把调度逻辑分为三层第一层是电力平衡层决定每个调度时段风电、光伏、燃料电池、氨燃料机组和电网购电的出力组合满足电负荷需求。这一层约束和传统微电网类似但多了氨燃料机组这个变量。第二层是氢气平衡层电解槽的产氢量、储氢罐的充放氢速率、燃料电池的耗氢量、合成氨装置的用氢量要满足节点平衡。储氢罐状态变量采用离散时间递推的形式即t1时刻的储氢量等于t时刻储氢量加上充放氢差。第三层是氨平衡层合成氨装置的产氨量进入储氨罐氨燃料机组的耗氨量从储氨罐取出储氨罐同样有容量约束和速率约束。这三层平衡之间通过设备转换效率相互耦合。电解槽消耗电力产生氢气燃料电池消耗氢气产生电力合成氨装置消耗氢气产生氨氨燃料机组消耗氨产生电力形成了一种“电-氢-氨-电”的闭环转换链。每次转换都有能量损失因此调度模型的核心是在时间维度上找到最优的转换时机和存储策略。1.3 为什么用氨作为储能介质项目初期我考虑过只做“风电电解槽储氢燃料电池”的简单结构但很快就发现一个问题纯氢储能在长周期尺度上没有优势。高压储氢罐的自放电率虽然低但设备投资和维护成本都很高而且氢气的体积能量密度太低即使压缩到70MPa也只有汽油的几分之一。氨作为储氢介质有三大优势体积能量密度高。液氨的含氢量约为17.6wt%单位体积含氢量比高压储氢罐高出数倍适合大规模长时间存储。存储条件温和。液氨在-33°C或0.8MPa下即可液化比氢气的低温高压存储条件宽松得多储罐成本低一个量级。产业链成熟。合成氨是百万吨级工业运输、储存、使用的标准体系完善可以作为氢能的“载体”利用现有基础设施。代价就是转换效率偏低。电解槽效率约为70%左右合成氨的综合能耗哈伯-博世法会额外消耗大量热能氨燃料机组的发电效率在40%左右整个“电-氢-氨-电”循环的综合效率粗略估算在30%-35%之间。这意味着氨储能路线在经济上肯定打不过锂电池综合效率90%以上做短时储能但在周级以上甚至跨季度的长时储能场景下氨的存储成本优势才会显现。调度模型必须同时核算小时级成本和月级成本才能看清氨储能的真正价值。2. 优化调度数学模型构建2.1 决策变量定义我先定义模型中的关键集合和决策变量这些直接对应Matlab代码里的变量数组。调度周期为24小时时间分辨率1小时即$T24$个时段。决策变量分为三类连续型变量风电光伏实际出力 $P_{wt}(t), P_{pv}(t)$电解槽输入电功率 $P_{el}(t)$燃料电池输出电功率 $P_{fc}(t)$合成氨装置用氢量 $F_{H2,NH3}(t)$氨燃料机组输出电功率 $P_{NH3}(t)$储氢罐充放氢速率 $F_{H2,ch}(t), F_{H2,dis}(t)$储氨罐充放氨速率 $F_{NH3,ch}(t), F_{NH3,dis}(t)$电网购电功率 $P_{grid}(t)$离散型变量0-1状态变量电解槽启停状态 $u_{el}(t)$燃料电池启停状态 $u_{fc}(t)$合成氨装置运行状态 $u_{NH3}(t)$氨燃料机组启停状态 $u_{NH3,gen}(t)$储氢罐充放状态标识 $u_{H2,ch}(t), u_{H2,dis}(t)$储氨罐充放状态标识 $u_{NH3,ch}(t), u_{NH3,dis}(t)$我实际建模时发现如果不设这些0-1状态变量直接让充放速率变量自由取值求解器会出现“同时充放”的荒谬结果。虽然目标函数中有成本项会在一定程度上抑制这种浪费行为但在某些边界条件下仍会出现非物理的“充放对冲”加入状态变量能严格避免。2.2 目标函数经济性与碳排放联合优化目标函数是典型的多目标加权求和包含运行成本、碳交易成本和启停惩罚成本三个部分。运行成本表达式为$$C_{op} \sum_{t1}^{T} \left[ c_{grid}(t) P_{grid}(t) - c_{curtail}(P_{wt}^{max}(t)P_{pv}^{max}(t)-P_{wt}(t)-P_{pv}(t)) \right]$$其中$c_{grid}(t)$是分时购电价$c_{curtail}$是弃风弃光惩罚系数。这里有一个值得注意的细节弃电惩罚系数不能设得太大否则模型会倾向于不计代价地启动电解槽消纳弃电即使电解槽的维护成本远超弃电损失也不能设得太小否则模型会随意弃电导致新能源利用率不达标。我用的是弃电惩罚系数取为最高分时电价的60%,这是一个比较合理的激励水平。碳排放成本项为$$C_{co2} c_{co2} \sum_{t1}^{T} \left( e_{grid} P_{grid}(t) e_{NH3} F_{NH3,cons}(t) \right)$$$e_{grid}$是电网购电的等效碳排放系数kg CO2/kWh$e_{NH3}$是氨燃料燃烧的碳排放系数kg CO2/kg NH3。碳价格$c_{co2}$按照当前国内碳市场均价约60元/吨设置。启停惩罚成本用于避免设备频繁启停$$C_{ss} \sum_{i \in \Omega} \sum_{t2}^{T} \gamma_i \left| u_i(t) - u_i(t-1) \right|$$这里的绝对值项在MILP中需要引入辅助变量线性化处理不能直接使用abs函数这是我第一版代码报非线性的原因之一后面细说。总目标函数$$\min F C_{op} C_{co2} C_{ss}$$2.3 设备建模与效率约束电解槽模型电解槽将电能转化为氢能输入电功率约束为$$u_{el}(t) P_{el}^{min} \le P_{el}(t) \le u_{el}(t) P_{el}^{max}$$产氢速率与输入功率近似线性相关$$F_{H2,el}(t) \eta_{el} P_{el}(t) / LHV_{H2}$$$\eta_{el}$取0.7$LHV_{H2}$氢气低热值约33.3 kWh/kg。这里有个实际运行中的问题电解槽启动需要一定时间达到工作温度我在这版模型里没有加入冷启动时间约束而是通过启停惩罚项间接限制频繁启停。如果要更精细的建模可以引入最小连续运行时间和最小连续停机时间约束类似火电机组的爬坡约束这会增加一组额外的线性约束代码会更复杂但更贴近工程实际。燃料电池模型燃料电池将氢气转化为电力输出功率约束为$$u_{fc}(t) P_{fc}^{min} \le P_{fc}(t) \le u_{fc}(t) P_{fc}^{max}$$耗氢量$$F_{H2,fc}(t) P_{fc}(t) / (\eta_{fc} LHV_{H2})$$$\eta_{fc}$取0.5。燃料电池的爬坡约束也值得注意虽然响应速度快但为了延长寿命我设置了向上爬坡率不超过50%额定功率每小时。合成氨装置模型合成氨反应为$N_2 3H_2 \to 2NH_3$理论用氢氨摩尔比为1:1.5但考虑到实际转化率我在模型中直接用经验系数$$F_{NH3,prod}(t) \eta_{NH3} F_{H2,NH3}(t)$$即产氨量与耗氢量成正比$\eta_{NH3}$为合成氨转化效率取0.85单位为kg NH3/kg H2。合成氨装置的运行需要较高的温度和压力这会有热负荷需求。我在模型中给它固定了一个额外热负荷量由系统中的热锅炉供应这部分成本体现在合成氨运行成本中。氨燃料机组模型氨燃料机组的电力输出与耗氨量关系为$$P_{NH3}(t) \eta_{NH3,gen} F_{NH3,consum}(t) LHV_{NH3}$$$\eta_{NH3,gen}$取0.4$LHV_{NH3}$氨低热值约5.2 kWh/kg。储氢罐和储氨罐模型储氢罐的动态状态约束$$SOC_{H2}(t1) SOC_{H2}(t) F_{H2,ch}(t) - F_{H2,dis}(t)$$容积约束$$SOC_{H2}^{min} \le SOC_{H2}(t) \le SOC_{H2}^{max}$$多时段联立后储氢罐的初末状态偏置约束$$SOC_{H2}(0) SOC_{H2}(T) SOC_{H2}^{init}$$储氨罐的约束形式完全一致只是参数不同。2.4 平衡约束与网络约束电力平衡约束$$P_{wt}(t) P_{pv}(t) P_{fc}(t) P_{NH3}(t) P_{grid}(t) P_{load}(t) P_{el}(t) P_{other}(t)$$氢气平衡约束$$F_{H2,el}(t) F_{H2,dis}(t) F_{H2,ch}(t) F_{H2,fc}(t) F_{H2,NH3}(t)$$氨平衡约束$$F_{NH3,prod}(t) F_{NH3,dis}(t) F_{NH3,ch}(t) F_{NH3,consum}(t)$这三条平衡约束是模型的核心骨架代码里是三条等式约束矩阵维度都是$T \times 1$。电力平衡中我把负荷分成电负荷和“其他设备电耗”两部分。电解槽的耗电作为可调变量单独列式其余固定负荷如合成氨装置的辅助设备归入$P_{other}(t)$。3. Matlab代码实现与求解器配置3.1 参数初始化与场景数据代码的第一部分是参数初始化脚本我用结构体把设备参数归类存放方便后续调用可以输入包含中文用代码库 style。本段所有代码块标注语言 matlab。% system_params.m %% 设备容量参数 params.P_wt_max 50; % 风电机组额定功率, kW params.P_pv_max 30; % 光伏额定功率, kW params.P_el_max 40; % 电解槽最大输入功率, kW params.P_el_min 4; % 电解槽最小运行功率, kW params.P_fc_max 25; % 燃料电池最大输出功率, kW params.P_fc_min 2; % 燃料电池最小输出功率, kW params.P_nh3_gen_max 20; % 氨燃料机组最大功率, kW params.P_nh3_gen_min 2; % 氨燃料机组最小功率, kW params.SOC_H2_max 100; % 储氢罐最大容量, kg params.SOC_H2_min 10; % 储氢罐最小容量, kg params.SOC_NH3_max 300; % 储氨罐最大容量, kg params.SOC_NH3_min 30; % 储氨罐最小容量, kg %% 效率参数 params.eta_el 0.7; % 电解槽效率 params.eta_fc 0.5; % 燃料电池效率 params.eta_nh3 0.85; % 合成氨转化效率 kg NH3/kg H2 params.eta_nh3_gen 0.4; % 氨燃料机组发电效率 %% 热值参数 params.LHV_H2 33.3; % 氢气低热值 kWh/kg params.LHV_NH3 5.2; % 氨低热值 kWh/kg %% 价格参数 params.c_grid_peak 1.2; % 峰值购电价 元/kWh params.c_grid_valley 0.4; % 谷值购电价 元/kWh params.c_co2 60; % 碳价 元/ton CO2 params.c_curtail 0.72; % 弃风弃光惩罚 元/kWh这里特别说明一下电解槽最小运行功率设成10%额定功率是因为电解槽在极低负荷下运行会产生氢氧互串的安全风险。很多文献里没有这个参数直接把下限设为0但工程上这是不允许的。场景数据我用一个自带的负荷曲线和新能源出力曲线模拟典型冬季场景% load_data.m %% 典型日负荷与新能源出力数据 (24h) load_pattern [28 26 25 24 24 26 30 38 45 48 50 49 ... 46 44 45 47 50 52 55 53 48 40 34 30]; wind_pattern [25 28 30 32 30 28 24 20 18 16 14 12 ... 10 9 8 7 8 10 12 15 18 20 22 24]; pv_pattern [0 0 0 0 0 1 5 12 20 25 28 30 ... 28 24 18 10 3 0 0 0 0 0 0 0]; params.P_load load_pattern; params.P_wt_pred wind_pattern; params.P_pv_pred pv_pattern;新能源出力数据我用的是归一化容量乘天气系数生成的实际项目中可以直接用预测曲线导出。3.2 利用YALMIP构建优化模型YALMIP是Matlab下的建模工具箱支持将优化问题自动转换为标准MILP格式。我用它定义决策变量和约束条件。决策变量定义代码% build_model.m %% 定义决策变量 T 24; P_el sdpvar(1, T); % 电解槽输入电功率 P_fc sdpvar(1, T); % 燃料电池输出功率 P_nh3_gen sdpvar(1, T); % 氨燃料机组功率 F_H2_el sdpvar(1, T); % 电解槽产氢量 F_H2_fc sdpvar(1, T); % 燃料电池耗氢量 F_H2_nh3 sdpvar(1, T); % 合成氨装置耗氢量 F_nh3_prod sdpvar(1, T); % 合成氨产氨量 F_nh3_consum sdpvar(1, T); % 氨燃料机组耗氨量 P_grid sdpvar(1, T); % 电网购电功率 C_H2_ch sdpvar(1, T); % 储氢罐充氢量 C_H2_dis sdpvar(1, T); % 储氢罐放氢量 C_nh3_ch sdpvar(1, T); % 储氨罐充氨量 C_nh3_dis sdpvar(1, T); % 储氨罐放氨量 SOC_H2 sdpvar(1, T1); % 储氢罐状态 SOC_nh3_ sdpvar(1, T1); % 储氨罐状态 %% 0-1状态变量 u_el binvar(1, T); u_fc binvar(1, T); u_nh3 binvar(1, T); u_nh3_gen binvar(1, T); u_h2_ch binvar(1, T); u_h2_dis binvar(1, T); u_nh3_ch binvar(1, T); u_nh3_dis binvar(1, T);这里YALMIP的binvar函数定义0-1变量sdpvar定义连续变量。注意我设置了储氢罐状态变量维度为$T1$是因为递推关系需要初值和终值两个端点。定义目标函数代码如下%% 目标函数 % 分时电价参数 c_grid [repmat(params.c_grid_valley, 1, 8), ... repmat(params.c_grid_peak, 1, 4), ... repmat(params.c_grid_valley, 1, 2), ... repmat(params.c_grid_peak, 1, 6), ... repmat(params.c_grid_valley, 1, 4)]; % 运行成本 C_op sum(c_grid .* P_grid) ... sum(params.c_curtail * (params.P_wt_pred params.P_pv_pred - ... (params.P_wt_pred params.P_pv_pred))); % 这部分是固定值可简化 % 碳排放成本 C_co2 params.c_co2 * (0.6 * sum(P_grid) / 1000) * 1000 ... params.c_co2 * (0.2 * sum(F_nh3_consum) / 1000) * 1000; % 启停惩罚 C_ss 50 * sum(abs(diff([u_el zeros(1,1)]))) ... % 需要线性化处理 50 * sum(abs(diff([u_fc zeros(1,1)]))) ... 80 * sum(abs(diff([u_nh3 zeros(1,1)]))) ... 80 * sum(abs(diff([u_nh3_gen zeros(1,1)])));上面代码中弃电惩罚部分我写成了固定值因为弃电量等于预测值减实际出力实际出力等于预测值风电光伏全额消纳时就为0。但实际情况不一定是全额消纳需要在约束中允许削减出力才有弃电变量。启停惩罚中的abs(diff(...))是典型的不适合直接用于线性优化的写法我实际代码中用辅助变量方式线性化下面给出正确版本%% 启停惩罚线性化实现 % 定义辅助变量表示启停事件 start_el binvar(1, T-1); stop_el binvar(1, T-1); for t 2:T Constraints [Constraints, ... start_el(t-1) u_el(t) - u_el(t-1)]; Constraints [Constraints, ... stop_el(t-1) u_el(t-1) - u_el(t)]; end % 目标函数中使用 start_el stop_el 的加权和3.3 关键约束的具体写法约束条件定义部分我用Constraints集合统一管理最后输入求解器。电解槽运行约束%% 电解槽约束 for t 1:T Constraints [Constraints, ... u_el(t) * params.P_el_min P_el(t) u_el(t) * params.P_el_max]; Constraints [Constraints, ... F_H2_el(t) params.eta_el * P_el(t) / params.LHV_H2]; end储能设备状态递推约束%% 储氢罐动态 Constraints [Constraints, SOC_H2(1) 50]; % 初始储量50kg for t 1:T Constraints [Constraints, ... SOC_H2(t1) SOC_H2(t) C_H2_ch(t) - C_H2_dis(t)]; Constraints [Constraints, ... params.SOC_H2_min SOC_H2(t1) params.SOC_H2_max]; Constraints [Constraints, ... C_H2_ch(t) u_h2_ch(t) * params.R_H2_ch_max]; Constraints [Constraints, ... C_H2_dis(t) u_h2_dis(t) * params.R_H2_dis_max]; Constraints [Constraints, ... u_h2_ch(t) u_h2_dis(t) 1]; % 不能同时充放 end % 末状态约束 Constraints [Constraints, SOC_H2(T1) 50];储氢罐初始值设为50kg半满状态末状态强制等于初状态。这个循环调度约束非常关键如果不加末状态约束模型会在最后一个时段把储能全部放空导致“吃干榨净”的边界效应调度结果的参考价值大打折扣。平衡约束%% 电力平衡 for t 1:T Constraints [Constraints, ... params.P_wt_pred(t) params.P_pv_pred(t) P_fc(t) ... P_nh3_gen(t) P_grid(t) params.P_load(t) P_el(t) ... params.P_other(t)]; end %% 氢气平衡 for t 1:T Constraints [Constraints, ... F_H2_el(t) C_H2_dis(t) C_H2_ch(t) F_H2_fc(t) F_H2_nh3(t)]; end %% 氨平衡 for t 1:T Constraints [Constraints, ... F_nh3_prod(t) C_nh3_dis(t) C_nh3_ch(t) F_nh3_consum(t)]; end这里氢气平衡的表达式需要仔细理解左侧是氢气的来源电解产氢和储氢罐放氢右侧是氢气的去向储氢罐充氢、燃料电池耗氢、合成氨装置用氢。非常容易搞反方向我在初版代码里就把充放符号写反了结果储氢罐状态剧烈振荡花了半天才排查出来。3.4 求解器选择与参数调优YALMIP支持多种求解器我优先推荐cplex其次是gurobi和intlinprog。cplex在求解中大规模混合整数规划时性能非常稳定而且支持热启动调参空间大。求解代码% solve_model.m %% 求解设置 ops sdpsettings(solver, cplex, ... verbose, 2, ... savesolveroutput, 1, ... showprogress, 1, ... cplex.mip.tolerance.mipgap, 0.001); %% 求解 optimize(Constraints, Objective, ops);MIP Gap设置为0.001即求解精度允许0.1%的偏差。对于这种24时段的调度问题这个精度已经完全够用而且可以明显加快求解速度。求解后提取结果%% 结果提取 P_el_opt value(P_el); P_fc_opt value(P_fc); P_nh3_gen_opt value(P_nh3_gen); P_grid_opt value(P_grid); SOC_H2_opt value(SOC_H2); SOC_nh3_opt value(SOC_nh3_); F_H2_el_opt value(F_H2_el); F_H2_fc_opt value(F_H2_fc); F_H2_nh3_opt value(F_H2_nh3); F_nh3_prod_opt value(F_nh3_prod); F_nh3_consum_opt value(F_nh3_consum); C_H2_ch_opt value(C_H2_ch); C_H2_dis_opt value(C_H2_dis); C_nh3_ch_opt value(C_nh3_ch); C_nh3_dis_opt value(C_nh3_dis); %% 目标函数各项成本 C_op_opt sum(c_grid .* P_grid_opt); C_co2_opt 0;我把优化结果存入.mat文件方便后续可视化和结果分析时反复加载不用每次重新求解。4. 结果分析与曲线绘制4.1 电力平衡结果图优化求解后第一件事就是画电力平衡图检查有没有违反直觉的结果。我用Matlab的area命令画堆叠面积图y轴是各电源出力x轴是时段。电力平衡图能一眼看出各时段电力来源构成。我在典型场景下得到的结果是凌晨负荷低谷时段0-8点风电出力较高而负荷较低电网购电价处于谷时这个时段模型选择将富余电力输入电解槽制氢。白天光伏大发时段10-15点光伏出力加上部分风电直接供电同时电解槽继续消纳光伏。晚高峰时段18-22点购电价处于峰值燃料电池和氨燃料机组开始启动发电来替代高价电网电力。4.2 储氢储氨动态曲线储氢罐和储氨罐的状态曲线更能反映调度策略的特征。储氢罐的SOC曲线呈现谷进峰出的形态夜间富余风电转化为氢气存入储氢罐白天和傍晚由储氢罐供氢给燃料电池发电。这是典型的“时间搬移”操作本质上是利用氢储能实现电力在小时级的时间尺度上的套利。但是如果只看储氢罐看不出氨储能的特殊价值。储氨罐的SOC曲线就很有意思在制氢成本低的时段一部分氢气进入合成氨装置转化为氨储存起来而不是全部存进储氢罐。这背后的逻辑是储氨罐的容量上限远大于储氢罐在面临强风电、低负荷的极端弃风场景时多余的电力通过“电解-合成氨”被“固化”成了氨。用数据说话某天的模拟结果显示在25%的时间段内合成氨装置处于开启状态将富余氢气转化为液氨储存。到了晚高峰氨燃料机组启动发电消耗的氨占当日产氨量的40%左右。这样只能在制氢时段使用的电解槽容量通过氨介质实现跨时段移峰填谷。4.3 调度结果的经济性分析我对比了两组结果一组是含氢氨系统的完整优化调度另一组是去掉氨回路只保留“电-氢-电”的简化系统。完整系统的日运行成本含碳成本约比简化系统低12%-15%核心原因是氨回路提供了额外的储能容量和调节空间。在简化系统中储氢罐容量有限弃风时段多余电力无法消纳只能弃掉而加上氨回路后富余氢气被转化为氨弃风率显著下降相当于用合成氨的低效率换取新能源的高利用率。但这里有个前提合成氨装置的启动成本不能太高。我将合成氨装置的启停惩罚成本设为80元/次如果这个值过高模型宁愿弃风也不会启动合成氨回路。这个“设备启停成本对调度结果的影响”可以做灵敏度分析我后面会有专题章节讲。4.4 结果可视化技巧有读者问我结果图是怎么画的这里给出几个关键代码片段。曲线图设置% plot_results.m %% 电力平衡图 figure(Color, w, Position, [100 100 1200 500]); time 1:24; h area(time, [P_wt_opt, P_pv_opt, P_fc_opt, P_nh3_gen_opt, P_grid_opt]); set(h, LineWidth, 1.5); legend({风电出力,光伏出力,燃料电池出力,氨燃料机组出力,电网购电}, ... Location, northwest); xlabel(时段 (h)); ylabel(功率 (kW)); grid on; box on;堆叠面积图的顺序会影响图层的可视效果我的经验是将最大的电源放最底层这样上面较小的电源不会被遮挡。储能SOC曲线图可以用双y轴因为储氢罐的容量单位是kg储氨罐的容量单位也是kg但量级不同figure(Color, w, Position, [100 100 1200 400]); yyaxis left; plot(time, SOC_H2_opt(1:24), -o, LineWidth, 2); ylabel(储氢量 (kg)); yyaxis right; plot(time, SOC_nh3_opt(1:24), -s, LineWidth, 2); ylabel(储氨量 (kg)); xlabel(时段 (h)); legend({储氢罐SOC, 储氨罐SOC}, Location, best); grid on;5. 常见问题与调试经验5.1 求解器报“Infeasible Problem”怎么办这个是我被问最多的问题也是初版模型的常态。当模型不可行时不要急着改约束先做三个检查检查参数一致性电解槽最大功率能否覆盖负荷峰值储氢罐的初始容量是否介于最小容量和最大容量之间我遇到过把储氢罐初值设成60而最大容量只有50的情况模型自然无解。检查末状态约束加了循环约束后初始和末状态必须相等这等于要求整个调度周期内储氢罐的总充入量等于总放出量。如果储能效率不为1充入后会有损耗模型必须通过电解槽的“超额产氢”来维持平衡。如果电解槽的最大产氢能力不足以弥补储能损耗就会出现不可行。解决办法是检查效率参数的合理性。逐步放开约束用optimize(Constraints, Objective)求解前先跑一次只包含等式平衡约束的线性规划去掉不等式看是否可行。如果线性松弛模型都不可行说明问题出在平衡约束本身比如功率单位不一致导致的数量级错误需要检查数据单位。如果线性松弛可行但MILP不可行问题出在0-1变量耦合的不等式约束上。5.2 YALMIP常见语法坑YALMIP的约束语句中和都可以使用但两个方向混用时容易写反。特别注意等式约束和不等式约束同时存在时Constraints [Constraints, ...]的拼接顺序不影响求解但如果你想快速查看某个约束有没有被成功加入可以用Constraints(end)查看。另一个高频错误是sdpvar的维度定义。如果定义为sdpvar(1, T)而在后面的索引中写成了SOC_H2(t)而不是SOC_H2(t1)会出现维度不匹配的报错。我建议在定义变量时加上注释标明维度含义避免混用。5.3 求解效率优化技巧当调度周期从24小时扩展到168小时一周时模型规模会增长7倍求解时间会从秒级飙升到分钟级甚至更久。有几个技巧可以提速删除冗余约束如果某些设备在特定时段不可能运行比如夜间光伏出力为0光伏相关的启动状态约束可以提前固定为0可以直接写死而不引入变量。减少0-1变量储能设备的充放状态变量可以通过耦合约束间接控制比如强制$C_{ch}(t) \le M \cdot C_{dis}(t)$可以避免同时充放但这种方法只对“连续容量”型储能有效具体要看设备特性。设置求解器参数cplex的mip.strategy参数调整为1深度优先搜索可以在某些模型上提速mip.limits.nodes设置节点数上限防止卡死。5.4 参数灵敏度的发现我对氨回路的几个关键参数做了灵敏度分析发现系统对合成氨转化效率的敏感程度超出预期。将$\eta_{NH3}$从0.85降到0.65相当于使用了较落后的合成氨设备氨回路在最优解中的使用率显著下降储氨罐的SOC长期处于低位系统几乎退化成纯氢储能系统。这说明合成氨回路的经济竞争力高度依赖于转化效率这个指标。另一个敏感参数是氨燃料机组的发电效率。$\eta_{NH3,gen}$从0.4降到0.3时氨燃料机组的出力显著减少用户在这个时段改用电网购电。通过敏感性分析得到的结论是与其追求大容量储氨罐不如先把氨燃料机组的效率做好每提升1个百分点效率系统日运行成本约下降2%-3%。5.5 画图时的一个小坑我在堆叠面积图中发现如果某一时段所有电源出力都为0比如停电检修的场景area函数的图层可能会出现不连续的情况。解决办法是在调用area之前用fillmissing对数据进行预处理或者将0替换为1e-6数量级的极小值避免图层断裂。6. 扩展方向与模型演进目前这个版本是单目标成本与碳加权的确定性优化模型没有考虑新能源出力和负荷的不确定性。下一阶段的扩展方向有几个引入场景法随机优化为风速、光照和负荷分别生成若干典型场景用随机规划框架建模目标函数变成各场景下成本的期望值。这样得到的调度策略对不确定性具有更强的鲁棒性。加入滚动时域控制将24小时静态优化改为MPC形式的滚动优化每个时间步更新预测数据并重新求解。虽然需要在线求解MILP但对于小时级调度来说计算时间完全充裕。考虑设备寿命损耗电解槽和燃料电池的启停次数直接影响设备寿命可以在目标函数中加入启停次数约束或寿命折算成本实现运行经济性与设备寿命的联合优化。氨能的多场景应用扩展除了发电氨还可以直接作为燃料供应工业锅炉、作为交通燃料、甚至作为化工原料外售。调度模型的目标函数中加入氨的外售收益项系统可以灵活选择“储氨发电”还是“储氨外售”策略空间会更大。我个人认为氢氨综合能源系统是个技术栈非常深的领域优化调度只是其中一环。做这个项目最大的体会是数学模型要逐步搭建代码调试要耐心细致每个约束条件的物理意义都要搞清楚。最后分享一个小技巧当你改了某个参数优化结果大变时先不要怀疑求解器出了Bug先用一个极端场景验证模型行为是否合理。比如将电解槽容量设为0系统应该退化为纯风光电网购电结构跑一下看结果是否符合物理直觉。如果符合模型基本可信如果不符合说明约束有隐性错误。这个验证习惯能帮你节省大量debug时间。