☰
并网/离网风光互补制氢合成氨容量调度优化:Matlab+Cplex实战
2026/10/7 3:42:20 网站建设 项目流程

风光互补制氢再合成氨,这个链路这两年在我接触的复现咨询里出现频率很高。直接储氢又贵又麻烦,把氢转化成氨再卖,产业链路一下子就顺了。最近我完整复现了一套并/离网风光互补制氢合成氨系统的容量-调度优化分析,Matlab建模、Cplex求解。这套题目看起来就是一个标准的可再生能源系统优化,真正跑起来才发现,风光储的坑它一个不落,化工装置连续生产的约束又额外加了一堆麻烦。这篇文章把从建模到求解的完整思路、Cplex在Matlab里的落地写法、并网和离网两种模式的差异对比,以及我踩过的坑全部展开讲清楚。打算做风光氢储优化、综合能源系统仿真的研究生,或者需要在Matlab里调Cplex解混合整数规划的朋友,这篇应该能帮你省下不少试错时间。

1. 从风到氨:先搞懂这套系统在优化什么

1.1 风光出力、电解制氢、储氢与合成氨的耦合关系

这套系统的物理链路并不复杂,但每一步之间的耦合关系决定了优化模型的形态。

风电机组和光伏阵列发出来的电,优先供给电解槽制氢。电解槽产出的氢气有两条去向:一部分直接进氢气储罐缓冲,另一部分和空分装置提供的氮气一起进入合成氨装置,发生Haber-Bosch反应生成氨,氨再进液氨储罐存储,最终外售。电解槽还会副产氧气,多数模型里直接弃掉或作为副产品不计收益。

这里有个容易忽视的点:氢储罐和氨储罐承担的时间尺度完全不同。氢储罐平衡的是小时级别的风光波动和电解槽负载波动,比如夜里光伏停了、风速又低,就需要靠储罐里的氢维持合成氨装置的进料;氨储罐平衡的是日级别甚至周级别的产量与销售节奏。所以建模时,储氢和储氨的动态方程都必须写,不能只写一个总储罐糊弄过去。

另外一个关键点是,为什么终点选氨而不是直接卖氢?因为氢的储运成本太高,氨是成熟的化工产品,液氨的储运基础设施完善,终端可以作为燃料、化工原料或者再裂解成氢。所以这套系统本质上是一个"电力-氢-化工产品"的多品级能源枢纽,目标函数里必须体现产品销售收益,否则优化结果会偏向过度储电而忽视氢氨产量。

1.2 并网与离网:两种边界条件下运行逻辑的本质差异

并网模式和离网模式看起来只是"有没有电网交互功率"这一个约束的区别,实际上整个系统的设计哲学都变了。

并网模式下,电网是一个无限容量的"缓冲水池"。风光不够时,可以从电网买电维持电解槽和合成氨装置运行;风光多了,多余的电可以卖给电网。这时容量配置可以相对紧凑,因为缺电风险被电网化解了。调度策略的核心变成了"电价信号驱动"——电价低的时段尽量多用电制氢储起来,电价高的时段减少购电甚至反向售电。

离网模式下,系统是一个孤岛,任意时刻的功率都必须自平衡。风不吹、光不照的时段,只能靠储氢罐和氨储罐的存量硬撑。这意味着容量配置必须冗余,风机和光伏的装机要远超平均负荷需求,否则极端低出力时段合成氨装置就得停机。调度策略的核心从"经济性最优"变成了"供给可靠性和弃电最小化"——发出来的电尽量别弃,电解槽尽量满负荷。

这两种模式对应的是两个极端的产品属性:并网模式追求单位氨成本最低,离网模式追求生产过程的完全独立性和低碳属性。复现的时候,我建议把并网跑通之后再去掉电网交互变量跑离网,这样对比起来非常直观。

2. 容量-调度优化问题的数学化过程

