风光互补制氢合成氨系统容量-调度联合优化模型Matlab+Cplex实现

前阵子帮一个做绿氨项目的朋友复现了一套 并/离网风光互补制氢合成氨系统的容量-调度优化模型，Matlab 建模，调用 Cplex 求解。这套模型把风电、光伏的装机容量、电解槽规模、储氢罐大小、合成氨设备能力，跟全年逐时的运行调度一次性联合优化出来，目标函数是年化总费用最小。这类题目在近两年的绿氢绿氨项目前期论证里特别常见，网上同名标题不少，但真正能跑通、能解释清楚结果的不多。这篇就把我复现过程中从系统建模、数学公式推导、Cplex配置到结果解读的完整路径整理出来，适合正在做氢能方向优化、或者想把“容量-调度联合优化”写成可执行代码的读者。

在动手之前先明确一点：这本质上是一个优化问题，不是仿真问题。仿真是给定一个系统看它怎么跑，优化是在一堆可能的系统配置和运行策略里找出最便宜的那一套。我拿到需求后第一件事不是打开 Matlab，而是先把这个组合系统的物理逻辑盘清楚，再考虑怎么把它写成 Cplex 能啃的数学模型。

1. 项目背景与系统结构拆解

1.1 为什么是“风光互补制氢合成氨”

制氢合成氨的能量流其实很直白：风电和光伏发电，电进电解槽制氢，氢气再和氮气在合成氨装置里反应生成氨。反应式是最经典的 Haber-Bosch 工艺，N₂ + 3H₂ → 2NH₃，理论上每吨氨要消耗大约 0.176 吨氢，工程计算里常取 0.18 吨氢/吨氨。氮气来源一般不考虑，默认有空分装置供应，空分耗电可以按辅助负荷计入。

选择“风光互补”而不是单一可再生能源，本质上是在利用资源和负荷的时间互补性。光伏白天出力高但夜间为零，风电在夜间和冬季往往更充沛，两者叠加可以拉平全天候的出力曲线。这个特性对制氢系统来说是致命的：电解槽虽然可以频繁启停，但频繁启停影响寿命，且低负载率下效率明显下降。如果只用光伏，夜间大概率得停机，设备利用率低；如果只用风电，白天午高峰可能正好没风。风光互补能让电解槽在一天里维持更高的运行小时数，同样的制氢产量对应的电解槽装机可以更小，储氢罐容量也能降下来。

1.2 并网模式与离网模式的本质差异

并网和离网是这个项目里最重要的情景划分，两者的差别不是“有没有一根线连着电网”，而是系统的平衡机制完全不同。

并网模式下，电网就是系统的缓冲池。风光出力高了，多余的电可以卖给电网；风小光弱的时段，可以从电网买电维持合成氨装置连续生产。系统可靠性的压力被转嫁给了电网，所以容量配置可以相对“抠门”。

离网模式下，系统孤岛运行，所有电力必须来自本地风电光伏，所有缺额必须靠储氢罐或停机减产来消化。任何一个时段出现电力缺口，都可能引发制氢中断，进而影响合成氨的氢源供应。这种模式下系统必须有足够的冗余——更大的风电光伏装机、更大的储氢罐，甚至可能需要配套储能电池来平滑短时波动。

两种模式的对比可以用下面这张表概括：

这个对比在后面跑结果时会不断被验证，我复现时发现，离网模式下的最优容量往往要比并网模式大好几倍，这也是项目前期论证中最核心的经济权衡问题。

1.3 容量-调度为什么必须联合优化

这是这套模型的灵魂。传统的容量配置方法是给定风光资源曲线，估算一个“等效利用小时数”，然后反推装机容量和制氢设备规模。这种做法在单一电源加固定负荷的场景下勉强能用，但在风光互补、储氢罐、可调负荷耦合的场景里会严重失真。

原因在于：容量决策和调度决策本质上是互相嵌套的。风电光伏装多少，取决于电解槽能不能在风光出力高峰期把电吃掉；电解槽装多少，又取决于储氢罐能不能把白天多余的氢留到夜间用；储氢罐大小又反过来影响合成氨装置能不能全天候满负荷运行。这是一条闭合的因果链，任何一环单独定下来，其他环节的“最优值”都会变。

