1. 项目概述与整体设计思路
1.1 项目背景与系统构成
前阵子复现了一个典型的可再生能源制氢-合成氨容量配置项目,主题是“风光互补制氢合成氨系统的容量-调度联合优化,用Matlab调Cplex求解”。这类系统这几年在新能源圈子里讨论很多,本质上解决的是风光出力波动与化工连续生产之间的匹配问题。风电、光伏出力天然不稳定,电解水制氢相当于把电力转化为便于存储的氢能,氢再与空气中分离出来的氮气通过合成氨工艺生成氨,氨既是化工原料,又是相对易储运的绿色燃料。项目代码同时考虑了并网和离网两种运行模式,并网模式下系统可以与电网买卖电作为缓冲,离网模式下则完全自平衡。
系统的抽象结构并不复杂:左侧是风电场和光伏电站,中间是电解槽、储氢罐、合成氨装置(配套空分制氮),右侧在并网模式下多一个电网交互节点。所有设备都有容量上限,所有能量流都必须满足时序平衡。所谓“容量-调度优化”,就是把两个层面的问题放进同一个优化框架:容量问题决定风电装机、光伏装机、电解槽规模、储氢罐规模、合成氨装置规模各建多大;调度问题决定在每一个小时(或典型日)风电怎么分配、电解槽开多少功率、储氢罐什么时刻充放、合成氨装置保持多高负荷、并网模式买电还是卖电。
这两个问题单拆开看都不难,难在互相耦合。装几台电解槽取决于全年氢需求的预期,而氢需求又由合成氨装置的调度决定,合成氨装置的负荷又反过来受风光出力和储氢状态影响。最稳妥的做法就是把容量变量和调度变量丢进一个统一的优化模型里,一次性交给求解器处理。这也是标题里“容量-调度优化”要表达的核心意思。
1.2 为什么必须做“容量-调度”联合优化
很多人拿到这类题目,第一反应是先算容量再算调度,两步走。但两步走的坑在于:容量配置的结果依赖于对未来运行方式的假设,这些假设一旦和真实最优调度不一致,容量投资就是错的。举个具体例子,如果容量优化阶段假设电解槽全年满发,风光配比会倾向多装风电保证氢量;但实际调度里夜间风光出力低,电解槽根本达不到满发,储氢罐又不够大,最后氢供应不足,合成氨装置被迫降负荷,项目收益被高估。
反过来,只做调度优化不做容量优化,给定一个固定装机容量,调度就算做得再漂亮,也无法回答“这个站一开始该建多大”的问题,而建多大恰恰是投资决策里最敏感的变量。所以项目里常见的做法是单层大模型:容量决策变量是标量,调度决策变量是时间序列,它们一起作为决策变量,目标函数统一考虑年化投资成本加年运行成本,约束条件同时覆盖投资约束和时序运行约束,交给Cplex求解。这种方式建模直观,求解在MILP框架下能保证全局最优,比启发式算法稳定得多。
1.3 求解器选型:为什么是Cplex
做这类优化,求解器选型基本等于半条命。项目用的是Cplex,这在能源系统优化领域算是长期验证过的选择。Cplex对大规模线性规划(LP)和混合整数线性规划(MILP)的数值稳定性好,分支定界效率高,内置的预处理和切割平面技术很成熟。相比之下,遗传算法、粒子群这类启发式算法虽然不需要建模太严格,但无法保证收敛到全局最优,而且每次运行结果可能不一样。对于容量配置这种动辄上亿元的决策,一个“差不多”的解和一个全局最优解之间的差距,可能就是几千万元的偏差。
Cplex在Matlab里的接入方式有两种:一种是通过YALMIP建模工具箱,另一种是直接用Cplex提供的Matlab接口函数。我更推荐YALMIP,因为它的语法接近数学表达式写法,约束可以通过向量化表达式一次性传入,调试起来比手拼稀疏矩阵舒服得多。后续第3节我会把环境配置和核心代码结构都列出来。
2. 数学模型:如何把系统“翻译”成优化问题
2.1 目标函数的构成与量纲
先看目标函数。模型里的成本主要分三块:投资成本年化、运行维护成本、购电成本(并网模式)。收入主要两块:卖电收入(并网模式)和氨销售收入。写成数学表达式就是:
目标 = 年化投资成本 + 年运行成本 + 年购电成本 − 年卖电收入 − 年氨销售收入
投资成本年化,要把设备初投资乘上一个年化系数CRF。CRF的公式是 r(1+r)^n / [(1+r)^n − 1],r是折现率,n是设备寿命。这个系数经常被初学者漏掉——如果不做年化,把30年的设备投资直接和1年的运行成本加在一起,量纲就错了,优化结果会拼命缩小设备容量。比如说折现率8%、寿命20年,CRF大约是0.10左右,也就是1000万元的设备投资,折算到每年约100万元。代码里应统一按“万元/年”或“元/年”口径处理,及时把各设备初投资折算成年值。
氨销售收入的结构要写清楚:氨产量设为变量Q_am,单位吨,价格price_nh3按元/吨。合成氨的氢耗换算要用化学计量关系:生产1吨氨约消耗0.18吨氢。这个换算系数在代码里写成一个常数,后面会接进储氢罐平衡方程。如果这个系数搞错,模型的物料平衡就会失真,优化结果看起来正常,实际上氢耗和产氨根本不匹配。
2.2 容量层决策变量与约束
容量层的决策变量主要是各设备的额定容量:风电装机、光伏装机、电解槽功率上限、储氢罐容量、合成氨装置产能上限,离网模式再加一个蓄电池容量(如果有)。每个设备通常有上下限约束,比如当地可建设面积限制了风电最大台数、光照资源限制了光伏容量,资金预算也可以作为约束加进去。这些边界条件的取值看似不起眼,实际上对结果影响很大,建议做敏感性分析时优先扫这些参数。
容量层还有一个容易被忽略的环节:设备额定容量和运行功率的换算关系。比如电解槽额定功率10MW,功率调节范围一般是20%—100%,也就是最低运行功率2MW。这意味着如果风电总出力低于2MW,电解槽要么停机,要么靠电网或储能把功率补足,离网模式下无法工作。这类“最小负载率”约束是模型可行性的关键,很多复现失败都是卡在这一条。
2.3 运行层约束:四类关键约束逐一拆解
运行层约束是整个模型里最耗费精力的部分,按约束类别逐一说明。
第一类是功率平衡约束。系统在任何时刻都必须满足:风电出力 + 光伏出力 + 购电(并网)= 电解槽用电 + 合成氨装置用电 + 售电(并网)+ 蓄电池充放电修正(离网)。这里所有功率必须先归一化到“系统母线”口径,风电出力并不是全额给电解槽,有时候还要考虑主动弃风弃光,所以在变量里会保留限电空间。
第二类是电解槽运行约束。除了功率上下限,还要考虑爬坡约束和开关状态。如果引入二进制变量表示电解槽开停机,最小负载率的表达就用得上:P_el >= P_min * b_el,P_el <= P_max * b_el。引入整数变量的代价是模型变成MILP,求解难度上升,但换来的是更真实的运行逻辑。如果调度粒度是小时级,开关变量数量等于T乘设备数,T是全年8760小时的话,变量规模会很大,通常要用典型日压缩。
第三类是储氢罐时序约束。储氢量S_h2(t)的状态转移可以写成:S_h2(t+1) = S_h2(t) + eta_el * P_el(t) * dt − h2_am(t),其中eta_el是电解制氢转换效率,h2_am(t)是合成氨装置消耗氢气量,dt是调度时段长度,单位小时。注意储氢罐不能是空的,合成氨装置要连续运行,所以一般还要设置最低储量约束。另一个关键约束是终值约束:S_h2(end) == S_h2(1),这是典型日优化里的周期性条件,代表罐子不能凭空多出氢气,也不能把初始储存的白拿光,不然下一轮循环就失衡了。
第四类是合成氨装置运行约束。合成氨装置有最低负荷率限制,也有最大产能限制,单位时间产氨量决定了氢耗速率。在简化模型里,合成氨装置可以看成一个“功率-氢耗-产氨”耦合的模块,通过线性关系把电耗、氢耗和氨产量联动起来。氮气基本默认充足,不太需要关注空分系统;如果想更贴近工程,可以再加一个空分电耗常数,进功率平衡方程。
2.4 并网与离网模式:模型差异在哪儿
并网和离网不是简单删一个变量的问题。并网模式下,系统有电网这个“无限池子”,缺电就买、富电就卖,优化结果往往倾向于小储氢罐、大电解槽,因为电价低时可以直接从电网取电制氢,不需要为了应付天生波动而多投资储能设备。但并网模式有一个重要约束:同一时刻不能同时购电和售电,否则模型会在买卖价差之间做无意义的套利,计算出虚假收益。代码里可以用一个二进制变量b_grid实现:购电量 <= M * b_grid,售电量 <= M * (1 − b_grid)。
离网模式则要严格满足自平衡,功率平衡方程里没有电网项,风光出力必须实时等于负荷,唯一能跨时段搬移能量的是储氢罐(如果有蓄电池则也包括它)。离网模式下模型更容易不可行,原因很简单:连续无风的几个晚上,光伏不出力,风电又不够,电解槽被迫停机,合成氨装置还得稳定供氢,储氢罐里的氢可能撑不到天亮。如果模型算出来不可行,先别急着怀疑求解器,大概率是某天的风光出力低谷和电解槽最小负载率约束撞车了。
2.5 时间序列简化:典型日的选取与权重
全年8760小时建模,变量规模直接爆炸,尤其是二进制变量。常见工程化做法是先做典型日聚类:把风速、光照、温度、电价序列聚成若干个典型日,每个典型日按出现频次赋予权重,优化目标也按权重折算到全年。聚类方法用k-means就够,但要注意把风速和光照放在同一个样本里做标准化,不要单独聚类。聚多少个典型日,视精度要求而定,一般4到12个就够用。
我见过有人把8760小时全部塞进去,结果Cplex跑了两天还在0.5%的gap里打转,最后只能靠加求解时间上限强行截断,得到的解根本没有说服力。用典型日是为了快速得到稳定解,前提是论文或报告里要交代清楚聚类方法、聚类个数以及误差分析。这一段内容在复现项目里非常重要,直接决定最终模型的可解性。聚类完之后,每个典型日的权重乘以该日运行成本,最后累加就是全年的运行成本口径。
3. Matlab+Cplex实战:环境配置、代码结构与调参
3.1 环境配置:YALMIP和Cplex的版本匹配
先说环境配置。Cplex要能被YALMIP正常调用,核心是把Cplex的Matlab路径加进来,并且确保Cplex版本和Matlab版本兼容。比如Matlab R2023a配Cplex 12.10/12.12一般没问题,Matlab太新而Cplex太旧时,会因为缺少对应的动态库报“Unable to load mex file”。个人建议:优先用yalmiptest命令诊断,如果Cplex的状态是found,就说明路径加载成功。
YALMIP本身只是个建模层,不参与求解,但它要求Cplex的求解器文件在搜索路径里。配置完成后,在Matlab命令窗口跑一句:
yalmiptest输出里看到CPLEX-IBM并且status是found,才算配置完成。如果没找到,多半是路径没加对,或者addpath之后没保存path,换个工作区就丢了。建议把路径写进startup.m文件,一次性解决每次打开都要手动加路径的麻烦。
3.2 决策变量定义与目标函数的代码写法
项目里建议先定义容量变量(标量),再定义运行变量(向量),全部用sdpvar,二进制变量用binvar。一个典型的定义片段如下:
T = 24; % 单日调度时段 % 容量决策变量 C_wind = sdpvar(1,1); C_pv = sdpvar(1,1); C_el = sdpvar(1,1); C_h2 = sdpvar(1,1); % 储氢罐容量 C_am = sdpvar(1,1); % 合成氨产能上限 % 运行决策变量 P_wind = sdpvar(T,1); % 风电消纳功率 P_pv = sdpvar(T,1); P_el = sdpvar(T,1); % 电解槽输入功率 P_am = sdpvar(T,1); % 合成氨装置电耗功率 P_buy = sdpvar(T,1); % 购电功率, 并网模式 P_sell = sdpvar(T,1); % 售电功率, 并网模式 S_h2 = sdpvar(T,1); % 储氢罐储氢量 b_el = binvar(T,1); % 电解槽开关状态 b_grid = binvar(T,1); % 购售电互斥标志, 并网模式目标函数写成:
CRF = 0.08; % 按折现率和设备寿命计算 IC = CRF * (inv_wind*C_wind + inv_pv*C_pv + inv_el*C_el + inv_h2*C_h2 + inv_am*C_am); OC = om_ratio * (inv_wind*C_wind + inv_pv*C_pv + inv_el*C_el + inv_h2*C_h2 + inv_am*C_am); EC = sum(price_buy .* P_buy * dt) - sum(price_sell .* P_sell * dt); Rev = price_nh3 * Q_am; % Q_am由合成氨产量累加得到 obj = IC + OC + EC - Rev;这里Q_am不能凭空定义,需要从合成氨装置的产氨约束推出来。一般做法是把单位时间产氨量设计成变量Q_am_t(T,1),乘时段数累加得到年产量,再由化学计量关系把氢耗量H2_need_t送进储氢罐平衡方程。这样收入项和耗氢项就挂上了钩。做这一步时,建议一开始就算好单位换算,避免后面再对“元/MWh”、“元/吨”、“MW”三套单位体系重新对齐。
3.3 约束条件写入:矩阵化比for循环快得多
约束的写入位置在目标函数之后、optimize之前。Matlab新手常见做法是用for循环把每条约束逐条加进F = [F, expr];,这在约束数量少的时候没问题,但到了几百上千条约束,反复append会显著拖慢建模速度。更好的做法是尽量用向量表达式一次写完整段约束。
功率平衡约束:
F = [F, P_wind + P_pv + P_buy == P_el + P_am + P_sell];风电、光伏出力不超过可用资源:
F = [F, 0 <= P_wind <= avail_wind .* C_wind]; F = [F, 0 <= P_pv <= avail_pv .* C_pv];电解槽最小负载率与开关联立:
F = [F, P_el >= r_min * C_el .* b_el]; F = [F, P_el <= r_max * C_el .* b_el];储氢罐状态转移与容量限制:
F = [F, S_h2(2:end) == S_h2(1:end-1) + eta_el * P_el(1:end-1) * dt - H2_need(1:end-1)]; F = [F, S_h2 >= S_h2_min]; F = [F, S_h2 <= C_h2]; F = [F, S_h2(1) == S_h2(end)]; % 周期性约束并网模式的购售电互斥:
M = 1000; % 足够大, 但不要离谱 F = [F, P_buy <= M .* b_grid]; F = [F, P_sell <= M .* (1 - b_grid)];这里的M是经典坑位。M取1e6,模型数值条件会非常差,求解器容易报numerical trouble;M取小了,又会误伤正常解。实际项目中先用一个粗略的上限,比如最大负荷的1.5倍,够用即可,求解结果出来后再检查有没有触碰到M边界。经验数据是M不要超过问题正常量级的3倍。如果嫌求M烦,可以把P_buy和P_sell改为共享一个总交换功率上限的写法,省掉互斥二进制变量,代价是可能允许同买同卖的小额操作,但对大多数规划性研究,这个近似误差可以接受。
3.4 求解参数设置:mipgap、时间限制与诊断输出
最后调用optimize时,求解参数直接传进去:
ops = sdpsettings('solver','cplex', ... 'verbose',2, ... 'cplex.mip.tolerances.mipgap',0.01, ... 'cplex.timelimit',3600); result = optimize(F, obj, ops);mipgap设成1%,意思是求解器只要找到的整数解和最优下界之间相对差距小于1%就停止。这是大模型里非常实用的设置——死磕到0%的gap,瓶颈和收益往往不成比例。timelimit控制最大求解时间,避免模型卡在难以收敛的分支定界过程。verbose等级决定输出多细,调试阶段设2,跑大模型设0,减少控制台刷屏。
求解完之后,第一时间检查result.problem。0表示求解成功,1表示不可行,2表示数值问题,4表示求解中断。YALMIP的check(F)可以逐条检查约束的最大违反量,定位是哪一条约束出了问题。这些诊断方法是复现路上最值得掌握的工具。
4. 复现避坑指南:常见问题与排查方法
4.1 Cplex找不到或mex文件加载失败
这个问题排在首位,因为环境都跑不通,后面全白费。常见原因有三个:一是Cplex的Matlab接口路径不在搜索路径里;二是Matlab版本和Cplex版本不兼容;三是Windows/Linux下Cplex安装路径含中文或空格导致动态库加载失败。排查思路很简单:先执行yalmiptest看状态,再检查addpath是否正确。如果还不行,打开Cplex安装目录下的cplex/matlab/文件夹,看有没有对应平台的子目录,确认这个子目录在Matlab路径里。个人建议装Cplex时直接用默认安装路径,别放中文目录。
4.2 离网模式一求解就不可行,问题出在哪里
离网模式不可行,是这类项目里最高频的bug。多数情况下并不是模型写错了,而是物理条件本身就不允许。排查建议按顺序来:第一,看功率平衡。把风力光伏的可用出力乘容量序列调出来,逐时段算最大可发电量,再看电解槽最小负载率和合成氨装置基础电耗的总和,如果某个时段可发电量小于最小必须功率,那这个场景一定无解。处理办法有两种:配置储能平滑功率,或者允许一定比例切负荷量,在方案里通常表述为失负荷率约束。
第二,看储氢罐周期性约束。终值等于初值的要求很严格,如果初始储量设太低、合成氨装置又必须维持连续生产,模型可能无论如何都找不到周期内的可行状态转移。变通做法是把终值约束放宽成S_h2(end) >= S_h2_min_period,或者把初值设成和终值一致的未知变量,让求解器自己寻找平衡点。
第三,看电解槽最小负载率。前面讲过,P_el >= r_min * C_el * b_el这条约束在某些低出力时段会逼迫模型要么把电解槽开到最低功率、要么停机。如果忽略启停费用,模型可能频繁启停,工程上不现实但数学上可行。如果发现“可行但结果诡异”,先检查是不是缺了启停成本惩罚或最小连续运行时间约束。
4.3 MIP求解太慢,卡在1%的gap不动弹
Cplex对MILP的求解速度总体不错,但模型一旦上了几千个二进制变量,还是会慢。缓解手段按优先级排:第一选择是压缩时间序列,用典型日聚类代替8760小时;第二选择是减少二进制变量——电解槽开关变量如果对结果不敏感,可以换成功率下限约束的连续化近似;第三选择是设置mipgap和timelimit,接受次优但工程上够用的解;第四选择是warm start,用启发式解或上一轮迭代的解作为初始解喂给Cplex,大幅缩减分支定界搜索空间。
还有一点很少有人提:Cplex预处理阶段会翻来覆去检查约束,如果模型里冗余约束过多,预处理时间也会很感人。把明显重复的约束删掉,比如两个等式方程可以互相推出的那种,能让求解速度明显提升。实际调试时,可以先用一个小规模样例验证模型正确性,再放大到全规模,否则一边查正确性一边等求解,效率极低。
4.4 优化结果反常识,问题多半在目标函数或单位
复现完第一版模型,我遇到过“最优结果是啥也不建”,氨产量为零,设备全不投资。第一反应是模型错了,仔细排查后没发现问题,最后发现是氨销售收入价格给低了,导致卖氨收入覆盖不了投资和运行成本,模型理性选择就是不做。所以看到反常识结果,先别急着怀疑求解器,去检查经济参数:氨价、电价、设备造价、折现率,任何一个偏差都可能彻底改变最优配置。这类系统的经济性测算结果存在真实边界条件,参数差10%结论可能完全翻转。
另一个容易踩的坑是单位。风电出力在模型里通常用“占额定容量的比例系数”表示,avail_wind取值0到1之间,再乘容量才能得到MW。如果avail_wind误填成kW/ MW混合,结果直接废掉。建议建模最开始就把单位约定写清楚,所有变量统一用MW、MWh、吨、元三套基础单位,并且互算关系单独写在代码注释里。
4.5 结果分析:怎么从优化解里提炼有用结论
代码跑通拿到最优解之后,事情并没有结束。容量层面要给出各设备容量配置,运行层面要给出典型日调度曲线。我最常用的分析手段:一是画出典型日内风电、光伏、电解槽功率、储氢量、氨产量的时序曲线,观察风光资源互补性在什么时段体现;二是做敏感性分析,把氨价、碳价、设备造价在±30%范围扫一遍,看最优容量配置如何移动,这种图往往是论文或报告里最有力的证据;三是并网和离网两种模式综合对比,量化电网为系统带来的灵活性价值。通常结论是并网模式允许更小储氢罐和更大电解槽,项目总投资更低。
4.6 扩展方向:模型还能往哪些方向改
复现完这套代码,你可能还想加现实因素进去,常见扩展至少三个方向。第一个方向是加入碳交易机制或绿氨溢价,在目标函数中增加碳排放成本项,让系统贴近双碳政策下的真实经济环境。第二个方向是加设备退化模型,电解槽频繁启停影响使用寿命,在目标函数中加入启停成本或寿命折算,调度结果会明显更“温和”。第三个方向是引入需求响应或现货电价曲线,并网模式下购电策略就有了博弈空间,储氢罐可以和电网交互联动起来,变成一个真正意义上的能源枢纽调度模型。无论往哪个方向加,核心逻辑还是这套“容量-调度”联合优化的骨架,公式、代码、求解器基本都能复用。
我个人在实际复现过程中的体会是,这类项目最难的不是把代码跑通,而是把每个参数背后的物理意义和工程假设搞清楚。Cplex只在数学上保证最优,但模型本身是否真实反映系统运行规律,完全取决于建模的人。建议拿到代码后先跑通一个典型日的小规模案例,确认各类约束行为符合预期,再逐步放大时间尺度和设备种类。这样既能快速定位问题,也不会在错误模型上浪费大量求解时间。