风光互补制氢合成氨容量-调度双层优化与Cplex求解

最近在复现一篇关于风光互补制氢合成氨系统的容量-调度优化论文，用的求解器是Cplex，代码环境是Matlab。断断续续啃了两周，踩了好些坑，也把整个系统的建模逻辑捋清楚了。这篇文章就把这次复现的完整思路、模型构建、Cplex接入方式和实际调试中遇到的问题都整理出来，给准备做或者正在做类似方向的朋友做个参考。

这套系统值得复现的原因在于：它不是单纯的微电网容量配置，也不是简单的制氢调度，而是把"风电+光伏→电解制氢→储氢→合成氨"这条完整的能源化工链条放进一个优化框架里，同时考虑容量规划（装多少风机、光伏、电解槽、储氢罐）和运行调度（每个小时风电出力怎么分配、电解槽开多大、氢怎么用）。这类问题在综合能源系统研究和工程规划里非常有代表性，用Cplex求解是因为模型最终可以转成混合整数线性规划（MILP）来精确求解，相比启发式算法，最优性有保证。

无论你是做新能源消纳、绿氢化工、综合能源系统优化，还是单纯想看看Cplex在Matlab里怎么解这类大规模调度问题，这篇文章应该都能提供一些实打实的帮助。

1. 看标题之前，先搞懂这条能源链路的优化对象

风光互补制氢合成氨系统，说白了就是把不稳定的风光电转化为可储存、可运输的化学品。这个系统不是几个设备随便连起来，而是涉及到一条多级能量传递链路，每一级都有容量配置和运行约束的问题，优化起来远比单一的风光电站或者单独的制氢厂复杂。

1.1 从风光发电到合成氨的完整工艺流程

整个系统的能量流动是单方向的：风力发电机组和光伏阵列产生电能，电能通过整流变流设备输送给电解槽，电解槽将水电解为氢气和氧气，氢气经过纯化、加压后进入储氢装置，然后送入合成氨装置与氮气反应生成合成氨。最终产品是氨，它的好处非常多——常温常压下是液体，储存运输方便，既可以直接用作肥料原料，也可以作为一种氢的载体。这也是为什么很多风光资源丰富的地区都在布局"风光制氢合成氨"一体化项目，本质上是解决风光电的消纳和氢的储运难题。

这个系统里每一级设备都不是独立存在的，上一级的输出直接决定下一级的可用资源。比如电解槽的耗电量决定了它需要风光出力的支撑，储氢罐的容量决定了它能缓冲多少氢气的供需不匹配，合成氨装置的年运行小时数又反过来约束了系统的有效产出。优化的时候，任何一个环节的容量变化都会传导到整个系统，这就是典型的"牵一发而动全身"。

用一句话来概括这个优化问题的本质：在满足制氢和合成氨需求的前提下，用最小的年化总成本去配置各个设备的容量，并且给出每一小时的最优运行策略。

1.2 容量优化和调度优化分别回答什么问题

这里要区分两个层次的优化，很多第一次接触的人容易混。

容量优化（Capacity Optimization）回答的是"建多大规模"的问题。风电装机多少兆瓦、光伏装机多少兆瓦、电解槽额定功率多大、储氢罐能装多少立方米的氢气、合成氨装置的年产能是多少。这些决策变量是设备规格层面的，一旦定了，项目初投资就基本确定了。容量优化通常在年尺度或典型日尺度上进行，因为它关心的是长期规划方案。

调度优化（Scheduling/Dispatch Optimization）回答的是"每个时刻怎么运行"的问题。给定设备容量之后，在8760个小时中每个小时风电和光伏实际出力多少、电解槽输入功率多少、储氢罐充放氢多少、合成氨装置是满负荷跑还是降负荷跑。这些决策变量是运行层面的，对应的是系统每天的运行成本、购电成本、收益等。

