最近花了一周时间复现了一个“并离网风光互补制氢合成氨系统容量-调度优化”项目整套模型用Matlab调用Cplex求解器实现。这个方向最近在氢能圈子里讨论得很多但真正能一步跑通的代码并不多。我复现的这套模型同时考虑了容量配置和日内调度还区分了并网和离网两种运行模式基本覆盖了这类优化问题的全部核心难点。这篇文章把从建模、编码到调试的完整过程记录下来适合正在做风光储氢系统优化、合成氨路径经济性分析以及刚接触Cplex/Matlab联合求解的同学参考。1. 复现这个项目之前先搞清楚它在算什么1.1 风光互补制氢合成氨系统的能量链路先看物理链路风力发电机和光伏板发出来的电经过汇流和变压后一部分直接供给电解槽制氢一部分存到电池储能里还有一部分在并网模式下可以卖给电网。电解槽产出的氢气进入储氢罐之后与空分装置提供的氮气混合在合成氨工段里反应生成液氨。这条链路比单纯的风光制氢要复杂因为多了合成氨这个化工环节。合成氨对氢气的需求量不是固定不变的它取决于反应器的运行状态、负荷率以及氢氮比。所以系统优化不能只盯电力平衡还要做氢气的物料平衡和储氢罐的容量约束。换句话说发电侧、制氢侧、化工侧三个子系统之间通过“电-氢”这条主线耦合在一起任何一个环节的容量选小了整个系统的产能都会受限。我一开始以为这只是一个“风光配储能电解槽”的经典容量优化问题实际建模时才发现合成氨工段的引入让问题规模明显上升。它不光多了连续变量还多了设备启停、运行模式切换这类0-1变量必须用混合整数线性规划框架去描述。1.2 “容量优化”和“调度优化”到底在算哪两个问题容量优化解决的是“装多少”的问题风电机组装多少台、光伏板铺多少千瓦、电解槽配多大功率、储氢罐容积多少立方米、电池储能容量多少兆瓦时。调度优化解决的是“怎么运行”的问题在已经确定的设备容量下每个时刻风电和光伏出力怎么分配、电解槽负荷率调多少、储能充还是放、储氢罐什么时候存氢、合成氨工段开工多少。这两个问题不是独立的。容量配置决定了调度可行域的上限但容量大小本身又要靠调度运行结果来评价。如果把所有容量方案都建模成一个超大规模的单层优化问题变量和约束的数量会爆炸而且工程上很难让人理解结果。所以主流做法是双层优化外层迭代容量方案内层针对每个容量方案做全年或典型日的调度优化返回总成本再让外层判断哪个容量组合最优。我用一个生活化的类比来帮助理解买车的时候你要决定发动机排量容量问题但排量大了车就贵排量小了上坡肉还得看日常开车的路况来决定调度问题。你不能只说“我要大排量”因为还要算油费、保养费同样不能说“我要最便宜”因为你每天通勤的路况决定了真实能耗。这就是容量和调度要联合优化的原因。1.3 为什么非要用Cplex而不是手写线性规划这个模型里出现了很多整数变量电解槽的启停状态、储能充放互斥状态、合成氨工段的开停机状态甚至风电光伏的投建台数在离散化后也是整数。这是一个典型的MILP混合整数线性规划问题手写单纯形法或者拉格朗日松弛基本不可能在可接受时间内收敛到全局最优。Cplex是目前公认的MILP求解最稳定的商业求解器之一它对大规模整数规划的分支定界、割平面算法支持非常成熟。更关键的是Cplex在Matlab里有多种调用方式可以直接用Cplex的Matlab工具箱也可以通过Yalmip建模工具把优化模型“翻译”成标准形式再交给Cplex求解。我用的是YalmipCplex组合因为它能大幅减少矩阵拼接和约束书写的代码量后期调试时也能方便查看目标函数和约束。之前看过有人用fmincon直接跑这个模型其实不太合适。fmincon针对非线性问题而风光储氢氨优化的目标函数和大部分约束在合理假设下都是线性的非线性化反而会引入局部最优。Cplex的线性规划和整数规划求解器在这方面是碾压级的表现这也是我选它的根本原因。2. 数学模型目标函数与约束条件拆解2.1 目标函数费用最小到底包括哪些项这类系统优化的目标函数一般取“年等值总成本最小”就是把一次性投资成本折算到每年再叠加每年的运行费用。折算时需要确定设备寿命和折现率这个细节直接影响优化结果很多人忽略。常见的目标函数组成如下表所示成本分项含义单位示例对容量/调度的作用设备投资年值风电、光伏、电解槽、储氢罐、电池、合成氨装置等初始投资的等年值万元/年推动容量不要铺得过大运维成本各设备年运行维护费通常按投资比例或按发电量/产氢量计算万元/年影响高利用率和低利用率设备的选择购电费用并网模式下从电网买电的成本万元/年影响是否多配光储来减少购电售电收益并网模式下余电上网的收入万元/年鼓励在风光大发时段多发电惩罚成本弃风弃光、未满足供氢需求时的惩罚万元/年保证系统可靠性和资源利用率目标函数可以写成min设备年投资成本总和 年运行维护成本 年购电费用 - 年售电收益 年惩罚成本我在写目标函数时比较注意量纲统一。投资成本里的“万元”和购电单价里的“元/千瓦时”很容易混。建议提前约定功率单位用kW能量单位用kWh氢气量用kg时间用小时钱用万元。然后在Matlab里统一除以10000转换到万元避免数值量级差距过大导致Cplex数值问题。2.2 关键等式约束电平衡、氢平衡、氮氢比这部分是整个模型的核心骨架漏掉任何一个等式约束算出来的结果在物理上都是不可行的。电功率平衡是最直观的约束。每个时刻注入系统的功率必须等于流出系统的功率公式大致如下风电出力(t) 光伏出力(t) 储能放电功率(t) 电网购电功率(t)电解槽耗电功率(t) 储能充电功率(t) 电网售电功率(t) 弃电功率(t)这里我额外加入了弃电功率变量。离网模式下没有电网购售电项如果风光出力超过负荷和储能吸收能力就必须弃掉。如果没有弃电变量约束会无解所以宁可多设一个变量来放松平衡约束再用目标函数里的惩罚项去约束它。氢气物料平衡同样关键。电解槽产氢量(t) 储氢罐放氢量(t) 合成氨工段耗氢量(t) 储氢罐充氢量(t)。这里的充放方向和储能类似需要保证储氢罐的容量上下限。合成氨工段还有一个物料比例约束因为氨的分子式是NH₃理论上每1 mol氮气需要3 mol氢气。折算成质量流量时需要根据具体压力和温度下的转化率修正。论文里一般会给一个“氢氮比”系数通常在3:1附近但实际合成回路会补入循环气所以系统的净耗氢量会略高于理论值。这个系数我从文献里取值时格外小心因为不同压力和催化剂条件下的数据差别很大。2.3 不等式约束设备出力上下限、爬坡、储能SOC不等式约束描述设备的运行边界。电解槽的输入功率通常不能低于某个百分比比如30%也不能高于额定容量这个区间用变量上下界写进Cplex比单独加两个约束更高效但也需要明确告诉求解器哪些是硬限制。电池储能需要同时约束充放电功率和SOC状态。SOC递推关系是类似SOC(t1)SOC(t)充电功率(t)充电效率/电池容量-放电功率(t)/(放电效率电池容量)。这个公式本身是线性的但有充放电互斥的逻辑约束同一时刻要么充电要么放电。传统做法是引入0-1变量充电状态为1时放电功率强制为0反之亦然。储氢罐的约束逻辑也差不多不过储氢罐的单位是kg充放的是氢气流量动态方程比电池简单一些不需要考虑容量衰减但储氢罐的初末状态必须保证一致否则整个调度周期里的氢量会对不上。设备爬坡约束在风电制氢系统里经常出现电解槽和合成氨压缩机在短时间内的负荷调整幅度有限。大部分论文为了简化会忽略爬坡但实际运行中频繁快速升降负荷会损伤电极和催化剂。如果要严谨需要在相邻时段之间加差值上限约束。Cplex处理这些约束没有任何问题只是变量和约束数量会翻倍求解时间会上升。2.4 决策变量离散化为什么出现整数变量这个模型里的离散变量主要来自三个地方设备启停状态、储能充放互斥、离散的容量建设方案。电解槽和合成氨装置的低负荷运行区间通常意味着启停决策当设备负荷低于下限时必须启动0-1变量来强制停机启动后负荷又必须高于下限。储能充放互斥直接用两个0-1变量控制。容量建设的台数和套数更是本来就是整数比如风电装2台还是3台。引入整数变量会让问题变成MILP求解难度比纯LP上升一个量级。Cplex内部的分支定界算法会反复解LP松弛、裁剪整数变量组合。这里我的实际体会是整数变量数量直接决定求解效率。如果你的典型日有24个时段储能、电解槽、合成氨各引入一个0-1变量那一个调度子问题就有72个整数变量Cplex处理起来还好但如果你把全年8760小时一起算整变量就是几百倍内存都会爆所以一定要做场景缩减。3. MatlabCplex环境准备与数据预处理3.1 环境版本匹配Cplex、Matlab、Yalmip的兼容组合我踩过最大的坑之一是版本兼容。Cplex的Matlab接口不是自动生效的需要把对应版本的cplex文件夹加到Matlab路径并且安装好mex文件。推荐组合Matlab R2021b以上Cplex 12.10或12.9Yalmip最新版。不能直接用Cplex 12.6去配R2023a因为mex文件是针对特定Matlab版本编译的。Cplex安装目录下有一个/matlab文件夹进入后在Matlab命令行执行installmex或者自己addpath。Yalmip只是建模层它通过solve函数自动识别已安装的求解器检测命令是yalmiptest。如果执行sol optimize(...)时报错“Could not locate solver CPLEX”第一反应不应该是重装而是先运行which cplex看看Matlab是否真的找到了Cplex的mex接口。如果没有再到Cplex的安装包里重新编译mex文件。3.2 风光出力时序数据的来源与典型场景缩减这个模型对风光出力数据非常敏感。同一个容量配置下换一组风速和辐照度数据最优结果可能完全不同。复现时我建议先统一数据来源避免把时间耗在无意义的数据差异上。一篇论文里常用的数据来源包括风电场实际SCADA数据、气象站观测数据、开源再分析数据比如NASA POWER或者用威布尔分布和Beta分布分别模拟风速和辐照度。如果论文里没有给实际数据用模拟数据复现也可以关键是保证时序特性比如风光互补性。为了降低调度计算的规模我不会直接跑8760小时而是做典型日缩减。最常用的是K-means聚类把全年日曲线聚成4类场景代表春夏秋冬或者不同天气类型。每个场景赋予权重最后在目标函数里按权重累加。这样调度的变量数从8760降到96或144个时段Cplex求解时间从几分钟缩短到几秒而且结果基本能保持原问题的趋势。聚类后还要检查一件事相邻时段的风光出力变化是否剧烈。如果聚类中心曲线的波动太大可能不符合实际环境也会让调度结果出现频繁启停的假象。我一般会在聚类后做一次3点滑动平均去除毛刺。3.3 参数表设计把论文里的量纲统一起来参数表是整个模型能不能跑通的地基。我见过很多复现代码报错最后查出来是风电单位造价写成了“元/kW”但光伏单位造价写成了“万元/MW”量纲差了一万倍结果容量配置直接跑偏。我做参数表时统一用如下框架参数名称数值单位备注风电单位投资6500元/kW包含基础、塔筒、安装光伏单位投资3500元/kW包含支架、逆变器电解槽单位投资1200元/kW碱性电解槽电池储能单位投资1500元/kWh容量成本储氢罐单位投资300元/kg按容量成本合成氨装置投资4500元/kW(氨)按氨产量功率折算所有成本换算成年值时要加折现率r和设备寿命n年金系数APr(1r)^n/((1r)^n-1)。这个公式我在代码里单独写成一个函数避免在循环里重复手算。4. 容量-调度双层优化主程序实现4.1 双层框架外层容量枚举/启发式内层Cplex调度复现时我采用的是最稳妥的全枚举内层调度框架。容量组合网格化风电容量从0到上限按步长500kW取光伏容量按500kW取电解槽容量按1000kW取这样网格节点数不多。对每一个容量组合调用内层调度求解得到该方案下的年运行成本和购售电费用累加年投资成本得到总成本。遍历所有网格后选择总成本最小的组合。全枚举在大规模问题上效率不高但优势是结果的全局性有保障。当容量维度不超过3个时这个方法非常可靠。如果以后设备种类增加到5种那就要换成遗传算法或者粒子群来做外层优化但那时需要额外处理启发式算法的随机性问题。我建议先枚举复现论文确认结果能对得上再考虑算法升级。4.2 内层调度优化核心代码我用Yalmip建模关键代码如下可以配合注释理解% 假设已经读到场景数据 % T: 时段数, wind(T), pv(T), price(T) % 已知容量参数: cap_el, cap_bat, cap_h2, cap_pem % 决策变量 P_el sdpvar(T,1); % 电解槽输入功率 P_bat_ch sdpvar(T,1); % 电池充电功率 P_bat_dis sdpvar(T,1); % 电池放电功率 P_grid_buy sdpvar(T,1); % 购电 P_grid_sell sdpvar(T,1);% 售电 P_curt sdpvar(T,1); % 弃电 H_el sdpvar(T,1); % 电解槽产氢量 H_st_ch sdpvar(T,1); % 储氢罐充氢量 H_st_dis sdpvar(T,1); % 储氢罐放氢量 H_amm sdpvar(T,1); % 合成氨耗氢量 SOC sdpvar(T1,1); % 电池SOC M_st sdpvar(T1,1); % 储氢罐储氢量 x_el binvar(T,1); % 电解槽启停 % 约束 Constraints []; % 电平衡 for t 1:T Constraints [Constraints, wind(t) pv(t) P_bat_dis(t) P_grid_buy(t) P_el(t) P_bat_ch(t) P_grid_sell(t) P_curt(t)]; Constraints [Constraints, 0 P_el(t) cap_el]; Constraints [Constraints, 0.3*cap_el*x_el(t) P_el(t) cap_el*x_el(t)]; Constraints [Constraints, P_bat_ch(t) cap_bat*1]; Constraints [Constraints, P_bat_dis(t) cap_bat*1]; Constraints [Constraints, P_grid_buy(t) 2000]; Constraints [Constraints, P_grid_sell(t) 2000]; end % 氢平衡 for t 1:T H_el(t) P_el(t) * el_eff; % 电氢转换系数 Constraints [Constraints, H_el(t) H_st_dis(t) H_amm(t) H_st_ch(t)]; Constraints [Constraints, 0 H_st_ch(t) 100]; Constraints [Constraints, 0 H_st_dis(t) 100]; Constraints [Constraints, H_amm(t) 0]; end % 电池SOC递推 for t 1:T Constraints [Constraints, SOC(t1) SOC(t) P_bat_ch(t)*eta_ch - P_bat_dis(t)/eta_dis]; Constraints [Constraints, 0 SOC(t1) cap_bat]; end Constraints [Constraints, SOC(1) SOC(T1)]; % 初末一致 % 储氢罐递推 for t 1:T Constraints [Constraints, M_st(t1) M_st(t) H_st_ch(t) - H_st_dis(t)]; Constraints [Constraints, 0 M_st(t1) cap_h2]; end Constraints [Constraints, M_st(1) M_st(T1)]; % 目标函数购电费用 - 售电收益 惩罚 cost_annual_ope sum(P_grid_buy.*price_buy)*dt - sum(P_grid_sell.*price_sell)*dt; punish 1000 * sum(P_curt); % 弃电惩罚 Objective cost_annual_ope punish; % 求解 ops sdpsettings(solver,cplex,verbose,2); sol optimize(Constraints, Objective, ops);这段代码是核心骨架。实际使用时还需要把电/氢单位统一比如电解槽产氢量的系数要合理折算避免产氢量巨大而储氢罐容量过小导致无解。4.3 外层容量配置循环与结果收集外层循环我习惯写成三层嵌套for分别遍历风电容量、光伏容量、电解槽容量。每次内层求解都调用一次optimize这个函数在默认设置下不会保留上次的模型但要注意Yalmip变量之间不要重名冲突。每次得到调度结果后我会把最优目标值、投资成本、购售电费用、弃电率、电解槽利用率全部存进结构体数组。等枚举完成后用find找到总成本最小的索引并输出对应的容量组合。我还会额外画一张成本变化曲面图用surf看风电/光伏容量的二维成本地形这样能直观看到最优区域是否平坦。这个环节容易犯的错误是外层循环把Cplex求解器的残余变量和约束带到下一轮。建议每个内层子函数独立创建一个Yalmip变量空间或者在函数开头插入clear相关变量。4.4 并网与离网两种模式的统一建模技巧并网和离网最大的区别是电功率平衡方程里有没有电网交换功率。并网时多一组购电和售电变量目标函数多对应费用和收益离网时这些变量必须为0目标函数里的购售电项目消失。我实现时用一个isGridConnected标志位控制建模过程。当离网时购电和售电变量直接不用创建同时要增设弃电变量因为离网下弃风弃光几乎不可避免。离网模式下为了保证供电可靠性储能容量通常会变大电解槽容量也会适当降低避免在无风无光时设备闲置。还有一个小细节离网模式如果要保证连续多天运行光做单一典型日调度是不够的因为储能SOC的初末状态容易掩盖跨日能量亏空。我通常会给储能和储氢罐增加5%的末状态可调整范围或者引入一个长周期的调度校核否则高负载天数接续时系统很容易无解。5. 复现踩坑实录与调试技巧5.1 Cplex无法调用或license报错这个问题在第一次跑通之前几乎必然出现。常见报错有两种一是Undefined function cplex说明mex文件没配好二是License过期需要重新设置环境变量。第一个问题好解决Cplex安装目录下运行installmex然后addpath到Matlab工作区重启Matlab再检测。License问题主要发生在用试用版的学生机上此时不要报错就放弃先检查运行cplexlicsetup是否成功。Yalmip会在调用optimize之前调用solvesdp内部检测求解器。我遇到过Yalmip识别出Cplex但实际调用失败的情况原因是Matlab的临时目录被清理过Cplex需要写log文件但没有权限。解决方法是手动给Matlab设置可写临时目录或者把路径移动到新目录里重新配置。5.2 调度子问题出现infeasibleCplex返回不可行时第一步不是调参而是用Yalmip检查约束系统。可以用diagnostics optimize(...)返回状态信息。如果是不可行我会把目标函数改成0只跑可行性问题然后用“松弛”方式快速定位冲突逐个把可疑约束删除观察是否可行。最常见的原因是SOC初始值和递推公式的系数不匹配。比如电池容量5000kWh但单位电价时段内充电功率却设成8000kW虽然上限没写错但SOC递推会在几个小时内越过边界。解决办法是检查每个变量的上下界是否与递推系数自洽。另一个高发问题是储氢罐初末状态不一致如果初始储氢量设成满罐但末状态也必须回满且合成氨工段全程高负荷耗氢就可能无解。解决方法是给初末状态约束加一个允许偏差比如M_st(1) 0.9*M_st(T1)。5.3 结果中制氢量/合成氨量对不上这个现象很普遍根子在单位换算。电解槽产氢量和合成氨耗氢量经常一个用kg一个用t或者一个用Nm³一个用kg。Cplex本身不关心单位只关心数值但错一个数量级调度结果就会有大量弃氢或氢量缺额。我在代码里把所有流量统一换算成kg。电解槽产氢量输入功率(千瓦时)*制氢电耗的倒数制氢电耗通常是4.5-5.5 kWh/Nm³再乘上氢气密度0.0899 kg/Nm³。合成氨耗氢量按氨产量和氢氮转化比折算。这里我建议把转化过程写成一个函数并加上单位注释确保前后一致。5.4 代码运行速度优化当典型日时段数从24扩展到96时求解时间会暴增。Cplex求解MILP的瓶颈在整数变量搜索变量越少越快。我常用三个加速技巧一是给所有变量设置明确的上下界减少求解器的探索空间。比如电解槽功率上限直接设成容量值储能充电功率上限设成额定功率不要留空。二是把明显的等式约束固化成变量关系。例如电解槽产氢量和输入功率之间如果是固定效率就可以直接代换成线性表达式不单独设一个变量和一个约束这样Cplex内部的presolve能更快地化简。三是先用LP松弛求解观察结果。把binvar改成sdpvar运行一次连续性松弛看目标趋势是否合理。如果LP结果和目标论文差别很大多半是模型少了关键约束而不是求解器问题。6. 结果怎么看从优化曲线到工程决策6.1 典型日调度结果怎么看跑通后不要只盯着一个最优总成本要把调度曲线逐条拉出来看。最典型的图是24小时的风电出力、光伏出力、电解槽功率、储能SOC叠加图。通过这个图能直观看到白天光伏大时电解槽是不是满负荷夜间风电是否有出力低谷储能是在什么节点充电、什么节点放电。我复现后发现最优调度并不是让电解槽一直满负荷而是在风光资源充足时尽量多制氢把氢气存在罐里等到合成氨工段需要时再释放。这说明储氢罐的价值不仅仅在时移氢量更是碳酸缓冲平台。如果只看电平衡很难发现这一点。另外要重点检查弃电时段的分布。如果弃电集中发生在某些特定时段比如夜间风电大发而电网无法消纳说明容量配置里储能或电解槽容量偏低。如果弃电全年均匀分布说明系统整体装机偏高可以适当减小风光容量来降成本。6.2 并网模式与离网模式的容量配置差异并网模式下电网像是一个无限容量的“蓄水池”系统不需要配置很大的储能只需要保证购电成本和售电收益之间的套利平衡就行。离网模式下储能和储氢罐就成了唯一的能量缓冲容量会显著增大尤其是储氢罐它往往是离网系统中最贵也最关键的设备之一。我跑下来的典型结果是并网最优方案里风电和光伏容量偏小因为电网可以兜底离网最优方案里风电光伏容量偏大目的是提高自发自用比例但同时弃电率也升高。这个规律并不意外但它说明“并网/离网”这两个边界条件对容量配置的影响非常大如果你要做一个“并离网可切换”的系统就不能只取一个中间值设计必须分别计算两个典型的配置方案。6.3 从复现到改进灵敏性分析和随机优化扩展论文复现只是第一步真正想用它做工程决策一定要加灵敏性分析。我会对电价、氢气价格、设备投资、风光资源乘子逐个做±20%变化观察最优容量组合的变化路径。这种分析能告诉决策者如果未来电价下降系统应该多靠电网还是多装储能如果制氢设备降价应该把电解槽容量提高多少。另一个方向是把确定性模型升级成随机优化。经典做法是把风光出力用多个场景加概率权重代替单一典型日目标函数变成期望成本最小。Cplex完全可以处理多场景MILP但变量规模会成倍增加。我的实践经验是先用聚类缩减到4-6个场景再跑随机优化这样既保留了不确定性特征又不会让求解时间失控到小时级。这类系统的边界条件实在太复杂没有任何一种模型能完全覆盖现实中的所有非线性因素。但是基于Cplex的MILP框架至少能把电-氢-氨这个主要能量和物料链路算清楚。写在最后的调试习惯我做这类复现项目时收获最大的一点是永远不要急着写完整代码先把数学问题在纸上拆成“目标、等式约束、不等式约束、整数变量来源”四个板块。Cplex并不是魔法它只是在所有可行解里搜索最优。如果模型本身有物理漏洞再快的求解器也救不回来。调试顺序上我建议先跑一个简单场景把时间维度压缩到6小时把容量维度固定成一组小参数确认内层调度能够可行。然后在逐渐扩大时段和容量维度。每扩大一次都要对比上一轮的目标值和调度曲线确认没有跳变。这个过程虽然慢但能帮你建立对模型的直觉后面排查异常结果时会省下很多时间。最后再分享一个小经验保存每次的求解日志和参数快照。同一套模型换了输入数据之后的目标值变化往往能暴露出数据清洗问题。光靠眼睛看图表容易漏落到数字上才安全。