更准确地说，容量-调度联合优化是一个单模型两层的结构：容量变量是一次性决策，全年投入运行；调度变量是逐时决策，每个时段独立变化。如果把两层拆开求解，先定容量再跑调度，很可能得到局部最优——容量定小了，调度再怎么优化也补不回来。反过来，如果先假设一个调度策略再去定容量，又会因为调度策略本身不最优而高估容量需求。所以正确的做法是把容量变量和全时段的调度变量放在同一个优化模型里，让 Cplex 一次性求出联合最优解。

2. 数学模型构建思路

2.1 目标函数：年化总成本怎么拆

这个模型的优化方向是经济性最优，所以目标函数用“年化总费用最小”来定义。年化总费用包含三个部分：等年值投资成本、年运维成本、年运行成本（购售电净支出）。

投资成本必须先做等年值折算，因为风电、光伏这类设备设计寿命通常在20年以上，而电解槽寿命可能只有10年，合成氨装置又是另一套寿命。直接把总投资加在一起比较是不公平的，要用资金回收系数 CRF 把一次性投资摊到每年：

CRF = i × (1 + i)ⁿ / [(1 + i)ⁿ - 1]

其中 i 是折现率，n 是设备寿命。折现率一般取 6%～8%，风电光伏寿命按 20 年算，电解槽按 10 年，合成氨装置按 20 年。每类设备的投资乘上对应的 CRF，再加起来，就是每年要分摊的设备成本。这也是为什么同样一笔投资，电解槽看起来最贵，因为它的寿命最短，每年分摊比例更高。

目标函数的整体形式可以写成：

min f = CRF_w × c_w × P_w_cap + CRF_pv × c_pv × P_pv_cap + CRF_el × c_el × P_el_cap + CRF_h2 × c_h2 × S_h2_cap + CRF_nh3 × c_nh3 × P_nh3_cap + Σ(运维成本) + 购电成本 - 售电收益

这里我把每个设备的容量都单独列出来，分别是风电装机、光伏装机、电解槽额定功率、储氢罐容量、合成氨设备产能。每个容量乘上对应造价和 CRF，就是这项设备的年化投资。

2.2 变量设计与线性化处理

模型的变量分成三个层次：容量变量、调度变量、辅助的 0-1 变量。

容量变量是全局标量，决定系统规模。调度变量是每个时段一个值，比如各时段的风电消纳功率、光伏消纳功率、电解槽输入功率、储氢罐充放氢量、合成氨装置产氨量。0-1 变量用于刻画设备启停和购售电互斥逻辑，这也是模型变成混合整数线性规划（MILP）的根本原因。

变量清单如下：

线性化是建模的关键一步。原问题里如果采用变效率电解槽模型，制氢量和电功率之间会形成非线性关系，Cplex 处理起来非常吃力。我复现时做了工程简化：电解槽按额定效率恒定处理，制氢量等于电功率乘以效率再除以氢气低热值。合成氨装置的耗氢量按比例系数折算，也是线性关系。储氢罐的动态用容量平衡方程描述，同样是线性差分方程。

这样处理之后，整个模型就是一个标准的 MILP，可以放心交给 Cplex。实测下来，线性化带来的误差通常在工程可接受范围内，因为电解槽在 50%～100% 负载率区间的效率变化并不剧烈。

2.3 约束条件清单

约束条件是这个模型里最考验对系统理解的部分。逐条列一下我在复现时保证都在模型里的约束。

第一是电力平衡约束。每个时段，风电消纳、光伏消纳、购电功率之和，必须等于电解槽用电加辅助用电加售电功率。辅助用电主要是空分装置和合成氨装置的厂用电，我按产氨量乘一个固定系数计算。离网模式下，购电和售电两项强制为 0，电力平衡变成纯粹的本地功率平衡。

