拿到“风光制氢合成氨系统优化”这个题目很多人的第一反应是先去找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(1i)^n / [(1i)^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 元/kWh20 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.txtconfig.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, encodingutf-8) as f: cfg yaml.safe_load(f) P_wt pd.read_csv(data/wind_power.csv, headerNone).values.flatten() P_pv pd.read_csv(data/pv_power.csv, headerNone).values.flatten() price pd.read_csv(data/price.csv, headerNone).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, withinBinary, 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(ruleobjective_rule, senseminimize)其中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, rulepower_balance_rule) # 电解槽最低负载率 def el_min_load_rule(m, t): return m.P_el[t] (cfg[electrolyzer][min_load_rate] * cfg[electrolyzer][capacity_mw
拿到“风光制氢合成氨系统优化”这个题目很多人的第一反应是先去找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(1i)^n / [(1i)^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 元/kWh20 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.txtconfig.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, encodingutf-8) as f: cfg yaml.safe_load(f) P_wt pd.read_csv(data/wind_power.csv, headerNone).values.flatten() P_pv pd.read_csv(data/pv_power.csv, headerNone).values.flatten() price pd.read_csv(data/price.csv, headerNone).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, withinBinary, 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(ruleobjective_rule, senseminimize)其中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, rulepower_balance_rule) # 电解槽最低负载率 def el_min_load_rule(m, t): return m.P_el[t] (cfg[electrolyzer][min_load_rate] * cfg[electrolyzer][capacity_mw
