做综合能源系统优化调度，最怕的是什么？不是模型复杂，而是你辛辛苦苦建完模、跑完代码，结果却不符合物理直觉——明明有大量弃风，调度结果却在用电高峰期买高价电。我第一次在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(t+1) = 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(t+1) - 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, T+1); % 储氢罐存储量 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:T+1) == 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(T+1) == 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 求解器配置与求解策略：Yalmip+Gurobi最常见的黄金组合
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指定求解器，verbose=2可以输出详细求解日志，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在所难免，但沿着物理直觉去排查，一定比对着报错信息瞎猜要快得多。
最后分享一个小操作经验：在你终于把模型调通、结果也合理之后，把边界条件（比如储氢罐容量、电解槽效率）往两个极端各改一次再跑一遍。如果结果变化符合物理预期，说明模型稳健；如果结果突然乱跳，说明那组参数附近可能存在数值敏感区，写论文时最好避开，或者补充说明。这个习惯能帮你省下跟审稿人来回扯皮的大量时间。
