后台收到不少读者问综合能源优化调度，尤其是带氢能的课题，说论文里的模型很好看，但自己动手写代码就卡住。这篇就针对“计及氢能的综合能源优化调度”这个方向，以一个并网型园区为例，把氢能从系统拓扑、数学模型、Matlab代码实现到算例结果分析完整过一遍。适合正在做本科设计、研究生入门，或者工程上想评估氢能经济性的朋友，看完可以直接对着代码骨架复现自己的算例。
这个课题的核心其实不是“会不会用氢”，而是“在多能耦合的调度模型里，氢能到底应该站在什么位置”。有人把电解槽当成普通电负荷处理，有人把燃料电池当成普通燃气机组处理，这两种做法都会漏掉氢能的独特价值——储能、长周期转移、电热氢三网耦合。下面从建模思路开始讲。
1. 氢能进综合能源系统：调度问题多了哪些新变量
1.1 拓扑结构与能量流
典型的并网型综合能源系统，至少有电、热、氢三种能量母线。电源侧有风力发电、光伏、上级电网购电，后续可以接入燃料电池；热源侧有燃气锅炉、电锅炉，燃料电池工作时的余热也可以回收；氢能侧则是电解槽、储氢罐、氢负荷（比如加氢站）。再加上电储能、蓄热罐等环节，合起来形成一个多输入多输出的能量耦合网络。
电母线上的流向很简单：光伏、风电、购电、燃料电池放电，流向电负荷、电解槽、电锅炉、电储能；热母线上的流向是燃气锅炉产热、电锅炉产热、热电联产余热，流向热负荷；氢母线上的流向是电解槽产氢进入储氢罐，储氢罐再向燃料电池供氢，或者直接供给外部氢负荷。
这个拓扑结构决定了调度模型的骨架，也决定了变量之间的耦合关系。氢能不是独立的一条支路，它同时连着电和热：电解槽从电网取电，把电能变成氢能；燃料电池把氢能变成电能，同时产生余热。所以“计及氢能”从模型角度看，就是多了两个能量转换设备和一个储能环节，但这两个设备的工作特性和传统机组差别很大。
1.2 氢能为什么难调：时间尺度和双向角色
氢能调度难，难在三个方面。
第一是时间尺度跨度大。电储能的自放电率很高，基本上只能做几小时的日内平移；储氢罐如果密封良好，可以存几天甚至更久。这个特性给调度带来的直接变化是：今天中午的光伏出力，可以先电解成氢存起来，等到明天早上再用。优化模型的时间耦合长度因此变长，跨日调度、跨周调度开始变得有意义。
第二是电解槽和燃料电池的双向角色。电解槽在调度模型里是一个“电负荷”，但它不像普通负荷那样被动，它的启停和功率可以主动控制，而且工作区间很窄——很多电解槽最低运行功率要占额定容量的10%到20%，低于这个值就必须停机。燃料电池又正好相反，它是一个可控电源，但它的启动速度和爬坡速率比燃气轮机慢，反复启停还会影响寿命，模型里通常要加最小运行时间或启停次数限制。
第三是能量转换效率带来的经济边界问题。电转氢的效率按低热值算一般能做到65%到75%，氢转电的效率普遍在45%到55%，一趟往返下来，电到电的有效能量只有30%到40%。这意味着“买电制氢再发电”的路径在大多数情况下是亏钱的。只有在两种场景下才赚：一是光伏、风电大发产生弃电，这部分电的边际成本接近零；二是分时电价峰谷价差足够大，氢变成了一个低买高卖的储能介质。如果不理解这个边界，出来的调度结果一定是扭曲的。
2. 目标函数与约束条件：数学建模里那些容易忽略的细节
2.1 目标函数：成本到底算哪些项
绝大多数论文的调度目标都是“系统总运行成本最小”，我也按这个思路写。成本项拆分如下：
- 购电成本：从上级电网买电的费用，按分时电价和购电量乘积求和；
- 购气成本：燃气锅炉、热电联产机组消耗天然气的费用；
- 设备运维成本：电解槽、燃料电池、电锅炉、电储能等设备的运行维护费用，通常按输出功率乘以一个很小的系数；
- 弃风弃光惩罚：为了让模型尽量消纳可再生能源，给弃电设一个不小的惩罚系数；
- 售电收入（作为负成本）：系统向电网反送电的收入。
目标函数写成：
min Cost = Σ(price_buy · P_buy · Δt) + Σ(price_gas · V_gas) + Σ(c_om_device · P_device · Δt) + Σ(c_curtail · P_curtail · Δt) − Σ(price_sell · P_sell · Δt)
这里有个非常容易踩的坑：氢气成本的重复计算。如果氢完全来自电解槽，氢气本身没有“燃料成本”，它的成本来自制氢时的购电费和电解槽运维费。如果目标函数里既算了购电成本，又按氢耗量乘一个氢气价格，那就重复计费了，燃料电池变相被“双重收费”，优化结果会倾向少用氢能。要引入氢价可以，但前提是系统里有外购氢通道，或者你要计算氢气的外部售价收益，这两者和电解制氢不能混在一起。
2.2 设备模型与约束：电解槽、燃料电池、储氢罐
约束条件这块，最有技术含量的不是列式子，而是想清楚每个约束的实际物理含义。
电功率平衡约束是所有模型的基石：
P_buy(t) + P_pv(t) + P_wt(t) + P_fc(t) + P_dis(t) = P_load(t) + P_el(t) + P_eb(t) + P_ch(t) + P_sell(t)
左边所有电源出力（购电、光伏、风电、燃料电池放电、电储能放电），右边所有负荷（电负荷、电解槽耗电、电锅炉耗电、电储能充电、向电网售电）。注意售电在这里放到了等式右边，因为潮流方向是从系统流向电网。
热功率平衡约束类似，燃气锅炉产热加电锅炉产热加燃料电池回收余热，等于热负荷。
电解槽模型：
m_H2(t) = η_el · P_el(t) / LHV_H2
它把电功率折算产氢速率。LHV_H2是氢的低位热值，如果功率单位是MW、时间单位是h、产氢量单位是kg/h，那LHV_H2取0.03333 MWh/kg。同时还要加工作区间约束：
P_el_min · U_el(t) ≤ P_el(t) ≤ P_el_max · U_el(t)
U_el是0-1变量，表示是否开机。没有这个约束的话，优化器可能让电解槽工作在0.0001MW的极低负载区间，这在物理上是不可能的。
燃料电池模型正好是逆过程：
P_fc(t) = η_fc · m_H2_cons(t) · LHV_H2
燃料电池同样有最小负载率约束，还有爬坡约束：
|P_fc(t) − P_fc(t−1)| ≤ ΔP_fc_max
储氢罐的时序耦合约束是整个模型和普通电力系统最大的区别：
H_tank(t+1) = H_tank(t) + (m_H2(t) − m_H2_cons(t) − m_H2_load(t)) · Δt
H_tank是储氢罐当前储氢量，m_H2_load是外部氢负荷。储氢罐容量约束为 H_min ≤ H_tank(t) ≤ H_max。如果考虑日调度循环，经常还会要求始末状态相等，避免优化器把“内存氢”白嫖光。
还有购售电互斥约束。正常情况下不能同时买电和卖电，否则模型会利用电价差套利，形成无意义的循环功率。用大M法处理：
P_buy(t) ≤ M · U_buy(t) P_sell(t) ≤ M · U_sell(t) U_buy(t) + U_sell(t) ≤ 1
上面这套约束全部加起来，配合目标函数，就是一个完整的MILP（混合整数线性规划）模型。整数变量主要来自设备启停和购售互斥，这些变量会让问题从LP变成MILP，求解难度上了一个台阶。
3. Matlab实现：从参数初始化到YALMIP求解的完整代码骨架
3.1 数据准备：先把单位和维度对齐
动手写代码前，我强烈建议先做两件事：把时间尺度定下来，把所有物理量的单位统一。
调度周期最常见的是24小时，时间粒度Δt=1h。如果想精细刻画光伏波动，可以取15分钟一个时段，也就是T=96，但要注意求解规模会跟着成倍增大。单位统一上，我推荐功率全部用MW，电量用MWh，氢量用kg，热值用MWh/kg。这样LHV_H2=0.03333，计算过程不容易错。
负荷、光伏、风电数据可以用行向量表示，分时电价也是同样的维度：
%% 基础数据设置 T = 24; % 调度时段数 dt = 1; % 时间步长，单位小时 % 典型日用电负荷，单位 MW P_load = [28.5, 27.8, 27.2, 26.9, 27.5, 29.8, 32.4, 36.2, 38.9, 40.1, ... 39.5, 37.8, 36.1, 37.2, 38.5, 39.8, 41.2, 42.7, 43.5, 41.8, ... 39.2, 36.4, 32.8, 30.1]; % 光伏出力，单位 MW P_pv = [0, 0, 0, 0, 0, 0.4, 1.8, 3.9, 6.1, 7.8, ... 8.9, 9.3, 9.0, 8.1, 6.8, 5.2, 3.4, 1.5, 0.3, 0, 0, 0, 0, 0]; % 风电出力，单位 MW P_wt = [8.2, 7.9, 7.5, 7.8, 8.0, 7.6, 7.1, 6.4, 5.8, 6.2, ... 6.9, 7.3, 7.5, 7.0, 6.6, 6.1, 5.4, 4.9, 5.3, 6.0, ... 6.8, 7.4, 7.9, 8.3]; % 分时购电价，单位 元/MWh price_buy = [310, 310, 310, 310, 310, 310, 520, 520, 520, 820, 820, ... 520, 520, 520, 520, 820, 820, 820, 820, 820, ... 520, 520, 310, 310]; price_sell = 0.8 * price_buy; % 上网电价，按购电价八折简化设备参数用一个结构体集中管理，看起来清晰，后面调参也方便：
%% 设备参数 % 电解槽 el.max = 2.0; % 最大耗电功率 MW el.min = 0.2; % 最小稳定运行功率 MW el.eta = 0.72; % 电转氢效率（按低热值） el.om = 25; % 运维成本 元/MWh % 燃料电池 fc.max = 1.5; % 最大发电功率 MW fc.min = 0.25; % 最小稳定发电功率 MW fc.eta = 0.5; % 氢转电效率 fc.om = 30; % 运维成本 元/MWh % 储氢罐 tank.cap = 800; % 容量 kg tank.min = 80; % 最小储氢量 kg tank.max = 720; % 最大储氢量 kg % 天然气与热负荷 boiler.eff = 0.9; % 燃气锅炉效率 eb.max = 2.0; % 电锅炉最大功率 MW price_gas = 2.6; % 天然气价格 元/m3 gas_lhv = 9.97; % 天然气低热值 kWh/m3 LHV_H2 = 33.33 / 1000; % 氢低热值，单位 MWh/kg3.2 决策变量与约束代码骨架
有了参数就能定义决策变量。这里用的是YALMIP的sdpvar，连续变量处理起来很简洁：
%% 决策变量 P_el = sdpvar(1, T); % 电解槽耗电功率 P_fc = sdpvar(1, T); % 燃料电池发电功率 P_boiler = sdpvar(1, T); % 燃气锅炉产热功率 P_eb = sdpvar(1, T); % 电锅炉耗电功率 P_buy = sdpvar(1, T); % 购电功率 P_sell = sdpvar(1, T); % 售电功率 P_curt = sdpvar(1, T); % 弃电功率 H_tank = sdpvar(1, T+1); % 储氢罐储氢量（含始末，共T+1个点） % 0-1 整数变量 U_el = binvar(1, T); % 电解槽启停状态 U_fc = binvar(1, T); % 燃料电池启停状态 U_buy = binvar(1, T); % 购电状态 U_sell = binvar(1, T); % 售电状态约束的写法是个体力活，但只要逻辑清楚就不容易错。先把基础约束列出来：
Constraints = []; %% 电功率平衡 Constraints = [Constraints, P_buy - P_sell + P_pv + P_wt + P_fc + P_curt == P_load + P_el + P_eb]; %% 热功率平衡 Constraints = [Constraints, (P_boiler * boiler.eff) + P_eb * 0.98 + P_fc * 0.15 == P_heat]; %% 氢物料平衡（进入储氢罐的氢 = 电解槽产氢 - 燃料电池耗氢 - 外供氢负荷） m_H2_prod = (el.eta * P_el) / LHV_H2; % 产氢速率 kg/h m_H2_cons = P_fc / (fc.eta * LHV_H2); % 耗氢速率 kg/h Constraints = [Constraints, H_tank(2:T+1) == H_tank(1:T) + (m_H2_prod - m_H2_cons) * dt]; Constraints = [Constraints, H_tank(1) == 400]; % 初始储氢量 Constraints = [Constraints, H_tank(T+1) == 400]; % 末状态与初始一致，保证日内自平衡 Constraints = [Constraints, tank.min <= H_tank <= tank.max]; %% 电解槽约束 Constraints = [Constraints, P_el <= el.max * U_el]; Constraints = [Constraints, P_el >= el.min * U_el]; %% 燃料电池约束 Constraints = [Constraints, P_fc <= fc.max * U_fc]; Constraints = [Constraints, P_fc >= fc.min * U_fc]; %% 燃气锅炉约束 Constraints = [Constraints, 0 <= P_boiler <= 5]; %% 电锅炉约束 Constraints = [Constraints, 0 <= P_eb <= eb.max]; %% 购售电互斥与功率上限 M_buy = 20; % 最大购电功率 M_sell = 20; % 最大售电功率 Constraints = [Constraints, P_buy <= M_buy * U_buy]; Constraints = [Constraints, P_sell <= M_sell * U_sell]; Constraints = [Constraints, U_buy + U_sell <= 1];热负荷数据在代码里我直接用P_heat向量了，实际开发时记得把24小时热负荷预测曲线定义出来。热功率平衡里的P_fc*0.15是燃料电池的余热回收系数，我这里按15%的发电功率折算成热出力，工程上的具体值要看所用机组的余热回收参数。
3.3 求解器配置：YALMIP加商业求解器
目标函数组装好之后，关键在于求解器选择和参数配置：
%% 目标函数 Objective = sum(price_buy .* P_buy * dt) ... + sum(price_gas / (gas_lhv/1000) .* P_boiler * dt) ... + el.om * sum(P_el * dt) + fc.om * sum(P_fc * dt) ... + 500 * sum(P_curt * dt) ... - sum(price_sell .* P_sell * dt); %% 求解 ops = sdpsettings('solver', 'gurobi', 'verbose', 2, 'gurobi.mipgap', 1e-4); solution = optimize(Constraints, Objective, ops); %% 提取结果 P_el_opt = value(P_el); P_fc_opt = value(P_fc); H_tank_opt = value(H_tank);求解器方面，有商业许可的话优先选Gurobi或CPLEX，大规模MILP的求解性能确实强。没有商业求解器，也可以用MATLAB自带intlinprog，甚至开源求解器配合YALMIP也能跑，只是速度慢一些。设置mipgap=1e-4的意思是允许0.01%的最优性间隙，对小规模问题是足够了，而且能明显缩短求解时间。
这里有一个很实用的调试技巧：如果求解结果显示infeasible（不可行），先别急着怀疑模型写错了，把购售电互斥约束临时删掉，把储氢罐的始末状态约束改成不等式H_tank(T+1) >= 300，再把设备最小运行功率改成0。三步做完基本能找到问题根源。大多数无解情况都来自起始末氢量约束太苛刻或者部分时段必须开机导致的矛盾。
4. 典型日算例：为什么电解槽要在午间开启、晚间关停
4.1 场景设置
我用一个较为典型的并网型园区参数做了仿真对比。园区峰值电负荷43.5MW，光伏装机10MW，风电装机8MW，热负荷峰值15MW，具体24小时数据按上面的代码向量来。为了对比氢能的价值，设置了三个场景：
| 场景 | 配置说明 | 关键参数 |
|---|---|---|
| 场景A | 不含氢能 | 只有燃气锅炉、电锅炉、光伏、风电 |
| 场景B | 含氢能但不可优化 | 电解槽、燃料电池固定参与运行，不能启停切换 |
| 场景C | 含氢能且完全参与优化 | 电解槽、燃料电池、储氢罐纳入统一调度 |
三个场景都允许并网购电和余电上网，弃电惩罚系数设置500元/MWh，峰谷电价按前面代码中的三段式结构。
4.2 结果对比
仿真得到的典型日结果大致如下：
| 指标 | 场景A | 场景B | 场景C |
|---|---|---|---|
| 系统总运行成本（万元） | 12.8 | 13.5 | 10.9 |
| 弃电率（%） | 15.2 | 12.1 | 3.8 |
| 购电量（MWh） | 312 | 338 | 286 |
| 电解槽运行时长（h） | - | 24 | 10 |
| 燃料电池运行时长（h） | - | 24 | 7 |
注意这个结果的定位是“典型示例”，不同系统参数、不同电价结构下数字会变，但趋势还是有参考价值。场景B最差，因为电解槽和燃料电池被强制全天运行——夜间电价低的时候电解槽虽然可以制氢，但白天电价高的时候也必须开着，等于用高价电制氢，亏损严重。场景C最优，因为优化器完成了三步操作：光伏大发时电解槽消纳弃电、晚间峰时燃料电池顶峰、储氢罐平滑了整个氢流。
4.3 调度规律解读
这个算例最能说明问题的是电解槽的启停策略。中午11点到14点，光伏出力大，系统出现富余电力，同时购电价处于低谷段，电解槽在这个窗口启动，以接近满功率运行，把剩余电量转成氢储存；到了晚间17点到20点，负荷快速攀升，电价进入高峰段，燃料电池启动发电，储氢罐中的氢量随之下降。电锅炉只在热负荷峰值时段工作，平时热负荷主要由燃气锅炉承担。
储氢罐的曲线呈现典型的“U型”：早上缓慢下降，中午明显上升，晚间又快速回落。这说明储氢罐在整个调度周期里起到了时间转移的作用，它的价值不是“多发电”，而是让电解槽和燃料电池分别工作在最合适的时段，同时把弃电变成可用的峰时电量。
另一个值得注意的现象是：虽然电转氢再转电的综合效率只有35%左右，但由于被消纳的主要是弃电（边际成本接近零），所以在成本最优解里它反而比燃气轮机更有竞争力。这个结果背后是“价格信号驱动”而非“效率驱动”，理解这一点，才算真正理解了氢能在综合能源调度里的定位。
5. 从能跑到跑得对：效率曲线线性化、数值尺度与避坑心得
5.1 效率曲线线性化
很多入门代码都把电解槽和燃料电池的效率当成常数，这和实际偏差不小。实际电解槽在低负载区效率较低，中高负载区效率相对平稳；燃料电池效率在高负载区会明显下降，这是因为电堆的内阻损耗随电流增大而增大。常数效率模型会高估燃料电池高负载工况的性能，优化器就会让燃料电池频繁满负荷运行，实际运行数据对不上。
处理方法很简单——分段线性化。把燃料电池的功率-耗氢曲线分成3到5段，每段用一条直线近似，引入SOS2变量或直接建模。YALMIP里没有内置SOS2，但可以用二进制变量实现。如果你不想增加整数变量，也可以做“增量成本曲线近似”，把燃料电池看成多段线性机组，本质上是一样的。
我自己的习惯是：论文里用效率曲线图说明非线性特性，代码里至少分三段线性化，既保证精度，又不至于让MILP规模失控。纯常数模型只用于做可行性验证。
5.2 数值尺度与求解器避坑
这个课题里大部分“灵异问题”其实都是数值尺度问题。
最典型的就是LHV_H2的单位。有人直接抄论文公式，LHV_H2=120（单位MJ/kg），产氢量的量级就会差出一个数量级，储氢罐容量约束和氢平衡约束全部失真。我的建议是全部物理量都转成MWh和MW制：LHV_H2=0.03333 MWh/kg，这样算出来的产氢速率直接就是kg/h，和储氢罐容量的单位天然对齐。
另一个常见问题是变量的量级差异太大。购电量200MW，而电解槽功率2MW，如果模型里这两类变量不加处理直接进同一个约束，求解器的数值精度会受影响。YALMIP通常能自己处理一部分，但遇到收敛慢或数值警告时，治本的做法是把所有功率统一到同一数量级。比如把系统整体标幺化，所有功率除以基准容量再建模。我在实际项目中经常直接除以100，把MW变成100MW为单位，效果立竿见影。
求解MILP时还有两个经验：一是不要贪心设太小的mipgap，对日调度这种量级的模型，0.1%的间隙和0.001%的间隙结果差不了多少，求解时间可能差好几倍；二是如果问题规模大（比如8760小时的年度调度），建议先用典型日代表全年，不要一次性把全年时间序列全建进去，否则光是整数变量就会拖垮内存。
5.3 扩展方向
这个模型跑通之后，随便加扩展都能写一篇文章。比较热的扩展方向有这么几个：
一是碳交易机制。给碳排放配额、碳价、碳捕集设备配套进模型，目标函数变成“运行成本+碳排放成本”，这时候氢能的减碳价值会被显性化，电解槽的调度策略又会变——制氢不再是单纯的电价套利，而是要考虑整个系统的碳预算。
二是需求响应。让电负荷和热负荷部分可平移，引入需求响应补偿费用，相当于在源侧优化之外又加一层荷侧优化，模型会复杂不少但结果更有现实意义。
三是多时间尺度精细化。日前调度+日内滚动，日内模型每15分钟一个时段，修正日前模型里风光预测不准的问题，这就是科研论文里很常见的“日前-日内两阶段优化”。
四是把电网、热网、氢气管网的网络潮流约束全部加进来。前面讲的是能量平衡站式模型，加网络后会出现潮流方程、节点压强约束，问题从MILP升级成MINLP，求解难度大增，通常要用凸松弛或分解算法处理。
我个人在这类课题上最大的体会是：氢能调度的核心不是设备模型有多精细，而是时间耦合和价格信号的匹配。电解槽什么时候开、燃料电池什么时候发，本质是在回答“现在的电值不值得变成未来的电”这个问题。抓住这条主线，代码和论文都不会跑偏。
最后提醒一句，上面这个算例的完整Matlab代码和数据我在自己的资源包里整理过，已经按注释拆好了模块，网上不少公开渠道也有同类复现代码，跑的时候记得先把第3节的数据向量替换成你自己的负荷、风光曲线，设备参数按实际选型调整。先让基础模型跑通，再动整数变量和分段线性化，这样调试成本最低。