第二是风电光伏出力约束。每个时段的风电消纳功率不能超过该时段可利用风功率，光伏同理。可利用风功率由风速-出力曲线换算得到，光伏则由辐照度换算。消纳功率一定小于等于理论出力，多出来的部分可以弃掉，弃电不会产生成本但会浪费，所以模型会自动权衡“多装多吃”和“装多了用不完”的经济账。

第三是电解槽运行约束。电解槽输入功率有上下限，下限取额定功率的一定比例，比如 20% 或 30%，低于这个值必须停机。这一约束对应到二进制变量 u_el(t)。

第四是制氢量等式。制氢量等于电功率乘以效率除以氢气低热值，这个等式的物理含义是能量从电能转化为氢能，效率取 0.7 左右。

第五是储氢罐动态。储氢罐的氢量用差分方程描述：下一时段氢量等于当前氢量加产氢量减耗氢量。容量有下限和上限，下限保证合成氨装置不至于断供，上限就是储氢罐容量。这一约束是链接制氢和合成氨的纽带，也是整个系统的时间耦合核心。

第六是合成氨产量约束。每个时段的产氨量等于该时段耗氢量除以 0.18，即每生产一吨氨消耗 0.18 吨氢。同时保证全年或典型日的氨产量达到目标值，这是一个全局累加约束，直接决定了系统的产能下限。

第七是并网交互约束。购电功率和售电功率不能同时为正，用二进制变量做互斥限制，同时各自的功率上限按并网协议设定。

把以上约束全部写进去后，模型就是完整的。实际复现时你会发现，最难调的不是哪个约束写错了，而是这些约束之间的量纲和数值尺度差异。这个坑我留到第五章专门讲。

3. Matlab+Cplex求解实现要点

3.1 开发环境准备与Cplex配置

网上关于“Cplex配置到Matlab”的提问很多，我这次也把配置过程踩了一遍。最顺手的方式是用 Yalmip 作为建模中间层，Yalmip 负责把数学模型翻译成 Cplex 能识别的标准形式，你不需要直接面对 Cplex 底层的 C API 或者 mll 文件。

配置步骤大致如下：先从 IBM 官网下载对应平台的 Cplex 安装包，区域求解器用免费社区版一般就够跑教学和小规模算例。安装完成后，Cplex 会在安装目录里带一套 Matlab 接口文件。打开 Matlab，把 Cplex 的matlab子目录添加到 MATLAB 路径里，同时把 Yalmip 的目录也加进去。然后跑一个测试命令：

yalmiptest

这个命令会把 Yalmip 支持的所有求解器全部探测一遍，如果输出结果里 Cplex 那一行显示 OK，说明接口已经通了。我遇到过的最常见问题是版本不匹配，Cplex 12.10 配特别老的 Matlab 版本会有 JVM 相关报错，处理办法是查看 Cplex 官方支持矩阵，选一个与本地 Matlab 版本兼容的 Cplex 版本。

3.2 Yalmip建模代码骨架

这里给一段最简的骨架代码，能跑通就能看到模型的基本结构。以 24 小时为例，实际复现时按典型日或者全年 8760 小时扩展，代码逻辑完全一致。

