风光互补制氢合成氨这个方向近两年在能源系统优化里热度一直很高。我前阵子刚把一个并网/离网双模式下的容量-调度联合优化模型完整复现了一遍求解用的是Matlab调用Cplex跑通之后回头看看整个模型的搭建思路、约束处理方式、还有求解器调参的细节都比预想中要讲究。这篇文章就把整个过程拆开讲清楚从物理系统建模到Cplex落地求解再到并网离网结果对比适合正在做新能源制氢、电-氢-氨耦合系统规划课题的研究生也适合想用MatlabCplex做优化调度的工程师参考。1. 这个优化问题到底在求什么1.1 系统的物理架构与能量流先把这个系统的家底摸清楚。风光互补制氢合成氨系统物理上包含五个核心单元风电机组、光伏阵列、电解水制氢装置、储氢罐、氨合成塔。在并网模式下还会和外部电网发生功率交换所以电网联络线也算一个边界环节。能量流的方向其实很清晰风、光发电给电解槽供电电解槽产出的氢气进入储氢罐缓冲储氢罐出来的氢气与空分装置提供的氮气按3:1比例进入氨合成塔生成氨产品。中间每一步都有能量损失和物料损耗这些细节最终都会转化为模型里的效率系数和物料平衡约束。系统调制的核心矛盾在于时间尺度的不匹配。风光出力的波动是分钟到小时级别的电解槽动态响应相对快但氨合成塔属于化工装置要求进料流量平稳不能频繁启停也不适合深度变负荷。储氢罐就是用来解耦电解制氢和氨合成这两个动态特性差异巨大的环节。1.2 容量优化和调度优化为什么要联立这个问题拆成两层来看会更清楚。容量优化解决的是一次投多少钱、装多大设备的问题决策变量包括风机装机台数、光伏装机容量、电解槽额定功率、储氢罐有效容积、氨合成塔设计产能以及离网模式下需要配套的储能电池容量。这些变量决定系统的初始投资属于长期决策。调度优化解决的是设备装好之后每天怎么运行的问题决策变量包括各时段电解槽功率、储氢罐充放氢量、氨合成塔负荷率、并网模式下买卖电功率、离网模式下储能电池充放电功率。这些变量决定系统的运行成本和收益属于短期决策。容量和调度之间是强耦合的。你装多大的电解槽直接决定了调度层能安排多少电解功率反过来如果调度层发现某个容量配置下运行成本特别高上层就应该调整容量方案。所以必须把两层问题嵌套在一起求解也就是常说的容量-调度联合优化。我这次复现采用的是让Cplex一次性求解整个混合整数线性规划模型没有做外层启发式和内层LP的嵌套迭代这样处理的好处是能保证全局最优坏处是变量和约束规模明显变大对建模规范性的要求更高。1.3 为什么选Cplex而不是其他求解器做这类中大规模MILP问题Cplex几乎是我第一个想到的求解器。首先是求解性能确实扛得住带几百个二进制变量的MILP模型它通常在几分钟内就能收敛到1%以内的最优间隙其次Matlab工具箱里直接带cplex接口和YALMIP、CVX这种建模层的配合也成熟。当然也有替代方案。Gurobi在纯LP和MILP上跟Cplex不相上下而且许可证对学生更友好如果模型是凸的用CVX或者YALMIP内置的Sedumi、SDPT3也可以但处理大规模混合整数问题时收敛速度明显不如商业求解器。另外国内常用的一些开源求解器如SCIP、CBC能力也不弱但在这个模型规模下求解时间会拉得比较长。我的建议是如果学校或课题组有Cplex授权直接用Cplex如果没有优先试Gurobi代码改动很小实在没有商业求解器可用再考虑SCIP但要做好求解时间翻好几倍的心理准备。2. 数学模型怎么建目标函数与约束拆解2.1 目标函数设计目标函数采用全生命周期成本最小化把投资成本通过资本回收系数折算到每年的等年值再加上年运行成本形成一个可比的经济性指标。投资成本部分涵盖五项风机、光伏、电解槽、储氢罐、氨合成塔。这里有一个关键的折算处理每项设备单位投资成本乘以其容量再乘以资本回收系数CRF即[ C_{inv} CRF \times \sum_{i \in S} (c_i \times Cap_i) ]其中CRF的计算公式是常见的那套基于折现率和项目寿命期的年金公式。如果项目寿命是20年、折现率取6%CRF大约在0.0872左右。这个系数要把折现率取准了取低了会低估投资的真实成本取高了会过度惩罚初投资导致优化结果偏向于少装设备。我在这里踩过坑后面会细说。运行成本分两块并网模式下与电网交互的购电费用减去售电收益离网模式下主要是弃风弃光损失对应的虚拟成本。有储氢罐以后系统实际上有一定的跨时段能量转移能力所以运行成本的计算必须按小时展开不能简化成单个典型日的总量。如果合成氨产品有出售收益目标函数也可以改成净成本最小化即总成本减去氨销售收入。取决于你研究的是系统经济性最优还是给定氨产量下成本最低两者建模方式略有不同。我这里采用的是后者即固定年氨产量在这个前提下做容量和调度的联合优化。2.2 关键约束条件约束条件是这套模型的灵魂也是复现时最容易出错的地方。我把核心约束按物理环节拆开来梳理。第一条是各时段功率平衡约束。并网模式下每小时满足风电出力、光伏出力、电网购电功率、储能放电功率之和等于电解槽耗电、电网售电功率、储能充电功率之和。注意离网模式下购售电项全部置零储能的作用就被放大了。功率平衡约束必须严格满足这是整个调度可行域的地基。第二条是电解槽运行约束。电解槽不是任意功率都能运行它有一个最小运行负荷率一般碱性电解槽在20%到40%之间PEM电解槽可以低到5%。这个约束用一条不等式就能表达电解功率要么为0要么位于最小和最大出力之间。这就引入了二进制变量也是模型变成MILP的原因。如果忽略最小负荷率优化器会让电解槽以极低的功率运行这在物理上是不可能实现的。第三条是储氢罐动态平衡约束。储氢罐在每个时段末尾的储氢量等于上一时段末尾的储氢量加上本时段产氢量减去本时段供氨合成塔用氢量同时考虑电解槽的效率折算和储氢罐自身的损耗。这个约束是跨时段耦合的它把调度问题从单时段优化变成了真正意义上的动态优化问题。第四条是氨合成塔运行约束。氨合成塔有两个关键限制设计产能上限和最低稳定运行负荷率。因为合成塔要求连续稳定运行基本上一旦开机就不能停停机成本极高所以模型里通常会要求氨合成塔在整个时间范围内保持最低负荷以上运行这样储氢罐就必须保证持续供氢。2.3 非线性项的线性化处理风光出力、电解槽效率这些物理关系天然都是非线性的但Cplex求解MILP的强项在线性规划松弛非线性约束要么手动线性化要么引入分段线性近似不能让求解器直接处理非线性项。风电出力跟风速的关系是典型的分段函数低于切入风速不发电在切入风速和额定风速之间按三次方多项式上升在额定风速和切出风速之间保持额定超过切出风速停机。我对风速数据进行预处理把每个时段的风电最大出力直接作为输入参数传入模型这样就不需要在该时段内引入非线性约束这是一个很实用的简化技巧。电解槽的制氢效率随运行功率变化严格来说是一条非线性曲线。我采用分段线性化的方式把电解槽功率分成三段每段定义一个等效产氢效率然后引入二进制变量保证只能落在其中一段用大M法表达分段激活逻辑。这样模型就变成了MILPCplex可以直接处理。3. Matlab调用Cplex实现容量-调度联合优化3.1 开发环境配置与Cplex接入先把环境搭好。我用的是Matlab R2022a搭配IBM ILOG CPLEX Studio 12.10操作系统是Windows 10。安装CPLEX时需要注意务必勾选Matlab Integration组件否则之后还要手动配置路径。安装完成后在Matlab里设置路径。最简单的方式是把CPLEX的Matlab接口目录添加到搜索路径这一步我通常会写进启动脚本addpath(C:\Program Files\IBM\ILOG\CPLEX_Studio1210\cplex\matlab\x64_win64);配置好后验证一下是否生效在命令行输入cplex如果弹出CPLEX对象的帮助信息说明接口已经打通。这一步卡住的人不在少数常见的原因是路径写错或者Matlab版本跟CPLEX不兼容。CPLEX 12.10官方支持的Matlab版本是R2017b到R2022a太新的Matlab版本偶尔会出现mex文件不兼容的情况这时候要么换CPLEX版本要么在配置里选择使用旧版Matlab兼容模式。3.2 数据准备与参数预处理模型跑得对不对很大程度取决于输入数据准备的精细程度。风资源数据采用某地全年的逐小时平均风速序列光伏数据采用对应年份的逐小时太阳辐照度序列这两个时间序列是全年8760小时的。直接把全年8760个小时全部放进优化模型变量规模和求解时间都会失控所以需要进行典型日选取。我采用的是K-means聚类方法用风速、辐照度、气温三个特征把全年数据聚成4类每类取聚类中心最邻近的一天作为典型日并统计该类的天数占比作为权重。这样模型从8760小时压缩到96个小时4个典型日乘以24小时求解规模大幅下降但代表性损失不大。设备成本参数我用的是目前行业内比较认可的一组数据风机单位投资7000元每千瓦光伏4000元每千瓦碱性电解槽5000元每千瓦储氢罐按储氢量计约600元每标准立方米氨合成塔按年产能分摊后大约每吨氨800至1200元。运维成本按设备投资的固定比例取2%。电价数据方面并网模式参考国内某省级电网的一般工商业分时电价峰时段1.2元每千瓦时、平时段0.75元每千瓦时、谷时段0.35元每千瓦时同时设置上网售电电价0.4元每千瓦时。分时电价的引入会直接影响调度层的决策特别是电解槽的启停安排谷电阶段会主动加大电解功率把氢气存起来供白天氨合成塔使用。3.3 核心代码框架与求解流程模型代码我按功能模块组织数据加载模块、参数设置模块、变量定义模块、约束构建模块、求解模块、后处理模块。目标函数的构建是核心需要把所有变量的成本系数汇总成一个向量f然后传给Cplex。变量的排序要谨慎我在代码注释里详细记录每个变量在决策向量中的索引位置避免后续构建约束矩阵时串位。实际调试中经常出现约束对不上变量的问题排查起来非常痛苦所以从一开始就要把索引管理好。约束矩阵的构建采用逐条添加的方式。先初始化稀疏矩阵然后按约束类型逐行填入系数。这里有一个经验性的建议约束矩阵用稀疏格式存储否则96个时段加上几十个变量虽然内存不至于爆掉但求解器预处理效率会明显下降。求解代码的核心部分大致是这个流程% 创建Cplex模型对象 cplex Cplex(WindPV_H2_NH3); % 设置目标方向为最小化 cplex.Model.sense minimize; % 目标函数系数向量按连续变量在前、二进制变量在后的顺序排列 cplex.Model.obj f; % 不等式约束矩阵 Aineq变量系数矩阵 % 注意同时设置 lhs 和 rhs实现上下界约束 cplex.Model.A A_ineq; cplex.Model.lhs lhs_ineq; cplex.Model.rhs rhs_ineq; % 等式约束 cplex.Model.A [cplex.Model.A; A_eq]; cplex.Model.lhs [cplex.Model.lhs; b_eq]; cplex.Model.rhs [cplex.Model.rhs; b_eq]; % 变量边界 cplex.Model.lb lb; cplex.Model.ub ub; % 变量类型C代表连续B代表二进制 cplex.Model.ctype ctype; % 求解器参数设置 cplex.Param.timelimit.Cur 1800; cplex.Param.mip.tolerances.mipgap.Cur 0.005; cplex.Param.mip.strategy.startalgorithm.Cur 4; cplex.Param.threads.Cur 8; % 执行求解 cplex.solve(); % 获取结果 if cplex.Solution.status 101 || cplex.Solution.status 102 x_opt cplex.Solution.x; obj_opt cplex.Solution.objval; else error(求解未成功状态码%d, cplex.Solution.status); end状态码101和102分别对应最优解和达到最优间隙的可行解这两个状态都可以接受。求解完以后把决策变量按之前的索引规则还原成各个设备的容量和逐小时功率数据再跑一个校验脚本把所有约束重新检查一遍确认无违规量。3.4 求解性能调优的实用策略Cplex求解MILP的速度跟模型质量高度相关同样的数学问题建模方式不同求解时间可能差出一个量级。我调试过程中最见效的一个优化手段是给二进制变量设置合理的初始解。提前跑一遍启发式策略把容量变量固定到几个合理的经验值上然后只优化调度层得到一个初步调度方案再把这个方案中的二进制变量取值作为MIP start传给Cplex。这个操作可以显著削减分支定界树的搜索空间。另一种策略是适度放松MIP间隙容忍度。学术研究追求1e-4的间隙当然好但工程实践中0.5%到1%的最优间隙已经足够支撑决策分析。在代码里设mipgap为0.005求解时间通常会从几十分钟缩短到几分钟而目标函数值的差异往往不到0.3%。还要注意大M参数的取值。处理电解槽要么关要么在最小功率以上运行这类逻辑约束时大M的取值不能太大大M太大容易导致数值病态求解器在判可行性时会出问题大M太小又有可能把可行域错误地截断。我的经验是M取该设备全年最大运行功率对应的影响量级乘以1.5就够了。4. 并网 vs 离网容量配置结果与调度规律4.1 典型日的选取与验证典型日聚类完成后先做一步验证工作把聚类得到的四个典型日按权重叠加算出来的全年总发电量和全年总制氢量跟直接用8760小时逐时数据计算的结果做对比。偏差在5%以内说明这四个典型日确实能代表全年的资源特性。我复现时最终选出的四个典型日大致对应冬季大风日、春季过渡日、夏季高温高辐照日、秋季多云日。冬季大风日风速高但辐照弱光伏出力小夏季光伏出力大但风速相对较低。这种差异化组合正好能覆盖风光的互补特性也让优化器有足够的场景去权衡容量配比。4.2 两种模式下的最优容量配置对比求解完成后最直观的结果是并网和离网两种模式下系统最优容量配置的差异。为了对比统一我在两种模式下都设定相同的年氨产量目标。表并网与离网模式下最优容量配置对比设备并网模式离网模式风电装机容量兆瓦5876光伏装机容量兆瓦4266电解槽额定功率兆瓦5162储氢罐有效容积万标准立方米12.821.5储能电池容量兆瓦时无38并网模式下系统大幅依赖电网购电来补充短时功率缺额风电和光伏的装机容量相对较小储氢罐也不需要做得很大因为电网本身就是一个无限的缓冲池。离网模式下风电光伏装机容量分别提高了31%和57%储氢罐扩容68%这是因为离网系统必须靠源端冗余和储氢来应对风光出力的连续波动没有任何外部支撑。从经济性指标来看离网模式的单位氨生产成本比并网模式高约15%到20%。多出来的成本主要来自更大的初投资以及为了满足零外购电约束而必然出现的弃风弃光电量。这一结果符合预期但量化出来还是很有参考价值的如果项目选址附近有稳定电网接入条件优先选并网模式单位成本明显更低只有定位在零碳示范项目或者偏远无电网地区才值得考虑离网方案。4.3 调度策略与经济性分析调度结果的规律比容量结果更有意思。并网模式下电解槽的运行曲线明显跟着电价走谷电时段电解功率爬升到接近额定功率峰电时段则降功率甚至停机氢气需求缺口依靠储氢罐中预先存好的氢气来补足。这实际上就是利用氢储能的跨时段转移特性做电价套利。离网模式下没有电价信号可以跟随调度策略完全由风光出力和储氢罐状态决定。典型大风日的调度曲线显示风电出力高峰时段电解槽满负荷运行同时储能电池充电风电出力低谷时段储能电池放电维持电解槽基本负荷氨合成塔始终保持不低于70%的设计负荷运行。储氢罐在这个场景里扮演的角色就是能量缓冲池风大的时候囤氢风小的时候放氢保证合成氨连续生产。经济性分析方面把成本拆开看并网模式下运行成本中购电费用占比接近40%这是降本的关键突破口离网模式下投资成本占比超过75%运行成本主要是设备折旧摊销边际运行成本很低多产一吨氨的增量成本远低于并网模式。所以长期来看如果风光设备成本持续下降离网模式的经济性劣势会逐步收窄。5. 调参、踩坑与排查实录5.1 求解失败与可行性问题复现过程中最难受的问题不是求解慢而是模型不可行。Cplex返回infeasible的时候第一反应不要怀疑求解器先怀疑自己的约束建模。我遇到的一次典型问题是氨合成塔的最低运行负荷约束和储氢罐的初始储氢量约束冲突了。模型假设第一天8点开始运行储氢罐初始储氢量设得偏低而氨合成塔要求全天候最低负荷运行导致晚上氢气供应不足约束冲突。排查方法是用Cplex的冲突定理解析功能让求解器自动找出一组不可行约束子集很快就能定位到问题所在。cplex.conflict.minConflict(); cplex.conflict.get();这个功能在调试初期特别实用但需要提醒的是对大规模模型启用冲突检测会显著增加求解时间建议单独跑一个调试版的模型文件来做冲突分析不要在生产版本里开这个功能。5.2 数据量纲和精度问题的几个典型坑量纲问题是新手最容易踩的坑我也栽过一次。风速数据用的单位是米每秒功率用兆瓦储氢量用公斤而电价用的单位是元每千瓦时几个量纲混在一起约束矩阵的系数跨度从1e-3到1e5数值范围差异过大Cplex内部处理时容易出现数值警告。解决方案是统一量纲体系所有能量统一用兆瓦时氢量统一用吨成本统一用万元。功率平衡约束里所有项都变成兆瓦物料平衡约束里所有项都变成吨这之后系数范围就基本控制在1e-2到1e3之间了数值上的报警信息也消失了。另一个坑是储氢罐动态方程里的时间步长。我们用的都是1小时步长但电解槽制氢量的计算涉及连续时间到离散时间的转换。用1小时步长的时候某时段电解槽以50兆瓦功率运行对应的氢气产量是该时段总制氢量这没有歧义。但如果后期把模型改成15分钟调度就有必要核查效率系数是否需要按时间步长重新换算稍不留神就会导致全年的氢气平衡不闭合。5.3 性能瓶颈与后续扩展建议模型规模增大后求解时间会非线性上升。96个时段的基础模型Cplex大约在5分钟内跑到0.5%间隙如果把典型日扩展到12个求解时间会涨到40分钟级别。如果后续想把8760小时全部纳入模型做精细化分析建议采用Benders分解或者拉格朗日松弛把原问题拆成容量主问题和调度子问题用迭代的方式收敛而不是让Cplex硬扛。从模型扩展角度看几个方向值得考虑。第一把电解槽的动态响应特性加入约束即电解功率的爬坡速率限制这会增加约束数量但更贴近实际第二考虑氨合成塔的启停成本在目标函数里引入固定启动费用这样调度优化会在维持低负荷运行和停机再启动之间做更真实的权衡第三加入可再生能源出力不确定性的随机优化或鲁棒优化扩展这是目前学术研究的热点方向也是工程实际中必须面对的问题。我这次复现的全部代码和数据文件按照模块化方式组织每个模块都有详细的注释跑完一遍以后最大的感受是这类问题的难点从来不在求解器本身而在于把物理过程转化为数学约束时做的那些简化和取舍。每一个线性化近似、每一个典型日聚类、每一处效率系数的取值都会影响最终容量配置的经济性和可行性。建模的时候多想一步调试的时候就能省十步。