企业官网建设流程全解析

我刚做完一个含氢气氨气的综合能源系统优化调度项目，用的Matlab代码实现，跑了各种工况，把这套系统的调度逻辑和代码实现细节整理一下。

先说结论：氢氨综合能源系统不是简单把电解槽、储氢罐、合成氨装置堆在一起，真正的核心难点在“多时间尺度耦合”和“跨介质能量转换”的建模与调度。纯电、纯气系统做优化调度相对成熟，但加入氢和氨之后，能量存储的时长尺度、转换环节的效率损失、以及化工过程的安全约束都会显著改变调度策略。

我这次在Matlab里实现了完整的优化调度模型，求解器用的YALMIP调用cplex，目标函数包含运行成本、碳排放成本和设备启停惩罚，约束条件覆盖电、氢、氨三种介质的平衡约束和设备运行约束。下面把整个项目的建模思路、代码实现和踩坑经验完整分享出来。

1. 完整系统架构与能量流分析

1.1 氢氨综合能源系统的基本组成

这套系统包含的核心设备有：风电机组、光伏阵列、电解槽、储氢罐、燃料电池、合成氨装置、储氨罐、氨燃料机组、以及常规电负荷和热负荷。

能量流动路径有三条主线：

• 电力线：风电和光伏直接供电，富余电力驱动电解槽制氢，燃料电池和氨燃料机组在电力不足时反向补电。
• 氢能线：电解槽产氢后进入储氢罐，一部分氢气直接供给燃料电池发电，另一部分进入合成氨装置作为原料。
• 氨能线：合成氨装置产出液氨存入储氨罐，氨燃料机组根据调度指令将氨转化为电力，实现跨季节储能。

这三条能量流在时间尺度上完全不对称。电的传输是瞬时的，储能只有几个小时到几天；氢的存储周期可以达到数周；而氨作为化学储氢介质，存储周期可以达到数月甚至跨季度。这就导致调度模型决策变量的时间窗需要拉得很长，否则无法体现氨储能的战略价值。

1.2 “电-氢-氨”多级耦合的调度逻辑

传统综合能源系统的优化调度通常只考虑电力平衡和热力平衡，本质上是单一时间断面内的资源配置问题。但氢氨系统引入了速率型设备（电解槽、燃料电池）和容量型设备（储氢罐、储氨罐），设备的输入输出关系变得复杂。

我在这套模型中把调度逻辑分为三层：

第一层是电力平衡层，决定每个调度时段风电、光伏、燃料电池、氨燃料机组和电网购电的出力组合，满足电负荷需求。这一层约束和传统微电网类似但多了氨燃料机组这个变量。

第二层是氢气平衡层，电解槽的产氢量、储氢罐的充放氢速率、燃料电池的耗氢量、合成氨装置的用氢量要满足节点平衡。储氢罐状态变量采用离散时间递推的形式，即t+1时刻的储氢量等于t时刻储氢量加上充放氢差。

第三层是氨平衡层，合成氨装置的产氨量进入储氨罐，氨燃料机组的耗氨量从储氨罐取出，储氨罐同样有容量约束和速率约束。

这三层平衡之间通过设备转换效率相互耦合。电解槽消耗电力产生氢气，燃料电池消耗氢气产生电力，合成氨装置消耗氢气产生氨，氨燃料机组消耗氨产生电力，形成了一种“电-氢-氨-电”的闭环转换链。每次转换都有能量损失，因此调度模型的核心是在时间维度上找到最优的转换时机和存储策略。

1.3 为什么用氨作为储能介质

项目初期我考虑过只做“风电+电解槽+储氢+燃料电池”的简单结构，但很快就发现一个问题：纯氢储能在长周期尺度上没有优势。高压储氢罐的自放电率虽然低，但设备投资和维护成本都很高，而且氢气的体积能量密度太低，即使压缩到70MPa也只有汽油的几分之一。

氨作为储氢介质有三大优势：

• 体积能量密度高。液氨的含氢量约为17.6wt%，单位体积含氢量比高压储氢罐高出数倍，适合大规模长时间存储。
• 存储条件温和。液氨在-33°C或0.8MPa下即可液化，比氢气的低温高压存储条件宽松得多，储罐成本低一个量级。
• 产业链成熟。合成氨是百万吨级工业，运输、储存、使用的标准体系完善，可以作为氢能的“载体”利用现有基础设施。

代价就是转换效率偏低。电解槽效率约为70%左右，合成氨的综合能耗（哈伯-博世法）会额外消耗大量热能，氨燃料机组的发电效率在40%左右，整个“电-氢-氨-电”循环的综合效率粗略估算在30%-35%之间。

这意味着氨储能路线在经济上肯定打不过锂电池（综合效率90%以上）做短时储能，但在周级以上甚至跨季度的长时储能场景下，氨的存储成本优势才会显现。调度模型必须同时核算小时级成本和月级成本，才能看清氨储能的真正价值。

2. 优化调度数学模型构建

2.1 决策变量定义

我先定义模型中的关键集合和决策变量，这些直接对应Matlab代码里的变量数组。

调度周期为24小时，时间分辨率1小时，即$T=24$个时段。

决策变量分为三类：

连续型变量：

• 风电光伏实际出力 $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_{t=1}^{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_{t=1}^{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_{t=2}^{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}(t+1) = 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, T+1); % 储氢罐状态 SOC_nh3_ = sdpvar(1, T+1); % 储氨罐状态 %% 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定义连续变量。注意我设置了储氢罐状态变量维度为$T+1$，是因为递推关系需要初值和终值两个端点。

定义目标函数代码如下：

%% 目标函数 % 分时电价参数 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(t+1) == SOC_H2(t) + C_H2_ch(t) - C_H2_dis(t)]; Constraints = [Constraints, ... params.SOC_H2_min <= SOC_H2(t+1) <= 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(T+1) == 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(t+1)，会出现维度不匹配的报错。我建议在定义变量时加上注释标明维度含义，避免混用。

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时间。