容量优化和调度优化之间存在嵌套关系：不同的容量配置下，同一个调度策略产生的运行成本是不同的，而运行成本又会反过来影响最优容量选择。这就是为什么说这是一个"容量-调度"双层耦合优化问题。在数学上，常见的处理方式是用双层优化结构，外层负责容量，内层负责调度，外层每给定一组容量方案，内层就求解一次最优调度并返回运行成本。

这里是典型的设置示例（容量变量记作X_cap，调度变量记作X_op，外层目标函数是年化投资+年运行成本的最小化，内层目标函数是给定容量下的最小化运行成本）：

外层容量优化： min (年化投资成本 + 年运行成本) s.t. 容量取值约束

内层调度优化（在外层给定容量后调用）： min 年运行成本 s.t. 逐时功率平衡、设备运行约束、储氢动态约束

容量和调度不是分开拍脑袋定的，而是在同一个目标函数下互相制约、共同寻优。这也是这类系统分析比单层优化要麻烦的核心原因。

1.3 为什么选Cplex而不是遗传算法

这个问题在复现论文时经常被问到。很多文献里做容量优化用的是粒子群、遗传算法或者NSGA-II，因为容量变量和调度变量混在一起形成非线性、非凸问题时，数学规划工具不好直接处理。但这篇标题明确写的是"Cplex求解"，说明它走的是精确求解路线——把问题建模成MILP，然后用Cplex的分支定界算法去求全局最优解。

用Cplex的优势非常明显：它有全局最优性证明，只要模型建得对，收敛到的就是数学意义上的最优解，不需要像启发式算法那样反复调种群大小、变异率、交叉率，也不用担心收敛到局部最优。代价是建模时必须做线性化或者近似处理，把非线性关系转成线性约束和整数变量，模型规模上去之后对内存和求解时间的要求也高。

以我的复现经验来看，如果系统规模控制在单节点、数十个设备变量、8760个时段，Cplex直接求解MILP是完全可行的；如果再加设备数量或者把时间分辨率细化到15分钟，就需要慎重考虑模型规模和求解性能了。后面具体说怎么控制规模。

2. 容量-调度双层优化模型的建模逻辑与数学表达

要做优化，第一步就是把实际问题翻译成数学语言。这个翻译过程看起来简单，实际上很多细节处理不好，最后求解出来的结果根本不能用。下面按我的建模顺序拆开讲。

2.1 外层容量规划的决策变量和目标函数

先列出外层的决策变量集合。以我复现的模型为例，容量变量一共5组：

• 风电装机容量 P_w_cap（单位MW）
• 光伏装机容量 P_pv_cap（单位MW）
• 电解槽额定功率 P_ele_cap（单位MW）
• 储氢罐最大储氢量 V_h2_cap（单位kg或Nm³）
• 合成氨装置产能 P_amm_cap（单位t/h或kg/h）

外层目标函数是最小化年化总成本，包含三个部分：

第一是年化投资成本。设备的一次性建设投资不能直接加进年成本里，需要乘以资金回收系数（CRF）换算成年均等值支出。CRF的计算公式是：

CRF = r × (1 + r)^N / ((1 + r)^N - 1)

其中r是贴现率（一般取6%到8%），N是设备生命周期（风电和光伏一般取20年，电解槽取15年，合成氨装置取20年）。年化投资成本就等于每个设备的单位投资成本×容量×CRF，加总之后得到总投资年均支出。

假设风电单位投资成本k_w、光伏k_pv、电解槽k_ele、储氢罐k_h2s、合成氨装置k_amm，那么年化投资成本表达式为：

C_inv_annual = CRF_w × k_w × P_w_cap + CRF_pv × k_pv × P_pv_cap + CRF_ele × k_ele × P_ele_cap + CRF_h2s × k_h2s × V_h2_cap + CRF_amm × k_amm × P_amm_cap

第二是年运行维护成本。一般按设备投资的一定比例估算，比如风电运维费取总投资2%左右，光伏取1.5%，电解槽取3%到5%。

第三是年运行成本，这部分得靠内层调度优化返回：包括从电网购电的费用、弃风弃光的惩罚成本（如果有）、设备启停损耗等，再减去卖电收益（并网场景）。

