☰
风光制氢合成氨系统优化:模型原理与Pyomo工程化实现
2026/10/26 7:31:09 网站建设 项目流程

拿到“风光制氢合成氨系统优化”这个题目,很多人的第一反应是先去找Python代码,然后直接把论文里的公式往里面套。但真正花过几周时间在“论文复现”上的人会告诉你,这个项目最折磨人的点从来不是求解器跑得慢,而是模型本身横跨了电力、氢能和化工三个领域,每个领域的文献都有自己的语言体系和习惯假设,第一步根本不应该是写代码,而是把系统运行方式和优化目标彻底搞清楚。

这篇文章我会从复现论文的视角,把整个风光制氢合成氨系统的数学模型、Python工程化实现、常见的坑和结果呈现方式完整拆一遍。内容面向正在做这类课题的研究生、从事新能源制氢方案设计或投资测算的工程师,也适合想学习用Pyomo/PuLP求解多时段优化问题但一直没找到完整案例的Python用户。

1. 从一张系统架构图说起:这个优化问题解决的是什么事

1.1 风光如何变成氨:能量转换链与各单元作用

风光制氢合成氨系统,字面上就包含三件事:风光发电、电解水制氢、合成氨。完整能量链条大概是这样的:

风电和光伏发出交流/直流电,经过汇流母线后送入电解槽;电解槽把水拆成氢气和氧气,氢气进入储氢罐缓冲;储氢罐再按照合成氨单元的消耗节奏供氢;空分装置提供氮气,氢氮混合气在合成塔内高温高压下反应生成氨。产品氨可以作为化工原料或燃料出售,这也是整个系统的主要收入来源。

这套流程里隐含了两个必须同时满足的平衡:电力平衡和物料平衡。电力平衡要求每一时刻风电、光伏、电网购电的功率总和,等于电解槽消耗加上弃电功率;物料平衡则要求氢气的“产、储、用”在每一时刻都保持守恒。很多复现代码跑出无解或者结果明显不合理,根源就是只写了电力约束而漏了氢平衡,或者反过来。

各单元的角色也不一样。风光机组是“看天吃饭”的电源,出力曲线不可控,只能被动接受;电解槽是灵活负荷,可以通过调节输入功率改变产氢量,但有最低负载率限制,不能无限降载;储氢罐是缓冲器,它让制氢和用氢不必实时匹配,这是整个系统具备调节能力的关键;合成氨单元是典型的化工装置,催化剂和反应温度决定了它不喜欢频繁启停,在优化模型里通常要求保持连续稳定运行或至少不低于某个最小负荷率。

1.2 为什么说这是一个带时序耦合的优化问题

风光制氢合成氨和普通电力调度最大的区别,在于储氢罐把不同时刻的决策“黏”在了一起。今天多产出来的氢可以存到明天用,明天产少了可以消耗昨天存的氢,这意味着第t时刻的最优决策依赖于第t-1时刻的储氢状态,而t-1时刻的状态又依赖于更早的决策。这种跨时段耦合是典型的“带库存的生产计划”问题,比单纯求解每个时刻的静态优化要复杂一个量级。

你可以把它类比成一家靠天吃饭的食品加工厂:原材料供应取决于天气(风光出力),生产线有最低开工率要求(电解槽最低负载),仓库容量有限(储氢罐),下游门店要求你持续供货不能断(合成氨连续运行)。每过一段时间就要做一次排产计划,既要保证产线不空转,又要控制库存成本,同时还要保证下游不停工。把这种问题写成数学模型,就变成了带时序递推约束的优化问题。

在数学形式上,这种问题通常写作混合整数线性规划(MILP)。为什么会有整数变量?因为电解槽和合成氨装置都有“开/关”状态,最低负载约束必须借助01变量才能表达,比如电解槽要么停机,要么至少以20%额定功率运行,不允许在0到20%之间连续调节。这类约束一旦出现,问题就从LP变成了MILP,求解难度和求解时间都会显著上升。

1.3 复现范围界定:先跑通典型日再上全年

风光制氢合成氨的论文里,优化周期通常有两种选择:典型日(通常24小时)和全年(8760小时)。典型日的好处是计算量小、结果容易画图呈现、方便做敏感性分析,论文里最常见的就是选取夏季晴天、冬季阴天等若干典型场景分别优化。全年8760小时更多用于经济性评价,会涉及年发电量、设备利用小时数、年化成本和收益的计算。