2.1 决策变量、目标函数与资本回收因子的处理

这个题目叫"容量-调度优化分析",我复现时第一件事是判断它到底是严格的双层优化还是单层协同优化。如果论文用的是Cplex直接求解,绝大多数情况下是单层混合整数线性规划(MILP),也就是把容量决策变量(慢变量)和调度决策变量(快变量)放在同一个模型里一起优化。

容量层决策变量包括:

  • 风电机组装机台数(整数变量)
  • 光伏阵列容量(连续变量)
  • 电解槽容量(连续变量)
  • 氢储罐容量(连续变量)
  • 合成氨装置产能(连续变量)
  • 液氨储罐容量(连续变量)

调度层决策变量是每个时段的:

  • 风电机组出力、光伏出力
  • 电解槽输入功率与产氢量
  • 氢储罐的充放氢速率
  • 合成氨装置的氨产量与耗氢量
  • 并网模式下与电网的交互功率(购电为正、售电为负,或用两个变量分开)

目标函数是年化总成本最小,典型形式为:

min 年化设备投资成本 + 年运行维护成本 + 年购电成本 - 年售电收益 - 年氨销售收入

投资成本不能直接拿初始投资额相加,因为风电设备寿命20年、电解槽可能只有10年,资金有时间价值。这里必须引入资本回收因子(CRF,Capital Recovery Factor):

CRF = r(1+r)^n / ((1+r)^n - 1)

其中r是折现率,n是设备寿命。我复现时折现率取8%,风电机组寿命20年,CRF约0.1019;电解槽寿命10年,CRF约0.1490。这意味着同样是一万元投资,电解槽每年的成本负担比风机高出近50%。很多复现结果对不上,往往就是这里把年化和一次性投资搞混了。

另外注意数量级问题。投资成本通常是几百万到几千万,调度成本可能是几十块到几百块。目标函数里如果直接用原始数值,两者相差10的6次方以上,Cplex的数值稳定性会变差,求解速度明显下降。我习惯把所有成本统一到"万元"为单位再建模。

2.2 约束条件里最容易出错的几个环节

模型的正确性基本都在约束里,我按出错频率排个序。

电力平衡约束是所有约束的基础。每个时段都必须满足:

P_wt(t) + P_pv(t) + P_grid(t) = P_el(t) + P_curtail(t)

其中P_grid在离网模式下直接被赋值为0。这里有个细节:如果不加弃电变量,低出力时段模型会强行让电解槽降低功率甚至让光伏限电,但实际中弃电是物理存在的,所以P_curtail是一个非负变量,而不是等式约束里忽略的松弛项。

电解槽运行约束是第二个坑点。电解槽不是任何功率都能运行的,碱性电解槽的负载范围一般在20%到100%之间。如果你只写0 <= P_el <= Cap_el,模型会利用0到20%这段不合理的区间去"蹭"平衡,导致结果偏乐观。严格起见,需要引入开机状态二元变量u_el(t),写成:

P_el(t) >= 0.2 * Cap_el * u_el(t)P_el(t) <= Cap_el * u_el(t)

氢储罐动态约束是第三个坑。储罐状态方程是:

SoC_h(t+1) = SoC_h(t) + H_ch(t) * eta_ch - H_dis(t) / eta_dis - H_ha(t)

其中H_ha是合成氨装置每个时段消耗的氢气量。注意充氢和放氢不能同时进行,这需要一对二元变量。但我的经验是:如果储罐的充放效率不是模型的核心关注点,可以用"储罐净变化量 + 充放损耗系数"的方式简化,去掉二元变量,模型规模直接下降一大截。

合成氨装置的约束会让人非常头疼。Haber-Bosch反应是连续化工过程,催化剂床层温度和压力都不能频繁大幅波动,所以装置通常有最小负载率约束,比如0.3 * Cap_ha <= NH3_prod(t) <= Cap_ha。如果需要更严格,还要加最小运行时间和最小停机时间的约束,但这些都是整数变量,变量规模会爆炸。我复现时的折中方案是:先不加启停时间约束,只加最小负载率,跑通后再逐步加严。