外层约束主要是容量的上下限范围，比如风电装机不允许超过当地可开发资源的上限，电解槽功率和风光总装机之间可以有粗略的比例关系约束，储氢罐容量要覆盖合成氨装置连续运行的用氢需求等。

2.2 内层调度优化的决策变量与逐时约束

内层调度优化的核心是对每个时段t（t=1,...,T）做运行决策，决策变量包括：

• 风电上网出力 P_w_t
• 光伏上网出力 P_pv_t
• 电解槽输入功率 P_ele_t
• 储氢罐充放氢流量 H_in_t、H_out_t（一般统一成储氢量变化量）
• 合成氨装置耗氢量 H_amm_t
• 并网场景下还有购电功率 P_buy_t、售电功率 P_sell_t

内层的功率平衡约束是系统运行的核心：任何时刻风电和光伏的实际出力加上购电功率（如果并网），必须等于电解槽耗电、合成氨装置辅助用电、再加上对外售电。写成线性约束是：

P_w_t + P_pv_t + P_buy_t = P_ele_t + P_aux_t + P_sell_t

所有的功率变量都含时间下标，这一条约束在8760个小时里每个小时都要成立。这是整个模型里约束数量的大头。

电解槽运行约束也比较关键。电解槽不是任何功率下都能运行，它有最低运行负荷率，一般需要在额定功率的10%到20%以上才能稳定运行，同时不能超过额定功率。而且电解槽启动和停机的状态变化需要用二进制变量表示，否则无法约束"不能频繁启停"或者"最低运行时间"这类逻辑。引入二进制变量u_ele_t代表电解槽在第t时段是否开机，约束如下：

P_ele_min × u_ele_t ≤ P_ele_t ≤ P_ele_cap × u_ele_t

P_ele_min是电解槽的最小运行功率，P_ele_cap是额定功率。这组约束的作用是：如果u_ele_t=0，则P_ele_t强制为0；如果u_ele_t=1，则P_ele_t落在允许范围内。

储氢罐的动态平衡约束需要特别小心，它连接了制氢和用氢两侧。第t时段末的储氢量等于第t-1时段末的储氢量，加上本时段电解槽的产氢量，减去合成氨装置的耗氢量：

S_h2_t = S_h2_{t-1} + η_ele × P_ele_t / HHV_h2 - H_amm_t

其中，η_ele是电解槽效率，HHV_h2是氢的高位热值（用来把电功率折算成产氢量）。S_h2_t必须始终在0到V_h2_cap之间：

0 ≤ S_h2_t ≤ V_h2_cap

还有一个隐蔽的约束：一个完整调度周期的初始和末尾储氢量要一致（或者至少让末尾储氢量不小于初始值），否则优化结果会利用"把氢用光"来占便宜，导致结果在工程上不可行。这个问题我后面专门讲。

2.3 双层问题的求解策略：KKT转化还是迭代逼近

这里需要解释一下双层模型的求解方法。外层容量变量给一组值，内层调度问题就是一个带这些容量参数的MILP，可以直接用Cplex求。但外层容量变量要搜索最优值，不能随便枚举，因为连续变量组合空间巨大。

我复现时采用的是一种工程上非常实用的迭代求解策略：

1. 初始化一组可行的容量配置；
2. 将容量变量作为参数传入内层调度模型，用Cplex求解，得到该容量下的年运行成本；
3. 把内层求解得到的运行成本作为外层容量优化目标函数的一部分（等价于外层每评估一组容量就会调用一次内层求解）；
4. 外层用合适的搜索算法（可以用Cplex的外层MILP建模，也可以使用粒子群之类的元启发式算法驱动，但考虑到标题明确说Cplex求解，更稳妥的方式是把容量变量也建成MILP，和内层约束做适当整合后一并求解）；
5. 检查前后两次迭代的最优目标函数值差值，小于收敛阈值就停止。