我复现这类项目时有一个很重要的习惯:第一版代码永远只跑24小时,先把单位、量纲和模型逻辑全部理顺,确认目标函数值和论文对得上,再扩展成168小时或者全年。直接一上来就写全年8760小时的模型,但凡有一个单位换算错误,跑出来的结果会偏差得非常离谱,而且根本没法排查。

典型日数据从哪来?论文附录里经常会给关键时段的功率曲线,有些会给出典型日的合成氨负荷和电价。如果没有,可以用实际风电场光伏电站的实测数据做缩放,比如把某个风电场某天的出力按装机容量折算成你要用的风电规模。数据精度能反映日内的波动特性就够了,不需要追求和原文完全一致,因为复现的重点是验证模型逻辑,而不是复现数据的每一处细节。

2. 模型数学骨架:目标函数、决策变量与约束的取舍

2.1 目标函数:成本最小化模型里藏着哪些明细项

风光制氢合成氨优化论文里最常见的建模目标是最小化系统总成本,把卖氨收入写成负成本,就变成了一个“成本最小化”的形式。等你把各项拆开,会发现里面有四种性质完全不同的成本,单位还不统一,新手第一道门槛就在这里。

目标函数可以写成下面这种形式:

min C = C_inv + C_om + C_grid + C_punish - R_NH3

各项的含义和计算方式如下:

  • C_inv:设备投资年化成本。风电、光伏、电解槽、储氢罐、合成氨装置都有各自的单位容量投资(元/kW或元/kg),乘上容量得到总投资,再用资金回收系数(CRF)折算成年值。
  • C_om:运行维护成本,通常按设备投资的一定比例估算,常见取值2%到5%。
  • C_grid:购电成本,等于每个时段的电价乘电网购电功率再乘时段时长。
  • C_punish:弃风弃光惩罚项。有些模型把弃电设成硬约束不允许发生,有些模型允许弃电但加惩罚,让优化器尽量避免它。
  • R_NH3:卖氨收益,由氨产量乘氨价得到。由于它在目标函数里是减号,形式上相当于“负成本”。

资金回收系数CRF的算式是CRF = i(1+i)^n / [(1+i)^n - 1],其中i是折现率,n是设备寿命。这个公式很多做电力的人不太熟,但在做“投资+运行”联合优化时避不开。风电寿命通常取20年,光伏25年,电解槽10年,合成氨装置20年,每一项都要单独计算年化成本,不能笼统用同一个n。

典型日模型里还有一个非常容易踩的坑:投资年化成本是按“年”算出来的,而购电成本、卖氨收益是按“天”算的。如果直接把年值加进日优化目标里,量纲就乱了。解决办法是年化成本除以365,换算成日均投资成本。我见过有人复现论文时漏掉这个除法,目标函数里其他项都是几十万量级,投资成本一项直接几千万,优化器为了降低投资成本会把电解槽容量压到最低,结果和原文完全对不上。

做一个简单的数值估算。假设电解槽单位投资5000元/kW,折现率5%,寿命10年,则CRF约等于0.1295,每kW的年化成本约647.5元。100MW电解槽对应的年投资成本约6475万元,折算到每天约17.7万元。你可以用这类粗略计算快速验证模型输出的成本项是否在合理量级,很多明显错误的代码在这一步就能暴露。

2.2 关键约束与常见陷阱

约束条件是模型里最需要理清逻辑的部分。下面这8条是风光制氢合成氨系统优化里最核心的约束骨架,建议复现时先把它们全部写对,再根据论文的具体假设逐步扩展。

电功率平衡约束:

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

左边是供电源:风电、光伏、电网购电;右边是消耗源:电解槽输入功率、弃电功率。这条约束意味着任何时刻系统都不能缺电,多余的电只能弃掉。注意这里的P_wt和P_pv是预测出力,如果模型允许削减风电光伏(curtailment),那么进入系统的实际上是削减后的值,P_cur表示的就是被削减的部分。

电解槽运行区间约束:

δ_el(t)·P_min ≤ P_el(t) ≤ δ_el(t)·P_max,δ_el(t) = {0,1}

这条约束表达了“要么停机,要么在最低负载和额定功率之间运行”。P_min通常是额定功率的20%左右。有些论文里电解槽还允许短时过载,那是额外的假设,复现时以论文为准。

电解槽产氢与耗电关系:

H_prod(t) = P_el(t) / e_el

这里e_el是单位制氢电耗,单位是kWh/kg H2。注意这是全项目最容易埋雷的地方,详细分析放在后面踩坑章节,这里先记住一个关键点:e_el的数值是效率定义和热量基准共同决定的。

储氢罐递推约束:

S(t) = S(t-1) + H_prod(t) - H_amm(t)

这说明储氢量等于上一时刻剩余量,加上这一小时产氢,减去这一小时合成氨用氢。这是把所有时序决策联系在一起的“账本”约束,也是整个模型里最核心的动态方程。

储氢罐容量约束:

0 ≤ S(t) ≤ S_max

储氢量不能为负,也不能超过储氢罐容量。S_max的取值直接决定系统的跨时段调节能力,值得做敏感性分析。

储氢罐循环约束:

S(0) = S(T-1)

这条在典型日模型里非常重要,它要求优化结束后储氢量回到初始值,代表系统可以日复一日地滚动运行。如果不加这条,模型会白白消耗初始储氢,优化结果看起来成本很低,但实际工程中无法持续。

合成氨负荷约束:

γ·A_max ≤ H_amm(t) ≤ A_max

这里A_max是合成氨单元每小时的用氢上限,γ是最小负荷率,比如40%。这条约束表达了合成氨装置不能随意停机、不能低负荷运行的要求。部分论文会把它写成固定产量约束,也就是H_amm(t)恒等于A_max,相当于让合成氨全天满负荷运行,这是更极端的连续生产假设。

氢氨化学计量换算:

M_NH3(t) = H_amm(t) / μ

μ是生产1kg氨所需的氢气量,理论值是0.1776 kg H2/kg NH3,考虑过程损耗后通常取0.18~0.20。具体依据来自化学反应式N2 + 3H2 → 2NH3,按H2摩尔质量2.016 g/mol、NH3摩尔质量17.031 g/mol计算,3×2.016/(2×17.031)≈0.1776。这个换算系数在计算卖氨收益时直接决定收入量级,取错会导致收益偏差10%以上。

上面这些约束在代码里看起来都不复杂,但组合起来会让模型的决策变量和约束数量快速增长。24小时模型有几十个变量和上百条约束,8760小时模型则有几万个变量和几十万条约束,这时候求解器的选择就变得很重要。

2.3 一个典型日规模算例的参数设计

为了让后面的代码讲解不悬空,我给出一组常见的系统参数,你在复现时可以参考。这套参数不是从某篇特定论文摘出来的,而是综合了多篇文献里常见的中等规模配置。

项目取值
风电装机容量300 MW
光伏装机容量200 MW
电解槽额定功率150 MW
电解槽最低负载率20%
单位制氢电耗50 kWh/kg H2
储氢罐容量20 t
合成氨最大用氢速率1250 kg/h
合成氨最低负荷率40%
氨价2500 元/t
购电价分时电价,0.3~0.6 元/kWh

20 t储氢罐配合最大1250 kg/h的用氢速率,意味着即使没有光伏风电输入,系统也能靠储氢支撑大约16小时的满负荷合成氨生产,这个缓冲能力足以应对夜间的光伏缺失和风速低谷。这类参数的物理含义可以帮你判断优化结果是否合理,比如储氢罐容量如果只有5 t,系统会频繁限产;如果达到100 t,边际效益就会明显下降,这些规律都可以在代码里验证。

3. Python工程化复现:从公式到可运行代码

3.1 工程目录与参数配置:复现项目的第一步

论文复现不是把公式翻译成代码那么简单,它本质上是一个小型软件工程项目。如果第一版代码把所有数据和参数都硬编码在脚本里,那么每次跑敏感性分析都要从头找参数、改脚本,纯属浪费生命。我的做法是用下面的目录结构组织代码:

wind_solar_h2_nh3/ ├── config.yaml ├── data/ │ ├── wind_power.csv │ ├── pv_power.csv │ └── price.csv ├── src/ │ ├── __init__.py │ ├── model.py │ ├── solve.py │ └── plot.py ├── main.py └── requirements.txt

config.yaml存放所有经济参数、效率参数和设备参数;data目录存放风电、光伏、电价的时序数据;src目录里的model.py负责建模,solve.py负责求解和结果导出,plot.py负责画图;main.py是入口。这个结构的目的只有一个:改参数不用改代码。

配置文件里建议这样组织,把所有参数按设备分组,让人一眼能看出每项参数属于哪个模块:

electrolyzer: capacity_mw: 150 min_load_rate: 0.2 power_consumption_kwh_per_kg: 50 # 含辅助系统 storage: capacity_kg: 20000 ammonia: max_h2_kg_per_h: 1250 min_load_rate: 0.4 h2_per_nh3_kg: 0.18 price: nh3_yuan_per_ton: 2500 capex: electrolyzer_cost_per_kw: 5000 electrolyzer_life_years: 10 discount_rate: 0.05