还有一个化学计量关系必须写对。合成氨反应是N2 + 3H2 -> 2NH3,按质量算,1kg氢气最多生成约5.667kg氨。如果论文里给的转换系数和这个偏差超过10%,大概率是单位或者化学配比搞错了。

2.3 典型日筛选:全年8760小时怎么压缩

这个问题很多复现新手会忽略,直接拿全年8760小时跑,结果模型变量几十万,Cplex跑几个小时都出不来。

我的做法是先用K-means聚类筛选典型日。把风速、光照、温度(如果影响负荷)、并网模式下的电价作为聚类特征,把全年8760小时聚成若干个代表日,比如春夏秋冬各取一个典型日,或者聚类成12个典型日。每个典型日乘以对应的天数权重,就能近似代表全年。

这里有一个非常隐蔽的坑:如果只选典型日,合成氨装置跨日连续运行的特性就会被破坏。比如典型日1的最后一天晚上氢储罐是满的,但下一个典型日的初始状态怎么接?如果不做处理,模型会让储罐状态在每个典型日之间任意跳变,相当于免费获得了一个"跨日调节"能力,优化结果偏乐观。

解决方案有两种:一是用典型周(7天一组)代替典型日,代价是聚类维度和求解规模增大;二是保留典型日,但把储氢罐和氨储罐的初末状态约束成相等,或者加跨典型日的状态连接约束。我实际复现时用的是第二种,因为变量规模可控,结果和论文对得上。

3. Matlab调用Cplex求解:落地路径和代码骨架

3.1 环境配置:Cplex安装与Yalmip接入

Cplex是IBM的商业求解器,从官方渠道下载IBM ILOG CPLEX Optimization Studio并安装后,Matlab调用它有两条路。

第一条是直接调用Cplex自带的Matlab接口,比如cplexmilp函数,把所有约束转成矩阵形式A*x <= b和Aeq*x == beq传入。这条路的问题是:模型约束一多,手写矩阵根本维护不了,改一个约束索引全乱,调试体验极差。

第二条路是装Yalmip工具箱,用符号变量建模,把约束集合写成对象表达式,最后交给Cplex求解。Yalmip会自动转换成Cplex需要的标准形式。这条路是所有做MILP优化的人第一推荐的方式,没有之一。

安装步骤很简单:

  • 安装CPLEX后,把...\CPLEX_Studio221\cplex\matlab目录加入Matlab路径
  • 下载Yalmip源码,把整个yalmip目录加入Matlab路径
  • 在Matlab里运行yalmiptest,看到cplex那一栏显示成功就说明通了

版本匹配是个老坑。CPLEX版本要和Matlab版本大致兼容,特别是旧版CPLEX在新版Matlab上经常报"未能加载库"的错。出现这类问题不用慌,优先检查CPLEX版本是否支持当前Matlab,或者是不是把cplex的Java库路径漏配了。

3.2 Yalmip建模与Cplex求解的代码骨架

我给出一个实际用的代码骨架,核心结构可以直接套用:

%% 基本参数 T = 24 * 365; % 总时段,按小时计 dt = 1; % 步长1小时 CRF_wt = 0.1019; % 风机资本回收因子,8%折现率20年 CRF_el = 0.1490; % 电解槽资本回收因子,8%折现率10年 inv_wt = 450 * 1e4; % 风机单位投资,万元/MW % ... 其他设备参数 %% 典型日索引生成 % 通过K-means聚类得到典型日编号idx和权重weight % 实际运行时段数为 idx 展开后的 T_rep %% 决策变量 N_wt = intvar(1, 1); % 风机台数,整数变量 Cap_pv = sdpvar(1, 1); % 光伏容量 MW Cap_el = sdpvar(1, 1); % 电解槽容量 MW Cap_hs = sdpvar(1, 1); % 氢储罐容量 kg Cap_ha = sdpvar(1, 1); % 合成氨装置产能 kg/h P_grid = sdpvar(T, 1); % 电网交互功率 MW,正值购电负值售电 P_el = sdpvar(T, 1); % 电解槽输入功率 MW H_prod = sdpvar(T, 1); % 产氢速率 kg/h H_dis = sdpvar(T, 1); % 储氢罐放氢速率 kg/h H_ch = sdpvar(T, 1); % 储氢罐充氢速率 kg/h SoC_h = sdpvar(T+1, 1); % 氢储罐状态 kg NH3_prod = sdpvar(T, 1); % 氨产量 kg/h % P_curtail等变量略 %% 约束集合 Constraints = []; % 电力平衡 Constraints = [Constraints, P_wt + P_pv + P_grid == P_el + P_curtail]; % 电解槽功率上下限 Constraints = [Constraints, 0.2 * Cap_el <= P_el <= Cap_el]; % 产氢量等于电解槽功率除以单位电耗 Constraints = [Constraints, H_prod == eta_el * P_el / 39.4]; % 储氢罐动态 Constraints = [Constraints, SoC_h(2:T+1) == SoC_h(1:T) + H_ch - H_dis - H_ha]; % 储氢罐容量约束 Constraints = [Constraints, 0 <= SoC_h <= Cap_hs]; % 化学计量关系:1kg H2 -> 5.667kg NH3 Constraints = [Constraints, NH3_prod == H_ha * 5.667]; % 合成氨装置负载范围 Constraints = [Constraints, 0.3 * Cap_ha <= NH3_prod <= Cap_ha]; % 并网交互功率上限 Constraints = [Constraints, -P_sell_max <= P_grid <= P_buy_max]; %% 目标函数 % 年化投资成本 + 运维成本 + 购电成本 - 售电收益 - 氨销售收入 Objective = CRF_wt * inv_wt * N_wt ... + CRF_el * inv_el * Cap_el ... + ... + dt * sum(price_buy .* max(P_grid, 0)) ... - dt * sum(price_sell .* max(-P_grid, 0)) ... - dt * sum(price_nh3 .* NH3_prod); %% 求解设置 ops = sdpsettings('solver', 'cplex', 'verbose', 2); ops.cplex.mip.tolerances.mipgap = 0.01; sol = optimize(Constraints, Objective, ops); %% 结果读取 if sol.problem == 0 N_wt_opt = value(N_wt); Cap_pv_opt = value(Cap_pv); % ... 读取其他变量 else disp(sol.info); end

注意一个问题:目标函数里的max(P_grid, 0)不是线性表达式,Cplex处理起来会有麻烦。我实际建模时会把购电和售电拆成两个变量:

P_buy(t) >= 0, P_sell(t) >= 0, P_buy(t) - P_sell(t) == P_grid(t)

然后在目标函数里对P_buy和P_sell分别计价。这样模型保持线性,Cplex求解最稳定。

3.3 Cplex求解器关键参数与调优方向

Cplex默认配置追求的是"最优性证明",但对于工程优化模型,我们往往不需要1e-4那么极端的MIP gap。实测下来,把gap放宽到1%就能让求解时间缩短一个数量级,而优化结果几乎没差别。

我常用的几个关键参数:

参数推荐值作用
mip.tolerances.mipgap0.01允许1%的次优解,大幅提速
mip.tolerances.integrality1e-5控制整数变量容差,默认值即可
threads4或8并行求解线程,多核CPU明显加速
time3600设置求解时间上限,防止卡死
mip.strategy.startalgorithm4使用barrier算法起步,对大规模LP有帮助