还有一种理论更严谨的方式是KKT条件转化：当内层问题是线性规划时，可以写出内层问题的KKT条件，把内层问题整体等价为一组约束塞进外层的单层MILP中。这个方法的问题在于引入互补松弛条件后模型变成非线性，还得用大M法线性化，模型规模膨胀得很厉害，实际求解效率并不理想。而且如果内层是MILP（有整数变量），KKT条件就不成立了。

所以我的建议是：如果内层没有设备启停整数变量，KKT转化值得尝试；但只要有电解槽启停这些二进制变量，老老实实用迭代法或者直接构建单层MILP，反而更快。

3. Cplex接入Matlab的环境配置与常见坑

既然标题明确了Cplex求解，Matlab代码实现的前提就是把Cplex的Matlab接口调配好。这一步看起来简单，实际是很多初次接触的朋友卡壳最多的地方。这里把完整的配置过程和报错处理写清楚。

3.1 安装与路径配置的完整步骤

第一步是安装CPLEX Optimization Studio。IBM官网上可以下载最新版本，教育版用户可以申请免费学术许可证。社区版（Community Edition）也可以免费用，但它有变量和约束数量的限制（通常限制在1000个变量和1000个约束以内），复现这个例子的话模型规模很可能超限，建议直接用完整版或者学术版。

安装完成后，打开Matlab，通过下面的命令把Cplex的Matlab API路径加进去。以CPLEX Studio 12.10在Windows下的默认安装路径为例：

addpath('C:\Program Files\IBM\ILOG\CPLEX_Studio1210\cplex\matlab\x64_win64');

路径取决于具体的安装版本和目录。设置好之后验证一下：

cplex = Cplex('test'); cplex.solve();

如果能够正常创建并求解一个空模型，说明接口已经通了。

为了避免每次重启Matlab都要重新设置路径，建议用Matlab的预设路径功能把上面的目录永久保存，或者在项目启动脚本startup.m里写上addpath语句。

3.2 用Cplex类构建MILP模型的基本代码骨架

Cplex在Matlab中的API比较简洁，核心是用Cplex对象的方法添加目标函数、变量和约束。以一个小规模调度模型为例，代码骨架如下：

% 创建模型对象 cplex = Cplex('scheduling_opt'); cplex.Model.sense = 'minimize'; % 求最小值 % 目标函数系数向量 f（与变量顺序一一对应） cplex.Model.obj = f; % 变量下界、上界和类型 cplex.Model.lb = lb; cplex.Model.ub = ub; cplex.Model.ctype = ctype; % 'C' 连续变量，'B' 0-1变量，'I' 整数变量 % 线性约束: lhs <= A*x <= rhs % 将等式约束拆成 lhs == rhs 即可 cplex.Model.A = A_sparse; % 注意要用稀疏矩阵 cplex.Model.lhs = lhs; cplex.Model.rhs = rhs; % 求解 cplex.solve(); % 提取结果 x_opt = cplex.Solution.x; obj_value = cplex.Solution.objval; cplex.Status

对于8760个时段的调度模型，变量数量轻松超过数万个，约束矩阵用稀疏矩阵存储是必须的，不然内存直接爆掉。Cplex自己的底层引擎支持大规模稀疏线性规划，但Matlab端传入的稀疏矩阵格式要处理好，每个变量和每条约束的顺序要和目标函数系数严格对应。

3.3 配置过程中最常见的三类报错

我在配置时遇到并且帮别人解决过的报错基本集中在三类。

第一类是"未定义变量cplex或类Cplex"。这个报错九成是因为没加路径，或者路径加错了层级。CPLEX安装目录下包含了cplex/matlab/x64_win64、cplex/matlab/demo等多个子目录，要用x64_win64那个版本匹配Matlab位数。升级Matlab版本后可能还需要重新配置一次。

第二类是许可证相关的报错，比如"No valid CPLEX license"或者"CPLEX Error 1016"。这个通常是许可证服务没启动或者环境变量CPX_LICENSE_FILE没设置。如果用的是许可证服务器，需要在系统环境变量里指定"server@host";如果是本地节点式许可证（node-locked），安装的时候会自动配置，但也容易因为Matlab以管理员身份运行时权限路径不一致导致找不到。

