过去两年我一直在跟新能源制氢方向的项目打交道,最让我头疼的不是电解槽本身,而是它背后那一整套“风光怎么配、氢怎么储、氨怎么合成、电怎么调度”的系统级问题。最近正好复现了一篇关于并离网风光互补制氢合成氨系统的容量与调度优化论文,用的是Matlab + Cplex求解,折腾了将近三周才把整个模型跑通。这篇文章我就把复现过程里最核心的建模思路、求解实现和踩过的坑一次性讲清楚。
这篇内容适合三类人:一是正在做新能源制氢、绿色化工方向课题的研究生,二是想了解风光互补系统容量配置与日前调度如何联合优化的工程师,三是拿到相关代码但看不懂模型逻辑、想彻底搞懂每一步的人。我会从系统架构、数学模型、Cplex求解实现、再到复现过程的高频问题,尽量用最直白的方式把这件事讲明白。
1. 并离网风光互补制氢合成氨:复现之前,先把物理图景还原清楚
复现代码最忌讳的事就是一上来就盯公式。如果你不清楚这套系统在物理世界到底长什么样、各个设备之间怎么连接、能流的走向是什么,那后面无论是建模还是调参都会一头雾水。所以第一步我花了整整两天把这篇论文的物理拓扑还原出来。
1.1 风光互补到底“互补”在什么维度
这个系统的能量源头是风电场和光伏阵列。单独看风电,它的出力特点是夜间大、白天小,冬季大、夏季小,具有很强的季节性和“反调峰”特性;单独看光伏,则是白天有、晚上没有,夏季充沛、冬季匮乏。两者叠加之后,在日尺度上可以实现“白天靠光、夜间靠风”的接力,在季节尺度上也能部分对冲彼此的高峰和低谷。
但这里有一个非常关键的认知:风光互补不是简单把两者装机加在一起就行,而是需要在容量配置层面去权衡。比如某地区风速资源好但光照一般,那你多配风电、少配光伏才是经济上最优的;反之亦然。论文里的容量优化部分,本质上就是在“给定风光资源曲线、设备成本、电价机制”的约束下,找出一组最优的风电装机、光伏装机、电解槽规模、储氢罐大小、蓄电池容量,让整个系统的年化总成本最低。
1.2 系统拓扑:从风、光、氢、氨到电网的全链路能流关系
我复现的这套系统,完整的能流链路是这样的(我画过不下十遍,最后还是用文字描述最清楚):
- 风电场和光伏阵列发出的直流/交流电,经过变流器汇入直流母线或交流母线
- 电力优先供给电解槽生产氢气,多余的电给蓄电池充电,再多才考虑上网卖电
- 电解槽产出的氢气进入储氢罐,而不是直接送去合成氨
- 储氢罐的氢气按需送入合成氨装置,与氮气在高温高压下反应生成氨
- 如果系统并网运行,母线还通过变压器和电网相连,缺电时可以买电,富电时可以售电
- 如果系统离网运行,电网这条通路就完全切断,系统必须靠自身风光出力和储能装置维持能量平衡
这里面有一个很多复现者容易忽略的点:合成氨装置和电解槽在运行特性上有本质区别。电解槽适合波动性供电,功率可以从零平滑升到额定值;但合成氨回路一旦停车再启动,需要很长的升温还原时间,而且频繁启停会严重缩短催化剂寿命。所以在调度模型里,合成氨装置通常被建模成“要么不运行,要么以某一较高负荷率运行”,不能像电解槽那样随意调节。论文里把这个约束通过“最小运行状态”或“最小出力比例”的方式表达出来,复现的时候一定要留意。
1.3 并网和离网,不只是差了一条电网连接线
很多人把并离网理解成“有没有电网接入”,但在优化模型里,这两者的差别是结构性的:
并网模式下,功率平衡约束多了一个与电网交互的变量$P_{grid,t}$,它既可以是正数(买电)也可以是负数(售电)。同时,目标函数里会多出“购电成本”和“售电收益”两项。这相当于给系统提供了一个无限容量的“虚拟储能”,理论上只要电价合适,电网就可以平衡风光出力的任何缺口或盈余。
离网模式下,功率平衡约束里没有电网项,系统只能靠蓄电池和储氢罐来平抑波动。这时有两个硬性约束变得极其重要:一是全年任意时刻都不能出现功率缺额,否则系统无法运行;二是蓄电池的容量必须足够大,才能扛过“连续无风无光”的极端时段。从数学上看,离网系统的可行域比并网系统小很多,最优容量配置通常也偏“冗余”,这也是为什么离网制氢的综合成本往往高于并网制氢。
理解了这些物理图景,后面建模才有根基。我强烈建议所有复现者,先把论文里的系统示意图用自己的方式重新画一遍,把每个设备的输入、输出、中间缓冲都标清楚,再做建模。
2. 容量与调度的耦合逻辑:两个优化层级到底在优化什么
标题里“容量-调度优化”这六个字是全文的灵魂。容量优化回答的是“建多大”的问题,调度优化回答的是“怎么用”的问题。但这两个问题绝不是独立的,它们之间存在强烈的耦合关系。
2.1 容量层:决定设备装机规模的上层决策
容量优化的决策变量包括:风电装机容量$P_{wt}^{cap}$、光伏装机容量$P_{pv}^{cap}$、电解槽额定功率$P_{el}^{cap}$、储氢罐容量$V_{H2}^{cap}$、蓄电池额定容量$E_{bat}^{cap}$等。
这些变量对应的成本包括初始投资成本和固定运维成本。论文里通常会把它们折算成年值——比如把一次性投资乘以一个资本回收系数(CRF),再叠加每年的固定运维费,得到“年化投资成本”。
容量变量是全局性的,一年四季都不变,但它们的取值直接决定了下层调度是否有足够的“物理空间”来完成功率平衡和物料平衡。
2.2 调度层:在每个时段内决定设备的运行状态
调度优化解决的是“在已知容量配置和风光出力曲线的前提下,每小时该干什么”。常见的调度时间尺度是小时级,一个调度周期可以是典型日、典型周,或者直接是一整年8760小时。
调度决策变量包括:每小时电解槽的输入功率$P_{el,t}$、电解槽的启停状态$u_{el,t}$、每小时储氢罐的充放氢量、每小时合成氨装置的运行状态与用氢量、每小时蓄电池的充放电功率$P_{dis,t}$和$P_{ch,t}$以及SOC变化、并网模式下每小时的购售电功率等。
如果调度周期是全年8760小时,那决策变量的数量会非常庞大。以一个典型系统为例:假设有10类主要决策变量,全年8760小时就是超过8万个决策变量,再加上约束条件里的连续变量和0-1变量混合,模型的规模很快就会突破十万级。这就是为什么需要用Cplex这种商业求解器,而不是自己写启发式算法。
2.3 为什么必须把容量和调度放在同一个优化框架里
常见的简化做法是两步走:先假设一些典型运行场景,完成容量优化;再固定容量,做调度优化。但这样做有个致命缺陷——容量配置和调度策略是互相影响的。举个例子,如果储氢罐容量配得特别小,那么调度层就只能减少电解槽的运行时长,氢气产量下降,氨产量跟着下降,收益减少;如果储氢罐容量配得大,虽然投资成本上去了,但调度灵活性更高、氨产量更稳定。两步走的做法往往得到的是一个局部最优,而不是全局最优。
论文里采用的是一体化建模思路:把容量变量和调度变量放在同一个优化问题里,用混合整数线性规划(MILP)统一求解。这意味着容量变量一旦确定,调度策略会自动适配;调度策略又会反过来影响容量变量的最优取值,两者在求解器内部同时迭代收敛到全局最优解。这也是这篇论文最大的技术亮点。
3. 数学模型拆解:目标函数、约束条件与关键表达式的设计意图
接下来的部分是整篇博文的核心,也是复现过程中最需要仔细推敲的部分。我会按照“目标函数-约束条件”的逻辑依次展开,每讲一个公式都会解释为什么这么建模。
3.1 目标函数:综合年化成本最小化
大多数相关论文的目标函数都写成“年化总成本最小”,但也有些论文会加入碳排放惩罚项。我复现的这篇采用的是年化总成本,包含以下几部分:
- 设备年化投资成本:把风电、光伏、电解槽、储氢罐、蓄电池、合成氨装置的初始投资总额乘以资本回收系数
- 年运行维护成本:各设备按额定容量的某个百分比计取
- 购电成本与售电收益:仅并网模式存在
- 制氨收益:按氨产量乘以氨的市场价格,作为负成本计入目标函数
写成数学形式,目标函数是:
$$\min ; C_{inv} + C_{om} + \sum_{t}(c_{buy}P_{grid,t}^{+} - c_{sell}P_{grid,t}^{-}) - p_{NH_3} \cdot M_{NH_3,total}$$
注意这里有几处容易踩坑的地方:
第一,资本回收系数的计算公式是$CRF = \frac{r(1+r)^n}{(1+r)^n - 1}$,其中$r$是贴现率,$n$是设备寿命。不同设备的寿命可能不同,比如光伏25年、电解槽10年,复现的时候需要分别计算再累加,不能图省事用一个统一寿命。
第二,购电价格和售电价格往往不一样,而且很多地区实行分时电价。论文里给出的电价曲线可能是一天24小时的固定序列,复现时要用插值或映射的方式扩展到全年8760小时。
第三,氨的价格波动较大,论文里通常给的是某个基准价格,做灵敏度分析时才会变化。复现时先用基准价格跑通,再做价格敏感性测试。
3.2 功率平衡约束:并网与离网模式的统一表达
功率平衡是整个模型最核心的约束。在统一框架下,它写成:
$$P_{wt,t} + P_{pv,t} + P_{dis,t} - P_{ch,t} + P_{grid,t}^{buy} - P_{grid,t}^{sell} = P_{el,t} + P_{load,aux,t}$$
其中$P_{load,aux,t}$代表辅助负荷,包括合成氨装置的电力消耗、压缩机的电耗、控制系统用电等。辅助负荷通常被简化成一个固定值或者与合成氨产量成比例的变量,但在精细建模时可以拆解成多项。
离网模式下,$P_{grid,t}^{buy}$和$P_{grid,t}^{sell}$都强制为0。这里要注意的是,如果没有蓄电池,功率平衡在任何时刻都必须严格相等;有了蓄电池之后,等式右侧多了一项蓄电池的净充电功率,而蓄电池的SOC变化又和时间耦合在一起,这就把一个“单时刻平衡问题”变成了“跨时间耦合问题”。
我最初复现时用了一个愚蠢的写法:把蓄电池的SOC当成独立变量,没写SOC递推方程,结果模型算出来的结果非常离谱——SOC变成了无源之水。正确写法必须包含:
$$SOC_{t+1} = SOC_t + \eta_{ch}P_{ch,t}\Delta t - \frac{P_{dis,t}\Delta t}{\eta_{dis}}$$
以及边界约束$SOC_{min} \le SOC_t \le SOC_{max}$和首尾相接约束$SOC_0 = SOC_T$。首尾相接这步特别重要,它保证调度策略在循环运行时不消耗“电池初始电量”这个额外能量来源。
3.3 电解槽运行约束:比你想的复杂得多
电解槽不是简单的一个可变负载。真实碱性电解槽的运行受到以下限制:输入功率必须在最小技术出力$P_{el}^{min}$和额定功率$P_{el}^{cap}$之间;启动时需要消耗额外的能量进行预热;爬坡速率有限,不能瞬间从0冲到满功率;频繁启停会降低寿命,所以模型往往还会限制一天内启停次数。
在MILP框架下,电解槽约束通常引入一个0-1状态变量$u_{el,t}$,表示第t小时是否运行。核心约束为:
$$u_{el,t}P_{el}^{min} \le P_{el,t} \le u_{el,t}P_{el}^{cap}$$
这个双线性约束看起来简单,但它在Cplex里是作为一个“变量上下限绑定”的线性约束直接处理的。意思是:当$u=0$时,$P_{el,t}$被钳制在0;当$u=1$时,$P_{el,t}$可以在上下限之间自由变动。
如果论文还考虑了冷启动/热启动,就需要引入额外的状态变量$v_{el,t}$表示“第t小时是否发生了启动动作”,并把$v_t \ge u_t - u_{t-1}$这样的逻辑关系写进去。启动能耗可以建模成一个固定的能量消耗项$E_{start} \cdot v_t$附加在功率平衡的负荷侧。启停次数限制则用$\sum_{t}v_t \le N_{start,max}$表达。
3.4 储氢罐与合成氨装置的物料平衡约束
氢气的物料平衡是连接电力系统和化工系统的桥梁。每小时电解槽的产氢量、储氢罐的充放流量、合成氨装置的耗氢量,三者之间必须满足守恒关系:
$$V_{H2,t+1} = V_{H2,t} + \eta_{H2}P_{el,t}\Delta t \cdot \alpha - Q_{NH_3,t}$$
这里面$\alpha$是电解槽产氢率系数,$\eta_{H2}$是氢气压缩和储运效率,$Q_{NH_3,t}$是第t小时送往合成氨装置的氢气量。储氢罐同样有上下限约束和首尾相接约束。
合成氨装置的建模则是整个系统里最考验化工知识的部分。氨合成的反应式是$N_2 + 3H_2 \rightarrow 2NH_3$,也就是说每生成1吨氨需要大约0.176吨氢气。这个化学计量关系是硬约束,钢材无法突破。另外,合成氨回路有一个特点——单程转化率低,需要通过循环气把未反应的原料气送回反应器,这也意味着实际用氢量会比理论化学计量稍高一些,但大多数系统级优化论文会忽略这个细差,直接用理论计量比。
合成氨装置的最小运行负荷通常不为0,比如不能低于额定负荷的50%,否则反应温度维持不住。同样需要引入0-1变量$u_{NH_3,t}$和约束$u_{NH_3,t}Q_{NH_3}^{min} \le Q_{NH_3,t} \le u_{NH_3,t}Q_{NH_3}^{cap}$。
3.5 风光出力约束与资源数据预处理
风光出力不能超过当地资源条件下的最大可发功率。也就是说$0 \le P_{wt,t} \le P_{wt,t}^{available}$,且$P_{wt,t}^{available}$由风速曲线和风机功率曲线决定,$P_{pv,t}^{available}$由光照辐射强度和环境温度决定。
复现中最容易忽略的一步是:从原始气象数据到风机/光伏出力的转换。如果你直接用论文附录里给出的“归一化出力曲线”,那问题不大;但如果你自己用的是NASA或当地气象站的数据,就必须做两步处理:第一步把风速数据通过风机功率曲线转化为出力系数,第二步把光照辐射数据通过光伏效率模型转化为出力系数,然后再乘以装机容量得到实际出力上限。我见过不少人直接把风速除以额定风速当作风电出力系数,这会把大量高于额定风速的时段全部算成满发,结果容量优化严重偏向减少风电装机,模型结果完全失真。
4. Cplex求解的Matlab实现:从Yalmip建模到求解器配置
理论模型搭建完毕之后,就到了真正动手写代码的阶段。我使用的是Matlab R2022a + Yalmip工具箱 + IBM Cplex求解器的组合。这个组合在学术圈非常成熟,Yalmip负责把高层的优化模型自动转换成Cplex能理解的标准格式,省去了手写矩阵的大把时间。
4.1 环境准备与工具箱安装
需要准备三样东西:Matlab本体、Yalmip工具箱、Cplex for Matlab接口。Yalmip的安装最简单——去GitHub上下载最新版,把文件夹加入Matlab路径即可。Cplex则需要先安装IBM ILOG Cplex Optimization Studio,然后在Matlab里配置路径。
配置完成后,在命令行输入yalmiptest,如果能看到一系列...passed的提示,就说明环境没问题。这里有个细节:Cplex的Matlab接口版本必须和Matlab版本兼容,比如Cplex 12.10支持R2021a及以后版本,太老的Cplex版本在新版Matlab上会直接加载失败。
4.2 模型参数与决策变量的定义方式
参数定义阶段,我把所有系统参数集中在结构体变量里,方便统一管理和后续调试。比如:
% 系统参数定义(示例节选) para.dt = 1; % 调度时段步长(小时) para.horizon = 8760; % 调度周期(全年小时数) % 风电 para.wt_capex = 4500; % 风电单位投资成本,元/kW para.wt_om = 0.02; % 风电年运维成本比例 para.wt_life = 20; % 风电寿命,年 % 电解槽 para.el_capex = 2500; % 电解槽单位投资成本,元/kW para.el_min_ratio = 0.1; % 电解槽最小技术出力比例 para.el_eff = 0.60; % 电解槽系统效率(含整流) para.el_life = 10; % 电解槽寿命,年决策变量的定义使用Yalmip的sdpvar和binvar函数:
% 容量决策变量(连续) C_wt = sdpvar(1, 1, 'full'); % 风电装机容量 C_pv = sdpvar(1, 1, 'full'); % 光伏装机容量 C_el = sdpvar(1, 1, 'full'); % 电解槽额定功率 C_h2 = sdpvar(1, 1, 'full'); % 储氢罐容量 C_bat = sdpvar(1, 1, 'full'); % 蓄电池容量 % 调度决策变量(连续) P_el = sdpvar(horizon, 1); % 电解槽逐小时输入功率 SOC = sdpvar(horizon + 1, 1); % 蓄电池SOC在horizon+1个时刻的值 V_h2 = sdpvar(horizon + 1, 1); % 储氢罐氢量 % 调度决策变量(0-1) U_el = binvar(horizon, 1); % 电解槽启停状态 U_nh3 = binvar(horizon, 1); % 合成氨装置运行状态Yalmip的厉害之处在于支持这样直接声明带时间索引的变量向量,后续写约束时完全用变量名而非矩阵索引,可读性极高。
4.3 约束建模与目标函数写入
以功率平衡约束为例,在Yalmip里写起来非常直观:
% 电网交互变量(并网模式) P_buy = sdpvar(horizon, 1); P_sell = sdpvar(horizon, 1); % 功率平衡约束 Constraints = []; Constraints = [Constraints, P_wt + P_pv + P_dis - P_ch + P_buy - P_sell == ... P_el + P_load_aux]; % 购售电互斥约束(可选,避免同时买和卖) Constraints = [Constraints, P_buy >= 0, P_sell >= 0]; Constraints = [Constraints, P_buy <= M*U_grid, P_sell <= M*(1 - U_grid)];最后那个互斥约束用了大M法。虽然并网模式下同时买电和卖电在数学上不一定最优(因为买卖价差),但不加约束可能造成Cplex在边界处给出同时买卖的解,物理上不合理。加了互斥后模型规模会大一些,但求解更干净。
目标函数在Yalmip里用一个sum表达式就搞定了:
Objective = C_inv + C_om + sum(c_buy .* P_buy) - sum(c_sell .* P_sell) ... - price_nh3 * sum(Q_nh3 * m_factor);其中C_inv和C_om是容量变量的线性组合,写成cost_wt * C_wt + cost_pv * C_pv + ...的形式。
4.4 调用Cplex求解器:参数设置与效率优化
配置求解器时,我用的是optimize命令加sdpsettings:
ops = sdpsettings('solver', 'cplex', 'verbose', 2); ops.cplex.mip.tolerances.mipgap = 0.01; % 停止间隙 ops.cplex.timelimit = 7200; % 最长求解时间 ops.cplex.mip.strategy.search = 2; % 分支策略 ops.cplex.mip.limits.nodes = 5e6; % 分支节点上限 result = optimize(Constraints, Objective, ops);这里几个参数必须说一下:
mipgap决定求解精度。论文复现通常设置在1%~5%就足够了,设太小会让求解时间呈指数级增长timelimit是硬性保护。由于MILP求解的最坏时间复杂度是指数级的,不给时间上限的话可能好几天都跑不完- 如果不设置任何参数直接用Cplex默认设置,很可能遇到“内存爆炸”问题。Cplex默认会用完所有可用内存来存分支树,对于10万级变量的模型来说非常危险
4.5 并网与离网模式的代码切换设计
我把并网和离网做成了同一个模型的两种配置,通过一个mode标志位切换:
if strcmp(mode, 'grid') % 电网交互变量激活 buy_var = P_buy; sell_var = P_sell; else % 离网模式:强制买卖功率为0 buy_var = zeros(horizon, 1); sell_var = zeros(horizon, 1); Constraints = [Constraints, P_buy == 0, P_sell == 0]; end这样做的好处是,容量优化的结果可以直观对比两种模式的投资差异。我跑出来的典型结果是:在同样的风光资源条件下,离网系统需要的光伏和风电装机比并网多15%~25%,蓄电池容量更是要翻倍以上,才能保证全年供电的可靠性。这也解释了为什么当前大规模绿氢项目普遍优先选择并网制氢。
5. 复现过程中最容易踩的坑与排查思路
这一部分是全篇博文的精华。我把自己在复现过程中踩过的坑、浪费过的时间和排查思路完整记录下来,希望后来者能少走弯路。
5.1 问题一:模型求解时间过长,甚至跑不收敛
现象:第一次跑全年8760小时、变量数量超过10万个的模型,Cplex跑了4小时都没给出一个可行解。
根因分析:问题出在我把全年小时级调度当作一个整体MILP来解。决策变量数量大是表层原因,深层原因是模型里存在大量强耦合的二元变量(电解槽启停、合成氨启停)和连续变量(功率、SOC、氢量),导致分支定界树非常深。
解决路径:这个坑我花了整整三天才彻底走出来。最终用的是“时间分解迭代法”——先按季度把全年分成12个典型场景,每个场景挑一个代表性周期分别求解,再把结果合并计算整年的成本和产量。虽然严格来说这不是数学意义上的全局最优,但工程上足够接近,而且求解时间从4小时压缩到了20分钟以内。
经验:复现论文代码时,先跑一个100小时的小算例验证模型正确性,再用完整8760小时跑最终结果。一上来就跑全规模,一旦模型有bug,你连等待结果的时间都耗不起。
5.2 问题二:离网模式下模型显示不可行
现象:模式切换到离网后,Cplex直接返回“infeasible problem”,连一个可行解都找不到。
根因分析:我缩小了风光装机上限、蓄电池容量上限,再叠加了全年所有时刻功率严格平衡的约束,导致无解。这是离网系统天然的挑战——如果没有电网兜底,任何一段持续的无风无光期都会让问题不可行。
解决路径:排查链路是这样的——先用optimize检查infesibility,再用check命令逐条检查约束是否全部满足(虽然infeasible时没有变量值,但Yalmip会帮你定位是哪一组约束导致了冲突)。我发现是蓄电池SOC的首尾相接约束和全年电量下限约束同时存在,导致了一个“要么首尾相连,要么全年SOC不落地”的死结。后来我把SOC的全年电量下限去掉,只保留首尾相接,问题就解决了。
经验:离网系统的蓄电池SOC约束是模型里最敏感的约束。论文里写的“SOC在0.1到0.9之间”看起来只是上下限,但和首尾相接组合后,相当于强制全年净充电量为零。这在物理上没问题——因为蓄电池作为一个储能装置,不应当从外部获得净能量——但它的确极大收缩了可行域。遇到不可行时,优先检查SOC相关的三条约束(递推、上下限、首尾相接)是否冲突。
5.3 问题三:线性化处理不当导致结果失真
现象:电解槽效率、氢气压缩能耗等环节,我在初版模型里用固定效率处理,结果算出来的最优容量配置和论文结果相差超过20%。
根因分析:固定效率的线性模型忽略了电解槽在低负荷工况下效率会明显下降这一事实。真实碱性电解槽在20%负荷时效率可能只有满负荷时的85%左右。如果把这个非线性关系强行线性化为固定效率,就等于无偿给了系统一笔“效率补贴”,导致优化结果偏向于小容量高利用率的配置。
解决路径:论文里的处理方式是分段线性化。把电解槽运行范围分成3~5个区间,每个区间用一个线性效率近似。在Yalmip里可以用implies或者直接引入分段线性约束来实现。校验的方法很简单:用优化出的容量配置回代到真实非线性效率模型里,对比产氢量偏差,偏差超过2%就必须增加分段数。
经验:任何“看起来重要”的非线性关系,在MILP里都是靠分段线性化处理的。但分段的数目不是越多越好——段数越多,二元变量越多,求解时间越长。对电解槽效率这种曲线比较平滑的关系,5段足够;对光伏效率随温度变化这种近似线性的关系,2段就够。复现时先跑一次2段和5段的对比,如果最优解变化很小,就用少段数。
5.4 问题四:数值奇异导致求解结果跳动
现象:同样的模型跑两次,容量优化结果完全一样,但调度结果略有不同;更奇怪的是,某个约束变量的数值量级跨度从0.0001到1e6,Cplex求解过程中频繁发出数值警告。
根因分析:这是典型的“数值缩放没做好”问题。我的模型里,风电容量单位是kW(量级1e3),而储氢罐容量单位可能是kg(量级1e1),氨产量单位是吨(量级1e-1),三个变量的量级差了好几个数量级。Cplex内部用的是浮点数,矩阵里出现1e-6和1e6同框时,数值误差会被放大,导致求解器容忍度失效。
解决路径:把所有变量和参数统一到同一套量纲体系。我最终的选择是:功率统一用MW,氢气量用吨,氨产量用吨,时间步长用小时,价格用元/MWh和元/吨。这样所有变量的数值都在0.001到1000之间,数值稳定性好很多。改完之后,Cplex的数值警告基本消失,求解结果也稳定了。
经验:建模第一件事不是写公式,而是先把每个变量的单位定死。论文里经常混用kW/MW、kg/t、元/kWh/元/MWh,复现代码时必须全部换算好。尤其在用大M法时,M的值必须和周围变量的量级匹配,一个错误的量级会让你得到一堆数学上正确、物理上荒谬的解。
5.5 问题五:并网电价策略对结果影响被忽略
现象:我起初全篇用一个固定平均电价,结果容量优化明显偏向“少配储能、多依赖电网”。直到做灵敏度分析、按分时电价重新跑通之后,结果才合理。
根因分析:这是我对电价机制理解不够深入导致的。分时电价背后是“峰谷价差套利”逻辑——蓄电池和储氢罐在电价低谷时充电/产氢,在电价高峰时放电/售氢,这个套利行为直接影响最优的储能容量配置。固定电价把这个套利空间完全抹掉了,系统的储能规模自然会萎缩。
解决路径:完全复现论文的电价曲线,通常由峰平台谷四个时段构成。我在代码里用repmat把一天24小时的分时电价扩展到全年8760小时,还加入了节假日和周末的电价差异。重新求解后,蓄电池最优容量比固定电价情景下高出30%~50%,这更加接近论文报告的数值。
经验:电力市场的价格信号是容量优化的一大驱动力。做这类系统优化,电价曲线不仅仅是一个参数,它是和储能容量配置直接耦合的“隐藏变量”。如果论文里提到了分时电价,那这个电价必须作为时序数据传入模型,而不是用一个标量替代。
6. 结果验证与灵敏度分析:怎么判断你的复现是否正确
跑出优化结果只是复现的一半,另一半是验证——如何证明你的模型结果不是“错误代码运行出的漂亮数字”。
6.1 能量守恒检验:最基础的合理性校验
模型跑完后,我做的第一件事是全年能量平衡校验。把全年风电光伏总发电量、电解槽总耗电量、蓄电池总充电量、总放电量、电网总购电量、总售电量、辅助负荷总耗电量全部算出来,检查下面这个等式是否成立:
$$E_{wind} + E_{pv} + E_{buy} + E_{dis} = E_{el} + E_{load} + E_{sell} + E_{ch}$$
如果偏差超过1%,那一定是模型里某条约束写错了。我做这个校验时抓到过两个bug:一个是蓄电池充放电效率写反了(充电时用了放电效率),另一个是储氢罐的氢气“泄漏”——物料平衡约束的间隔口径和储能SOC递推差了1个小时,导致全年氢量总量不对。
6.2 边界场景法:逐步逼近极端工况
边界场景法的思路很简单——不断把参数推向极端,看模型的反应是否符合物理直觉。
- 把光伏资源全部置零,系统应该自动把光伏装机优化为0,而风电装机被迫增大以补偿发电缺口
- 把电解槽投资成本拉到极高,系统应该削减电解槽容量,甚至完全不做制氢转而直接买氨(如果允许的话)
- 把氨价格降到极低,系统应该放弃合成氨,只做纯制氢
- 把电价拉得极高,并网系统应该减少购电、增加自发电和储能规模
如果某个极端场景下模型给了一个违背直觉的结果,那大概率是约束条件里缺了一环或者参数定义错了。边界场景法可以在几分钟内暴露模型逻辑问题,比直接看完整结果高效得多。
6.3 与论文结果对比:差异在合理范围内才算复现成功
论文给出的最优容量配置和年化总成本是我校验的终极标准。复现结果和论文结果差异在5%以内,算是一次成功的复现;5%~15%的差异,可能是设备参数、电价数据、资源数据的取值不同导致的,需要仔细核对;超过15%的差异,基本可以确定你的模型结构或数据预处理出了系统性问题。
我在核对参数时发现,论文里“电解槽系统效率”一词有歧义——它可能指电解槽单体的电流效率,也可能指包含整流器、变压器、水泵在内的系统总效率。两者差值通常在5~10个百分点,对最优容量配置影响很大。复现者必须结合论文的上下文和引用文献推断出真正的定义,必要时对比多篇同方向论文取合理值。
6.4 灵敏度分析:验证系统应对不确定性的能力
复现完成后,我还额外做了一组灵敏度分析,这是论文正文没有的,但对理解系统很有帮助:
- 风光资源年际波动(好年份、正常年份、差年份)对最优容量影响
- 贴现率从4%调到8%时,对长寿命设备(光伏、风电)和短寿命设备(电解槽)配置比例的影响
- 氨价格变化对系统是否“值得建设”的经济阈值影响
灵敏度分析的价值在于:它告诉你最优解在一组参数下的“稳健性”。如果贴现率从4%变到6%,最优容量配置就剧烈变化,那说明模型对贴现率极度敏感,在真实世界做投资决策时必须非常谨慎地标定这个参数。这类信息在论文里往往不直接写,但对实际工程项目极有价值。
7. 扩展方向:这台模型的下一步还能做什么
复现完论文模型之后,我并没有停下来。在实际操作中我体会最深的一点是:一套风光互补制氢合成氨优化模型的价值,不在于它能算出一个漂亮的数字,而在于它能作为一个平台承载更复杂的真实问题。下面几个方向是我认为最值得尝试的扩展。
第一,把确定性优化升级为随机优化或鲁棒优化。风光的随机波动性是这个系统的天生属性,单条确定性出力曲线算出的最优容量配置,在真实场景中可能不足。可以考虑典型场景法(多个出力场景加权求和)或两阶段随机规划,前者简单直接,后者更严谨但求解规模大增。
第二,加入电解槽和合成氨装置的详细退化模型。当前模型假设设备在整个寿命期内性能不变,但实际上电解槽催化剂会衰退、质子交换膜会老化、氨合成催化剂会失活。把这些退化因素引入模型,可以让“到第几年该更换催化剂”成为模型的内生决策。
第三,把电网侧的约束引入模型,不只是无限容量的购售电。实际上,有些地区的电网接入容量有限制,有些地方有自消纳比例要求(比如绿电制氢必须保证一定比例的本地消纳率),这些约束会显著影响最优容量配置,值得作为条件约束加入模型。
第四,把碳排放约束纳入目标函数。随着碳交易市场在更多地区落地,碳排放成本会成为绿色氢氨项目经济性的一个重要变量。把碳价作为惩罚项加入目标函数,可以完整评估“绿色溢价”对项目可行性的影响。这种扩展在双碳背景下有很强的现实意义。
我在跑完整套复现和扩展之后的最大感受是:这套模型的方法论价值大于数值结果本身。容量-调度一体化优化本质上是一个“投资决策与运行策略协同优化”的问题,它的方法论可以套用到任何涉及长期资产配置和短期运行调度的场景——不只是制氢合成氨,还包括电化学储能电站的设计、综合能源系统的规划、甚至数据中心的新能源供电方案。能把这个模型的骨架搭明白,相当于掌握了一套可以迁移到多个新能源场景的核心分析能力。