在Yalmip中设置方式是ops.cplex.xxx = yyy,具体key和CPLEX文档一致。另外,如果求解器报告内存不足或者想中断,随时可以通过sol.info查看求解状态,Yalmip返回的sol.problem == 0表示成功找到最优解,其它值都有对应的错误码。

4. 并网与离网调度策略的对比实验

4.1 两种模式在调度层约束上的差异

把同一套代码拿去跑并网和离网,本质上只需要改三个地方。

第一,电网交互变量。并网模式下保留P_buy和P_sell,并加上交互功率上限;离网模式下直接赋零。我习惯在代码里用一个模式开关变量:

if strcmp(mode, 'grid') Constraints = [Constraints, P_buy - P_sell == P_grid_net]; Constraints = [Constraints, 0 <= P_buy <= P_buy_max]; Constraints = [Constraints, 0 <= P_sell <= P_sell_max]; else Constraints = [Constraints, P_buy == 0, P_sell == 0]; end

第二,目标函数中的购售电项。并网模式要加上购电成本和售电收益,离网模式直接把这两项删掉,只保留投资、运维和氨销售收入。

第三,关键的是容量约束的"隐性变化"。离网模式下,如果任何一个时段的风光出力为0且储氢罐放空了,合成氨装置就会断料。所以离网模式对储罐容量的需求远高于并网模式。如果不提高储罐容量上限,模型大概率infeasible。

这两类约束差异直接反映在最优调度行为上。并网模式的电解槽往往随电价波动:电价低谷时拉满功率制氢储氢,电价高峰时降低功率甚至停机。离网模式的电解槽则是跟随风光出力:风大光强时满负荷甚至让部分风光弃电,风小光弱时靠储氢罐维持合成氨装置的最低负荷。

4.2 典型结果解读:容量配置、成本拆解与储罐动作曲线

我跑出来的结果趋势基本符合预期,这里给一组典型对比:

指标并网模式离网模式
风电装机(相对并网)基准约1.6倍
光伏容量(相对并网)基准约1.4倍
电解槽容量基准约1.3倍
氢储罐容量基准约2倍以上
弃电率接近015%以上
年化总成本较低高30%左右
单位氨生产成本低高

这个结果逻辑上完全说得通。并网模式用电网兜底,容量冗余小,购电成本换投资成本;离网模式必须靠自己硬扛低出力时段,所以风电光伏装得多、储罐配得大,弃电率也下不来。

调度曲线方面,并网模式最典型的特征是电解槽功率和电价曲线呈镜像关系。谷电时段购电制氢,峰电时段卖电或者用储氢维持生产。离网模式的SoC曲线则明显跟着风光出力走,连续几天低温差天气时,储氢罐的SoC会一路下滑到下限附近,这时候如果风速预测没偏差,模型就会在目标函数里狠狠地惩罚容量不足。

这些结果也提醒了一个事:单纯比较"并网好还是离网好"没有意义,关键是看应用场景要的是什么。如果追求最低氨成本,并网模式赢;如果强调生产过程完全独立、不受外部电网约束,离网模式是唯一选择。很多论文做的是离网模型却拿并网电价做经济性论证,这属于逻辑不自洽,复现时要注意。

5. 复现复盘:最容易让结果翻车的五个问题

5.1 单位不一致导致的约束错乱

这个坑我几乎每次复现都会碰到一次。功率用MW,电量用MWh,氢储能单位又变成kg,氨产量单位是t/d,几个单位混在一起,约束的系数稍微错一个数量级,结果就跑偏了。

我的习惯是开工之前先把单位表列好:功率统一MW,时间统一h,能量统一MWh,氢气统一kg,氨统一kg。如果论文里氨产量用的是t/d,转成kg/h的系数就是除以24再乘以1000。电解槽产氢量的换算系数也要提前确认,用的是氢气的高热值还是低热值,直接影响电解效率的取值。这个系数差个百分之十几,优化结果容量配置就会明显不同。