第三类是模型求解时Matlab卡死或者内存溢出，这个不一定算"配置错误"，更多是模型规模问题。解决办法是检查约束矩阵的稀疏性，去掉多余的变量，适当把时间精度从1小时扩展到2小时或4小时来减少规模。还有就是用分段线性化代替非线性约束以减少整数变量。

4. 并网与离网场景的建模差异与结果对比

标题里"并_离网"三个字，意思是并网和离网两种运行场景都要分析。这两种场景虽然系统主体相同，但数学建模差别非常大，优化结果也会呈现出有意思的规律。

4.1 并网模式的功率平衡与购售电建模

并网模式下，系统与大电网之间存在能量交换。每个时段允许从电网购电，也允许向电网售电。这个看似简单的改动，实际会显著改变优化结果。

并网模式的功率平衡约束改为：

P_w_t + P_pv_t + P_buy_t = P_ele_t + P_aux_t + P_sell_t

新增的变量P_buy_t和P_sell_t都有边界约束，购电功率不能超过与电网签订的协议容量上限，售电功率也有限制，购电和售电不能同时发生，这个逻辑约束用二进制变量可以表达：

P_buy_t ≤ M_pbuy × u_buy_t P_sell_t ≤ M_psell × u_sell_t u_buy_t + u_sell_t ≤ 1

经济上讲，购电电价一般按峰谷平三段或者分时电价来计算，售电电价又分上网电价和市场化交易电价。引入购售电之后，系统实际上多了一个"电网储能"的调节手段——风光出力大发时可以把多余电卖出去，风光出力不足时可以从电网买电维持电解槽运行。这会直接影响最优容量配置：离网时需要靠增加储能或者提高风光装机来保证供电，并网时可以靠电网兜底，所以往往并网方案的风光装机容量和储氢容量会比离网方案小，但多出一项购电成本。

目标函数里多了一组购售电费用项：

C_grid = Σ_t (price_buy_t × P_buy_t - price_sell_t × P_sell_t)

分时电价政策对调度策略影响很大。我实测的结果是，在峰谷电价差足够大时，优化模型会自动选择在低谷时段多买电制氢、高峰时段少用电甚至卖电，电解槽的利用曲线明显跟着电价走。

4.2 离网模式的供电可靠性约束

离网模式完全不同。没有任何外部电网支援，任意时刻所有负荷只能由风电和光伏承担。功率平衡约束直接简化为：

P_w_t + P_pv_t = P_ele_t + P_aux_t

这看上去只是把P_buy_t和P_sell_t去掉的问题，但背后藏着极大的建模差异。离网模式下系统必须有足够的容量冗余应对风光出力的波动，否则就会出现供电不足。为了保证系统可行，需要引入可靠性约束，常用的做法是指定最大允许的失负荷概率（Loss of Power Supply Probability，LPSP），或者设置系统最小备用容量系数。

以LPSP约束为例，需要引入一个表示切负荷量的变量P_cut_t，当风光出力不足时允许部分负荷被切除，约束设置为：

Σ_t P_cut_t / Σ_t (P_ele_t + P_aux_t) ≤ LPSP_max

同时功率平衡变成：

P_w_t + P_pv_t + P_cut_t = P_ele_t + P_aux_t

这样模型在极端天气时段可以"断电"，但是全年累计的缺电比例不能超过设定值。实际操作时，这个约束还会导致模型中出现P_cut_t和P_ele_t的乘积项（如果电解槽功率也要跟着削减），需要线性化，或者直接假定电解槽是可以灵活调节的负荷，缺电时优先削减电解槽功率，这在实际控制系统里也是合理的。

4.3 两种场景下的容量配置对比

我复现并网和离网两种场景后，得到的结果规律性很强。以一套典型的风光资源和负荷参数为例（此处参数来自常见工程案例区间，具体数值根据当地资源重新标定），两种场景下的最优配置差异大致如下：