T = 24; % 时段数，典型日场景 % ----- 容量变量 ----- P_w_cap = sdpvar(1, 1); P_pv_cap = sdpvar(1, 1); P_el_cap = sdpvar(1, 1); S_h2_cap = sdpvar(1, 1); P_nh3_cap = sdpvar(1, 1); % ----- 调度变量 ----- p_w = sdpvar(1, T); % 风电消纳功率 p_pv = sdpvar(1, T); % 光伏消纳功率 p_buy = sdpvar(1, T); % 购电功率 p_sell = sdpvar(1, T); % 售电功率 p_el = sdpvar(1, T); % 电解槽输入功率 u_el = binvar(1, T); % 电解槽启停 s_h2 = sdpvar(1, T+1); % 储氢罐氢量 h_el = sdpvar(1, T); % 产氢量 h_nh3 = sdpvar(1, T); % 合成氨耗氢量 m_nh3 = sdpvar(1, T); % 产氨量 % ----- 参数 ----- eta_el = 0.7; % 电解效率 LHV_H2 = 120; % 氢气低热值 MJ/kg k_h2_nh3 = 0.18; % 吨氨耗氢量 t/t P_el_min_ratio = 0.2; % 电解槽最小负载率 Constraints = []; for t = 1:T % 电力平衡 Constraints = [Constraints, p_w(t) + p_pv(t) + p_buy(t) == ... p_el(t) + p_sell(t)]; % 风光出力上限 Constraints = [Constraints, 0 <= p_w(t) <= P_w_cap * avail_w(t)]; Constraints = [Constraints, 0 <= p_pv(t) <= P_pv_cap * avail_pv(t)]; % 电解槽功率上下限与启停 Constraints = [Constraints, P_el_min_ratio * P_el_cap * u_el(t) <= p_el(t)]; Constraints = [Constraints, p_el(t) <= P_el_cap * u_el(t)]; % 制氢量 Constraints = [Constraints, h_el(t) == eta_el * p_el(t) * 3600 / LHV_H2 / 1000]; % 储氢罐动态 Constraints = [Constraints, s_h2(t+1) == s_h2(t) + h_el(t) - h_nh3(t)]; Constraints = [Constraints, 0 <= s_h2(t) <= S_h2_cap]; % 合成氨产氢耗与产氨量 Constraints = [Constraints, m_nh3(t) == h_nh3(t) / k_h2_nh3 / 1000]; end % 目标函数：等年值投资 + 运维 + 购电 - 售电 Objective = CRF_w * c_w * P_w_cap + CRF_pv * c_pv * P_pv_cap + ... CRF_el * c_el * P_el_cap + CRF_h2 * c_h2 * S_h2_cap + ... CRF_nh3 * c_nh3 * P_nh3_cap + ... sum(om_w * P_w_cap + om_pv * P_pv_cap + om_el * P_el_cap) + ... sum(price_buy .* p_buy) - sum(price_sell .* p_sell); % 求解 ops = sdpsettings('solver', 'cplex', 'verbose', 2); ops.cplex.mip.tolerances.mipgap = 0.01; optimize(Constraints, Objective, ops);

注意代码里avail_w(t)和avail_pv(t)是 t 时段的风光可利用出力因子，由气象数据预处理得到，取决于风速、辐照度和对应机组的功率特性曲线。price_buy和price_sell是分时电价序列，并网模式用得到，离网模式这两个序列可以全置零。

3.3 大规模问题加速技巧

如果直接把全年 8760 小时全部放进模型，规模是很大的。每个时段十几个变量和约束，8760 个时段就是十万以上量级的约束和变量，再加上二进制变量，Cplex 分支树会非常庞大，从几分钟到几小时都有可能。

我复现时做的第一个加速处理是场景削减。风电光伏日曲线虽然每天都在变，但形态高度相似，比如夏季晴天、夏季阴天、冬季晴天、冬季阴天。用 K-means 聚类或者同步回代削减算法，把全年 8760 小时压缩成 8～24 个典型日，每个典型日加权表示该类型日出现的天数。这样模型规模缩小到原来的几十分之一，而成本估计误差通常能控制在 5% 以内。

第二个加速手段是给 Cplex 设置合理的收敛标准。工程测算不需要精确到小数点后五位，mipgap设置为 1%～2% 完全够用，求解时间可以缩短一个数量级。如果时间仍然紧张，可以再加cplex.mip.strategy.search = 2或者限制最大节点数。

第三个技巧是变量初始值。先把容量变量固定到一个合理的估计值，跑一次纯调度优化，把得到的储氢罐状态作为初始解喂给联合优化模型，Cplex 的热启动机制能显著减少分支下的搜索空间。这个操作在代码里只需要在optimize前把yalmip变量的初始值赋好即可。

4. 并离网结果对比与工程解读

4.1 参数与场景设置

我复现时用的是一组公开的西北地区气象数据，风速和辐照度按小时记录，覆盖一整年。电价按典型工业分时电价设置，峰时 0.9 元/kWh、平时 0.6 元/kWh、谷时 0.3 元/kWh，售电按上网电价 0.35 元/kWh 计。设备造价按当前市场水平：风电 4500 元/kW、光伏 3500 元/kW、电解槽 2000 元/kW、储氢罐按材料成本折算、合成氨装置按吨日产能计价。