5.2 储能二元变量与大M取值

储氢罐充放互斥约束,理论推导是标准的:

H_ch(t) <= M * u_ch(t)H_dis(t) <= M * u_dis(t)u_ch(t) + u_dis(t) <= 1

这里的M取值是个讲究活。M取太小会砍掉可行解,M取太大(比如1e9)会引发数值病态。工程上按储罐最大充放速率的两倍取就够了。比如储罐最大充氢速率是500 kg/h,M取1000,完全够用。

我个人的偏好是:除非论文明确要求储罐充放效率差异显著,否则我直接不看充放互斥约束,用净充放变量代替。因为在一个小时级别的时间尺度上,储罐同时充放本来就是物理上不可能的,但经济最优解几乎不会出现这种病态行为,删掉这组二元变量可以显著降低MILP规模。

5.3 求解性能瓶颈与MIP gap取舍

容量-调度优化模型最容易卡死的地方在于时间步长和二元变量数量的乘积。如果T=8760并且每个时段都有电解槽启停、储罐充放、合成氨装置开停三组二元变量,二元变量数量超过2.6万个,Cplex也需要相当长的求解时间。

我的处理顺序是先不加化工装置的启停约束,只保留最小负载率约束,阶段一快速跑通模型拿到可行解和成本量级;再逐步加严约束,观察求解时间变化。如果求解时间超过30分钟且gap还在5%以上,我一般会把mipgap调到0.03甚至0.05,先拿一个"工程可接受"的次优解,再针对性优化瓶颈约束。

5.4 离网模式下等式约束过紧导致的不可行

离网模式最常见的报错是solver returned infeasible。这时不要急着调参数,先检查是不是某些时段功率平衡等式和储罐上下限联合起来把可行域压缩成了空集。

具体场景是:连续几天低风速,光伏夜间又归零,电解槽没法产氢,储氢罐已经放空,合成氨装置的最小负载率必须维持,这时候约束组直接矛盾。处理办法有三个:一是调大储氢罐容量上限;二是降低合成氨装置最小负载率(但可行性存疑);三是允许在极端场景下合成氨装置停机,并引入停机惩罚成本。第三个方案其实最贴近物理实际,但需要加二元变量。

5.5 与论文结果对不上的排查顺序

复现论文最难受的就是"模型能跑但结果和论文差很多"。我的排查顺序固定如下:

第一步,核对目标函数的成本和收益项,看有没有漏掉或重复计入的项。这是最常见的偏差来源。第二步,核对化学计量系数和经济转换系数,特别是"1kg氢产多少氨"这种基础常数。第三步,核对约束方向,尤其是不等式符号写反或者上下界颠倒了,会导致完全不同的最优解。第四步,核对时间序列数据,比如典型日聚类是不是把某一段极端天气的代表日漏掉了。

如果以上都查完还对不上,那大概率是论文本身对某些约束做了简化,或者参数取的具体数值在附录里没写清楚。这时候我会反推:根据论文给出的容量配置结果,反解它的隐含成本参数,把数值调到自己模型里,看能不能逼近它的结果。这个办法虽然有点笨,但非常见效。

从实际工程角度看,这类系统的优化求解最大的价值不在于那个"最优容量数字"本身,而在于模型把不同时间尺度、不同物理过程的耦合约束完整体现出来之后,你能清晰地看到成本和可靠性的矛盾在哪里。Cplex在这里扮演的角色就是一个可靠的MILP求解内核,而Matlab和Yalmip负责把模型快速翻译成求解器能理解的语言。最后再提醒一句:跑这类模型,别一上来就追求全年8760小时无所不包,先跑典型日、先跑并网、先不加启停约束,全链路跑通之后再逐层加复杂度,这是能让你少掉不少头发的最实际建议。

需要专业的网站建设服务?

联系我们获取免费的网站建设咨询和方案报价,让我们帮助您实现业务目标

立即咨询