离网场景下，储氢罐实际上承担了跨日甚至跨季的缓冲作用——新疆、内蒙这类地区冬季风光出力低，连续几天阴天或者无风的时候，只能靠储氢罐里攒下的氢维持合成氨装置运行。所以模型会自动把储氢罐容量配得比较大。

并网场景下，电网承担了部分缓冲作用，储氢罐可以适当缩小。总成本往往比离网低，但这时候要小心一个"伪优化"陷阱：如果购电电价过低，模型可能会大幅缩小风光装机，导致系统名义上是"风光互补制氢"，实际大部分电力来自电网，这与项目初衷相悖。解决办法是给风光发电占比或者系统综合可再生能源利用系数加个下限约束，比如要求风光发电量占总用电量的比例不低于80%。

5. 复现过程中的关键细节与调试心得

前四章把模型和求解环境都讲清楚了，这一章专门写我在代码复现和调试过程中遇到的、网上资料比较少提及的细节问题和排错思路。

5.1 时间尺度的选择：全年8760小时还是典型日

调度优化最理想的情况是直接把全年8760个小时全部建模进去，这样风光出力的季节特性、连续多天的储氢动态过程都能如实反映。但8760个时段的MILP模型规模相当大，尤其是还有电解槽启停二进制变量（8760个二进制变量）和储氢动态约束时，Cplex求解时间可能会到几十分钟甚至数小时。

我建议的折中方案是先用典型日法做快速验证，再用全年数据做最终校核。典型日可以通过K-means聚类从全年风光出力数据里选出来，比如春夏秋冬各选3到5个典型日（或典型周），然后按聚类占比加权计算年运行成本。这样模型规模能压缩一个数量级，迭代调试效率高得多。

要注意的一点：如果用典型日代替全年，储氢罐的跨日动态约束必须谨慎处理，因为典型日拼接边界上储氢量可能不连续。最简单的办法是把边界条件设成储氢量在周期首尾等值，并且把典型日之间的过渡视为瞬态。

5.2 储氢罐动态平衡约束的边界处理

这是我调试过程中掉坑最深的地方。储氢罐动态约束里如果处理不好"末状态等于初状态"这个条件，优化器会给出一个看似成本很低、实际完全无法运行的方案。

举个例子，假如全年最后一个时段模型发现储氢罐里的氢没有用完，它会想办法在最后一个时段加大合成氨耗氢量把氢清空，从而减少储氢罐的容量需求。这个操作在数学上是合法的，但实际操作中，下一年还要继续运行，你不能每年年末把罐清空。所以约束必须写成：

S_h2_T = S_h2_0

即周期末储氢量必须等于周期初储氢量，这样系统才能循环运行。这个约束加上之后，储氢罐容量会合理增大，模型结果也真正可落地。

另一个细节是储氢量的量纲。电解槽产氢量的物理量和储氢罐容量直接相加减，必须统一单位。我用的单位是kg。电解槽输入功率P_ele_t（单位kW）乘以效率η_ele再除以单位电耗（kW·h/kg），就得到小时产氢量（kg/h）；合成氨装置的耗氢量则按单位氨产品耗氢量折算成kg/h。千万不要同时用Nm³和kg混合建模，量纲错乱排查起来极其痛苦。

5.3 求解速度优化与可行性调试

模型第一次求解时，我遇到的情况是Cplex运行超过2小时还没收敛。排查下来主要是三个原因导致的。

第一个是二进制变量过多。8760个时段的电解槽启停变量加上购售电状态变量，加起来接近两万个。解决办法是用Cplex求解器的参数控制收敛性：设置相对MIP gap（比如0.5%或1%），让求解器在满足精度要求时提前停止。实际工程应用中，0.5%的gap对容量规划决策完全够用。

cplex.Param.mip.tolerances.mipgap.Cur = 0.005; cplex.Param.threads.Cur = 8; % 开启多核