写注释时把单位写清楚,这条建议怎么强调都不过分。很多人复现完一个项目隔了几个月再回头看代码,最痛苦的就是想不起来某个参数到底是kW还是MW、是元/kWh还是元/MWh,注释里写清单位能省去大量返工。

3.2 核心建模代码逐段拆解

建模部分我使用Pyomo,它是目前Python生态里最成熟的优化建模库之一,语法清晰,社区文档丰富,支持求解器的种类也多。下面给出一个核心骨架示例,以24小时典型日为优化周期。

先配置数据读取和模型对象:

import yaml import pandas as pd from pyomo.environ import * with open('config.yaml', 'r', encoding='utf-8') as f: cfg = yaml.safe_load(f) P_wt = pd.read_csv('data/wind_power.csv', header=None).values.flatten() P_pv = pd.read_csv('data/pv_power.csv', header=None).values.flatten() price = pd.read_csv('data/price.csv', header=None).values.flatten() T = len(P_wt) # 折算系数:按50 kWh/kg电耗,1 MW功率产氢量 # 1000 kW / 50 kWh/kg = 20 kg/h per MW k_h2 = 1000 / cfg['electrolyzer']['power_consumption_kwh_per_kg']

这里的k_h2很关键,它表示每MW电功率对应的产氢速率是20 kg/h。后面所有涉及“电解槽电功率→产氢量”的换算都靠它,单位必须是kg/h。

接下来定义决策变量:

m = ConcreteModel() m.T = RangeSet(0, T - 1) # 连续变量 m.P_el = Var(m.T, bounds=(0, cfg['electrolyzer']['capacity_mw']), doc='电解槽输入功率 MW') m.P_grid = Var(m.T, bounds=(0, 50), doc='电网购电 MW') m.P_cur = Var(m.T, bounds=(0, None), doc='弃电功率 MW') m.S = Var(m.T, bounds=(0, cfg['storage']['capacity_kg']), doc='储氢量 kg') m.H_amm = Var(m.T, bounds=(0, cfg['ammonia']['max_h2_kg_per_h']), doc='合成氨用氢 kg/h') # 0-1变量 m.delta_el = Var(m.T, within=Binary, doc='电解槽开停')

注意这里储氢罐容量用kg作为单位,电解槽功率用MW,购电功率也用MW。为什么不用kW统一?因为输入数据文件里风电光伏通常是MW量级,电价是元/kWh,这就会出现一个单位换算陷阱:如果要算购电成本,需要用P_grid(单位MW)乘以1000变成kW,再乘以电价和时长。这个1000倍换算很容易漏,建议在规划成本时单独处理。

然后是目标函数。这里为了控制篇幅,我聚焦在购电成本和卖氨收益的部分,年化投资成本和运维成本在代码里以注释形式给出计算思路:

def objective_rule(m): grid_cost = sum(m.P_grid[t] * 1000 * price[t] for t in m.T) # 购电成本,元 # 投资年化成本折算到日,在上层计算后作为常量传入 # inv_cost_daily = sum(capex_i * capacity_i * crf_i / 365) # 运维成本按投资比例估算 # om_cost_daily = inv_cost_annual * om_ratio / 365 nh3_revenue = (sum(m.H_amm[t] for t in m.T) / cfg['ammonia']['h2_per_nh3_kg']) / 1000 * cfg['price']['nh3_yuan_per_ton'] # 注意这里把kg换算成吨再乘氨价,得到单位元 return grid_cost + inv_cost_daily + om_cost_daily - nh3_revenue m.objective = Objective(rule=objective_rule, sense=minimize)

其中inv_cost_daily和om_cost_daily需要在外层根据参数算好再传入模型,不要在目标函数内部重复计算。nh3_revenue的换算逻辑是:H_amm加起来是kg,除以h2_per_nh3_kg得到NH3的kg数,再除以1000变成吨,乘氨价得到元。这行代码里叠加了两个换算比例,最容易出错,建议每一层都加注释。

核心约束的写法如下:

# 电功率平衡 def power_balance_rule(m, t): return P_wt[t] + P_pv[t] + m.P_grid[t] == m.P_el[t] + m.P_cur[t] m.power_balance = Constraint(m.T, rule=power_balance_rule) # 电解槽最低负载率 def el_min_load_rule(m, t): return m.P_el[t] >= (cfg['electrolyzer']['min_load_rate'] * cfg['electrolyzer']['capacity_mw

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

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

立即咨询