特别说明一点：这些参数只是示例，不同地区、不同年份的造价差异很大。真正的项目论证应该用可研报告里的最新价格，这篇文章重点在于展示模型的逻辑和结果分析思路，绝对值不具备普适性。

4.2 并网模式的最优配置结果

在并网模式下，系统最优配置有明显的“电网依赖”特征。风电和光伏的装机不需要满足全部制氢负荷，因为缺电时段可以从电网买电。电解槽的规模主要取决于电价低谷时段能吸收多少电，同时考虑风光出力特性，储氢罐的容量只需要覆盖一天的缓冲即可。

我跑出来的典型结果呈现这种趋势：光伏配比要比纯离网模式高得多，因为白天光伏出力高峰和高制氢需求恰好重合，电价高的时段可以少买电或卖电给电网。风电装机相对保守，只需要在夜间和冬季补充出力。储能方面，储氢罐容量很小，因为合成氨装置可以通过购电维持运行，储氢罐只起短时调节作用。

并网模式的经济账也简单：购电成本是年度运行费用的主要组成部分，但设备投资省下来了。对项目方来说，只要能谈下比较便宜的绿电交易电价，并网模式的单位氨成本通常显著低于离网模式。

4.3 离网模式的最优配置结果

离网模式的容量配置会“豪放”很多。因为没有电网兜底，风电和光伏的装机必须覆盖最恶劣的气象场景，而不是平均场景。我复现时另一个很深的感受是：离网模式的结果对场景选取极度敏感。如果只用夏季典型日来削减场景，很可能得到一个在冬季连续无光的时段完全无法满足氨产量的配置。

在系统结构上，离网模式会出现几个明显特征：光伏配比下降，风电配比上升，因为风电在夜间能出力，对全天候制氢的贡献更稳定；储氢罐容量显著增大，往往需要满足合成氨装置 2～3 天的氢消耗量；电解槽额定功率与最大制氢能力挂钩，反复在满负荷和最低负荷之间切换的情况增多。

经济上，离网模式的单位氨成本会比并网模式高出不少，这是由两方面的原因造成的：一是装机容量增大导致投资成本上升，二是弃电率升高导致能量利用效率下降。风光装机翻倍之后，大量电力出现在低负荷时段，电解槽吃不下，只能弃掉。

4.4 对比结论如何指导工程决策

把两种模式放在同一张表里对比，结论就非常直观了：

工程决策的角度看，如果项目所在地电网接入条件好、电价平稳，并网模式是更经济的选择。如果项目定位是“零碳绿氨”出口，必须严格证明生产全过程的绿电属性，那就需要离网模式或者至少是“高比例自用电+极小比例外购”的准离网模式。

我个人的建议是：前期论证时可以同时跑并网和离网两个情景，把电价、设备造价都做敏感性分析。如果并网模式和离网模式的成本差距在 10% 以内，说明这个项目的核心成本驱动因素是设备投资和风光资源，而不是电网交互策略；如果差距超过 30%，说明绿氨的竞争力高度依赖低廉的电力和电网配套条件。

5. 常见问题与调试排查实录

5.1 Cplex安装与Yalmip调用的坑

最先遇到的坑几乎都集中在求解器配置环节。Yalmip 报 “No solver specified” 通常不是因为没有 Cplex，而是 Cplex 接口没被 Yalmip 探测到。解决方法是仔细检查 matlab 路径里是否同时包含了 Yalmip 和 Cplex 的 matlab 子目录，然后在命令窗口重新执行yalmiptest。

另一个高发问题是 Cplex 许可证失效。社区版免费 Cplex 支持的变量数有限，如果模型规模超出限制，求解过程会直接终止。这个问题的判断方法是看报错堆栈里有没有license相关的提示。解决办法是换正式授权的 Cplex 或者 Gurobi。