第二个是约束矩阵条件数太差。有些约束里出现几万倍差距的系数（比如投资成本是万元量级，储氢量是公斤量级），数值求解时容易出问题。解决办法是对模型做无量纲化或者量纲归一化处理，把所有成本项统一成同一量级，把功率、容量都折算到基准值上。

第三个是模型不可行。刚建完模型求解时，Cplex可能会报"Model is infeasible"。不要急着一行行查代码，直接用Cplex的IIS（Irreducible Inconsistent Subsystem）功能快速定位不可行约束集：

cplex.optimize(); if strcmp(cplex.Status,'infeasible') cplex.refineConflict(); disp(cplex.Conflict); end

冲突分析结果会直接告诉你哪几条约束同时满足不了，比如储氢罐容量上限和初始储氢量约束冲突，或者功率平衡约束与负荷上下限冲突，问题一下就定位到了。

5.4 复现论文时的参数来源与结果验证

复现阶段最容易出现的问题是参数张冠李戴。不同论文里的风电单位投资成本、电解槽效率、电价曲线差异非常大，直接用别人的参数会导致结果不可比。

我的做法是建立一个参数清单表格，每个参数标明来源（论文、工程报告、假设值），并记录敏感性分析范围。

关键参数的常见参考范围（具体值应以原文或实际项目为准）：

结果验证方面，我会检查几个工程常识来判断模型输出是否合理：电解槽年利用小时数是否落在1000到7000小时的合理区间；储氢罐平均储氢水平是否处在20%到80%的正常波动带；弃风弃光率是否与风光容量比例匹配。如果这些指标明显异常，通常说明模型约束或者参数有问题，得回头调试。

6. 一套可直接修改运行的Matlab+Cplex求解骨架

前面讲了大量理论和调试心得，最后给出一份可以在自己机器上跑起来的最小实现骨架。这份代码不是论文复刻的完整版，而是把容量-调度两层优化和Cplex接口的核心结构串起来，你拿到之后改参数、加约束就能用。

6.1 主程序结构与数据组织

主程序分三段：参数初始化、内层调度求解函数、外层容量优化驱动。为便于修改，我把参数全部集中在一个结构体里。

%% 参数初始化 params.T = 24; % 调度周期时段数，测试时先用24h params.dt = 1; % 时段长度(h) params.P_w_cap = 30; % 风电容量(MW) params.P_pv_cap = 20; % 光伏容量(MW) params.P_ele_cap = 25; % 电解槽额定功率(MW) params.P_ele_min = 0.1 * params.P_ele_cap; % 最小运行功率 params.V_h2_cap = 2e4; % 储氢罐容量(kg) params.eta_ele = 0.62; % 电解槽效率 params.HHV_H2 = 33.3; % 氢气高位热值(kWh/kg) 约39.4kWh/kg，这里按实际情况调整 params.Elec_per_kg = 4.5; % 电解水制氢单位电耗(kWh/kg H2, 约54kWh/kg，实际系统效率约70%) params.H2_per_NH3 = 0.178; % 合成氨单位耗氢(t H2/t NH3) params.NH3_rate = 0.5; % 合成氨装置每小时产氨能力(t/h) params.price_buy = [0.3, 0.6, 0.9, 0.3*ones(1,21)]; % 24小时购电电价(元/kWh) params.price_sell = 0.25 * ones(1,24); % 售电电价(元/kWh)

这里要说明一下，具体电耗数值要根据你实际的电解槽效率和单位换算重新标定。模型运行正确与否的关键在量纲统一。

6.2 内层调度模型的Cplex构建函数

内层调度就是给定容量参数后，用Cplex求解一个整数小时级的运行优化问题，返回最小运行成本和最优调度结果。

function [cost, x_opt] = solve_dispatch(params) % 创建Cplex对象 cplex = Cplex('dispatch'); cplex.Model.sense = 'minimize'; % 变量索引组织（以24时段为例） % 变量顺序：P_ele(1:24), P_w(1:24), P_pv(1:24), P_buy(1:24), P_sell(1:24), % H2_sto(1:25), H2_amm(1:24), u_ele(1:24), u_buy(1:24), u_sell(1:24) n_T = params.T; % 变量总数 n_var = 8*n_T + n_T + 2*n_T; % 这里按实际变量个数填写 lb = zeros(n_var, 1); ub = inf(n_var, 1); ctype = char(zeros(1, n_var)); f = zeros(n_var, 1); % 目标函数：购电费用-售电收益 % f中对应P_buy的系数 = params.price_buy * 1000 * params.dt % 对应P_sell的系数 = -params.price_sell * 1000 * params.dt % 按变量顺序填写目标系数、边界和类型 % 二进制变量用'B'，连续变量用'C' % 约束矩阵 A = []; lhs = []; rhs = []; % 逐时功率平衡约束 % P_w + P_pv + P_buy - P_ele - P_aux - P_sell = 0 % 风力与光伏的最大可用出力按各自时序曲线给定，写成<=约束 % 储氢动态约束 % S_h2(t+1) = S_h2(t) + 产氢 - 耗氢 % 产氢 = P_ele_t / params.Elec_per_kg % 耗氢 = params.H2_per_NH3 * 1000 * params.NH3_rate % 装配进cplex cplex.Model.obj = f; cplex.Model.lb = lb; cplex.Model.ub = ub; cplex.Model.ctype = ctype; cplex.Model.A = sparse(A); cplex.Model.lhs = lhs; cplex.Model.rhs = rhs; % 求解 cplex.solve(); cost = cplex.Solution.objval; x_opt = cplex.Solution.x; end

核心逻辑就这些，真正写的时候约束矩阵用循环逐时段添加，但为了提高构建效率，建议直接用块矩阵组装，别一条一条addRows，那样模型构建阶段会非常慢。

6.3 外层容量优化的迭代驱动与结果输出

外层驱动负责在给定的容量搜索空间内调用内层求解，并汇总年化总成本。我测试时用过最简单的方式是网格搜索加内层Cplex求解，把小规模问题跑一遍没什么问题。如果维度再高，就该用前面说的单层MILP或者元启发式外层驱动。

% 测试不同容量配置 cap_candidates = [20, 25, 30, 35, 40]; % 风电候选容量 results = []; for i = 1:length(cap_candidates) params.P_w_cap = cap_candidates(i); [op_cost, x_opt] = solve_dispatch(params); inv_cost_annual = CRF_w * k_w * params.P_w_cap + CRF_pv * k_pv * params.P_pv_cap + ... CRF_ele * k_ele * params.P_ele_cap + CRF_h2s * k_h2s * params.V_h2_cap + ... CRF_amm * k_amm * params.NH3_rate * 8760 * 0.001; total_annual_cost = inv_cost_annual + op_cost; results = [results; cap_candidates(i), total_annual_cost]; end % 输出成本最低的配置 [min_cost, idx] = min(results(:,2)); fprintf('最优风电容量: %.2f MW, 年化总成本: %.2f 万元\n', results(idx,1), min_cost);

跑通这个骨架后，就具备了继续扩展的基础：加入全年8760小时数据、更细致的电解槽运行模型、合成氨装置的启停策略、以及并网离网对比模块。每一步扩展都是对原有模型的小幅修改，不会伤筋动骨。

说到我自己的体会，做这类"容量-调度"耦合优化，最容易犯的错误不是约束不够，而是约束太多太死，导致模型不可行或者结果违背工程直觉。好的建模习惯是从简单模型出发，先跑通再逐步加细节，每加一条约束都重新验证结果的变化是否在预期范围内。另外写代码时统一量纲、统一单位、参数集中管理，能让你在后期排查时省下大量时间。

这篇文章把这次复现过程中的关键内容都覆盖到了：系统的工艺流程、双层优化建模逻辑、Cplex接入Matlab的方法、并网和离网建模差异、储氢动态约束的边界处理、求解性能调优，以及一个可以直接跑通的最小代码骨架。如果你正在复现类似的论文，照着这个思路走，应该能绕开不少我踩过的坑。