版本兼容这块，Cplex 每年发布一个大版本，对应的 Matlab 支持列表在安装目录的 README 里写得很清楚。我的实际经验是：Cplex 12.9～12.10 配 Matlab R2019b～R2022b 基本稳定；太新的 Matlab 版本偶尔会遇到 JVM 兼容问题，可以在 Matlab 启动时加-nojvm规避，但会牺牲一部分图形界面功能。

5.2 模型不可行的定位思路

求解结果直接告诉你 “infeasible” 是家常便饭。第一次跑离网模式时，我花了整整半天才找到问题根源。

不可行问题最直接的排查工具是 Cplex 的 IIS（Irreducible Infeasible Subset）诊断。在 Yalmip 里可以这样处理：

ops = sdpsettings('solver', 'cplex', 'cplex.iis', 1); [diagnostics, infeasible] = optimize(Constraints, Objective, ops);

Cplex 会返回一个极小不可行子集，指向导致冲突的那几个约束。我在复现时发现，离网模式下最常触发 IIS 的是氨产量目标约束：当连续设定的极端日没有风光出力时，无论储氢罐多大，模型都可能在某个约束组合上找不到解。

定位到源头后，我一般用两种方式处理。第一种是把氨产量约束从硬约束改成软约束，加一个缺产惩罚变量，让模型在极端条件下可以选择减产而不是直接无解。这样既保留了模型的可行性，又能从惩罚变量的最优值看出系统在什么时段、多大程度上存在供氢风险。第二种是检查初始氢量，储氢罐起始氢量不能设为零，否则第一个时段合成氨装置就没氢可用。

还有一个很容易忽略的点：储氢罐的容量单位是千克，制氢量的计算结果是千克，但产氨量的计算是吨，单位换算差了 1000 倍。我第一次跑的时候就是把m_nh3 == h_nh3 / k_h2_nh3 / 1000里的1000漏掉了，导致模型怎么算都不可行，查了半天才意识到是量纲问题。

5.3 求解时间过长与数值病态处理

在模型能求解但很慢的时候，瓶颈几乎总是出在二进制变量上。电解槽和合成氨装置的启停变量一旦和容量变量乘在一起，约束矩阵会变得非常难处理。

速度优化的第一原则是先查数值尺度。把年化投资成本用百万元做单位，功率用 MW，氢量用吨，储氢容量用吨，各变量的数量级尽量控制在 0.01～100 之间。Cplex 比较青睐这样的尺度，否则会出现大量的数值填充和退化，导致分支树爆炸。我试过模型中同时出现 0.001 和 1 亿这样跨度的数值，求解器光是预处理就花了十几分钟。

第二原则是严格控制 0-1 变量的数量。电解槽的启停变量可以保留，因为它对应真实物理约束；但像购售电互斥这种逻辑，如果分时电价设计得合理，购电和售电同时为正在经济上天然不成立，可以不加二进制变量，只在参数层面保证售电价低于购电价即可。

第三原则是分层求解。前面提到先把容量固定跑调度，得到初始解后再放开容量跑联合优化，这个方法在工程上非常实用。对于 8760 小时的大模型，这是从“跑不出结果”到“几分钟出结果”的关键操作。

最后分享两个小技巧

在这套模型的复现过程中，我最大的体会是：模型的难点不在数学形式，而在对物理系统的抽象程度。只有把电力平衡、氢物质平衡、设备运行边界这些物理规律写清楚，Cplex 才能成为一个真正有用的工具。

其中一个技巧是：场景削减之后一定要做全年 8760 小时的校核算例。削减后的典型日会把极端气象事件抹平，模型自动“乐观”。用同一个容量配置跑一遍全时序仿真，如果发现弃电率过高或者氨产量不足，就应该回到场景削减环节，把更极端的日场景加进去。

另一个技巧是：不要迷信某个求解器的默认参数。Cplex 默认的 mipgap 是 1e-4，对工程问题过于苛刻。设置mipgap = 0.01或0.02，往往能让求解时间从几小时降到几分钟，而目标函数值只差不到 1%。工程师的时间比那零点几个百分点的理论最优值值钱得多。