搞新能源优化研究的人，十有八九都绕不开"风光制氢合成氨"这个方向。我去年拿到一个类似课题，第一步就是要把一篇相关论文里的优化模型完整复现出来，用于后续的算法对标。当时想得简单：能量守恒加设备约束，列一个目标函数交给求解器跑就完事了。真正动手才发现，光是搞懂论文里模型假设和公式的映射关系，就花了大半个月，后面调约束、改参数又是一堆凌晨。

这篇文章不打算逐行讲某篇具体论文的代码，而是把"拿到这样一类论文，如何用Python把优化模型完整复现出来"的整个路线讲透。内容包含系统建模思路、求解器选型、数据场景生成、Pyomo代码实现、常见陷阱与调试技巧，最后聊一下从确定性模型往随机优化扩展的思路。如果你正在做新能源制氢、综合能源系统优化，或者需要复现类似的论文、扩展自己的模型框架，这篇文章应该能帮你把最痛苦的摸索期压缩掉很大一截。

1. 系统建模：先想清楚"优化什么"

1.1 风光制氢合成氨的系统架构与能量流

风光制氢合成氨，本质上是把不稳定的风光电转成稳定、可储存、可运输的化学品氨。系统一般由这么几块组成：风电场和光伏电站作为一次电源；电解水制氢单元把电能转成氢能；储氢装置用于缓冲氢气生产和下游氨合成之间的供需波动；空分单元提供氮气；合成氨装置在高温高压下完成氢氮反应生成氨；最后还有储氨和外运环节。很多研究还会考虑与电网交互，也就是允许系统从电网购电，甚至在某些时段向电网反送电，这部分会影响目标函数和约束条件的写法。

建模之前，我强烈建议先画一张能量流图，不要觉得这是浪费时间。这张图上，电源侧、变换侧、储能侧、负荷侧分得越清楚，后面变量定义就越顺。我自己复现论文时，先根据摘要里的系统简图画了一版，又对着正文里的公式补了一版，最终把变量归成三类：

• 容量配置变量：风机装机容量、光伏装机容量、电解槽额定功率、储氢罐容量、合成氨装置额定产量等。
• 运行调度变量：各时段风电/光伏出力、电解槽输入电功率、产氢速率、储氢罐充放氢速率、合成氨产量、购电功率、弃电功率等。
• 逻辑变量：设备启停状态、容量是否建设的0/1决策等。

把变量分类之后，再去看论文的目标函数和约束，你会发现思路清晰得多。这类论文的数学模型通常就是围绕这三类变量展开的，只是有的论文用单层优化，有的用配置层和运行层嵌套的双层规划，表述形式略有不同。

1.2 优化目标怎么定：成本、收益还是综合指标

风光制氢合成氨系统的优化目标，看过的论文里基本逃不出三类：

第一类是总投资和运行成本最小化。设备投资按等额年值折算，加上每年的运维成本、购电成本，构成年化总成本。这类目标常用在规划问题里，输出结果是一个最优容量配置方案。

第二类是项目净收益最大化。收益端包括售氨收入、可能的售电收入，成本端包括投资、运维、购电。这类目标更偏运营优化，输出的是给定容量下的调度策略。

第三类是综合目标或多目标优化。比如把可再生能源弃用率、碳排放成本引入目标函数，或者用加权法、ε-约束法把经济性和环保性同时纳入优化。这类模型求解复杂度更高，但更贴近工程实际。

复现时要特别注意论文里写的是"总成本最低"还是"净收益最大"，这两个写法在目标函数符号和约束处理上很容易搞反。我自己一开始就吃过亏，目标函数少加了一项储氢溢出惩罚，结果求出来的电解槽容量比论文结果大了将近一倍。后来逐行比对才发现是惩罚项漏掉了。

很多论文还会用平准化成本作为最终评价指标，公式大概是：

LCOA = 年化总成本 / 年合成氨产量

这个指标用来横向比较不同配置方案很方便，但注意复现时年化成本里是否包含储氢、储氨的投资，各篇论文口径不一样，直接比较时要先统一口径。

1.3 约束条件：这些细节不能漏

优化模型里约束条件是最容易出问题的地方。我复现时会把约束分成五组来写，每写一组就回去对一下论文公式：

功率平衡是最基础的约束，任意时段进入系统的电功率必须等于消耗的功率。如果允许弃电，要显式加入弃电变量，不然求解器为了让平衡成立可能强行压电解槽出力，结果不符合实际。

电解槽约束里有几个容易被坑的点。很多电解槽有最小稳定运行负载，比如额定功率的20%，低于这个值电解槽不能运行。这个约束通常用二进制变量表示：电解槽开启时，输入功率必须在最小负载和额定功率之间。如果漏掉这一条，模型可能会让电解槽在很低功率下运行，产氢成本异常低，优化结果和论文对不上。

储氢动态约束是连接电解槽和氨合成单元的桥梁。需要注意单位统一，产氢速率、储氢量、氨合成耗氢速率，三者单位必须一致。我第一次写模型时产氢用了kg/h，储氢用了MWh，氨合成耗氢用了吨，结果求解器报不可行，查了半天才意识到是单位问题。另外很多论文要求调度周期的初始储氢量和终止储氢量相等，这是为了保证调度方案的循环可行性，别忘了加。

合成氨部分的核心是化学计量关系。合成氨反应的理想化学计量是：

N₂ + 3H₂ → 2NH₃

按分子量换算，生产1吨氨大约需要消耗176公斤氢气。这个系数直接决定氨产量和耗氢量的换算关系，写错整个物料平衡就乱套了。

2. 模型求解路线与代码架构设计

2.1 为什么用Python：从MATLAB搬过来的真实感受

如果你以前用MATLAB做优化，转Python需要一点适应期，但一旦上手你会觉得回不去。论文复现不是一次性工程，过程中要反复换参数、换数据、做对比实验。Python的pandas加numpy组合处理时序数据非常顺手，读CSV、透视、重采样都是几行代码的事。Pyomo这类建模工具写出来的约束和变量跟数学公式几乎一一对应，调试时可以非常直观地检查模型结构。

另外Python社区优化生态非常全。线性规划和混合整数规划有多个开源求解器，非线性有IPOPT，启发式算法库也很丰富。即使论文里用的是很冷门的模型形式，基本也能找到对应的参考实现。我复现时用的建模框架是Pyomo，选择的理由是它比较通用，换求解器只需要改一个字符串，不需要动模型定义代码。

2.2 求解器选型：先判断问题类型再选求解器

先判断你的模型是什么类型。风光制氢合成氨优化模型常见的形态是混合整数线性规划，也就是MILP。为什么不是非线性？因为论文为了求解方便，通常会把电解槽的效率曲线做分段线性化，把储氢罐的压力-容量关系简化成线性形式，合成氨单元的运行区间也用线性约束表达。所以不要一上来就想着非线性求解器，优先把模型线性化，能用MILP就不要用MINLP，求解速度和稳定性差了一个量级。

常用求解器我整理了一个对比表：

就个人复现而言，开源求解器完全够用。我复现的那个论文模型属于中等规模，大概有几千条约束和几百个变量，用HiGHS跑几十秒就出结果。如果模型规模大、场景数量多，商业求解器会更省心，但先别急着上商业执照，把开源求解器用熟再说。

2.3 工程目录与环境配置

工程目录建议一开始就按可复现的方式组织，别把所有脚本塞进一个文件。我常用的目录结构如下：

project/ ├── README.md ├── requirements.txt ├── data/ │ ├── raw/ # 原始风/光数据 │ └── processed/ # 清洗后的数据 ├── src/ │ ├── __init__.py │ ├── data_prep.py # 数据加载与场景生成 │ ├── model.py # Pyomo模型定义 │ ├── solver.py # 求解与结果导出 │ └── plot.py # 可视化 ├── results/ └── configs/ └── default.yaml

环境依赖建议用requirements.txt锁版本。复现论文这种事，版本浮动有时候会让求解结果出现细微差别，对严肃的对比实验来说是致命的。

numpy>=1.21 pandas>=1.3 pyomo>=6.4 highspy>=1.5 matplotlib>=3.5 pyyaml>=6.0

pip安装好之后，Pyomo负责建模，highspy是HiGHS求解器的Python接口。把依赖固定住，剩下的时间全部花在建模和调试上。

3. 核心代码实现与关键参数

3.1 数据准备：没有论文数据时怎么生成场景

复现论文最卡人的往往不是模型，而是数据。很多论文不会公开风电场和光伏电站的逐时出力数据，这时候有两条路。一是找公开数据集，比如气象再分析数据或者一些开放的风光实测序列；二是根据论文里对资源条件的描述，用统计分布模拟出力曲线。

我当时用了模拟加局部敏感性验证的方法。风电按威布尔分布采样，光伏按白天正弦曲线叠加噪声，再加入一阶延迟模拟时序相关性。这样生成的数据虽然不能和论文原数据完全一致，但用于验证模型逻辑和结果趋势是足够的。

import numpy as np import pandas as pd np.random.seed(2024) hours = 24 * 365 # 一年逐时数据 time = pd.date_range("2023-01-01", periods=hours, freq="h") # 风电：威布尔分布采样后折算功率 wind_speed = np.random.weibull(2.2, hours) * 6.5 wind_capacity = 300 # MW wind_power = wind_capacity * np.clip(wind_speed / 12, 0, 1) # 光伏：正弦曲线加上噪声 pv_capacity = 200 # MW hour_of_day = time.hour.to_numpy() pv_power_raw = np.sin(np.pi * (hour_of_day - 6) / 12) pv_power = pv_capacity * np.clip(pv_power_raw, 0, 1) pv_power += np.random.normal(0, 8, hours) pv_power = np.clip(pv_power, 0, pv_capacity) df = pd.DataFrame({"wind": wind_power, "pv": pv_power}, index=time)

注意，模拟数据只用来调试模型逻辑，不能替代理论文里的算例验证。如果要做严肃的定量对标，还是要尽量还原论文使用的典型日或典型场景数据。

3.2 优化模型构建：Pyomo代码逐步说明

下面给出一个简化的24小时调度模型框架。它省略了一些工程细节，但骨架是正确的，后续可以往里加约束。先定义时段、参数和决策变量：

from pyomo.environ import * model = ConcreteModel() # 时段 T = range(24) # 示意数据：典型日风电/光伏出力（单位MW） wind_in = [110.2, 105.8, 98.5, 91.0, 85.2, 82.3, 88.4, 101.2, 118.5, 130.2, 141.0, 137.6, 125.4, 112.3, 105.6, 98.4, 90.1, 88.2, 86.5, 89.3, 95.7, 102.4, 108.9, 115.0] pv_in = [0, 0, 0, 0, 0, 0, 18.5, 62.5, 112.3, 155.2, 175.6, 183.4, 178.9, 155.6, 128.3, 86.5, 45.2, 12.3, 0, 0, 0, 0, 0, 0] # 决策变量 model.p_el = Var(T, within=NonNegativeReals, bounds=(0, 150), doc="电解槽输入功率(MW)") model.h2_sto = Var(T, within=NonNegativeReals, bounds=(0, 50000), doc="储氢量(kg)") model.f_h2 = Var(T, within=NonNegativeReals, doc="产氢速率(kg/h)") model.f_am = Var(T, within=NonNegativeReals, doc="合成氨产量(t/h)") model.p_grid = Var(T, within=NonNegativeReals, doc="电网购电(MW)") model.y_am = Var(T, within=Binary, doc="合成氨装置启停")

这里我故意把变量注释写得清楚，工程上建议每个变量都标注单位。接下来写目标函数和约束：

# 目标：最小化购电成本 + 储氢溢出惩罚（简化示例） def objective(model): return sum(0.5 * model.p_grid[t] for t in T) \ + 1e-3 * sum(model.h2_sto[t] for t in T) model.obj = Objective(rule=objective, sense=minimize) # 功率平衡约束 def power_balance(model, t): return model.p_el[t] + model.p_grid[t] == wind_in[t] + pv_in[t] model.power_bal = Constraint(T, rule=power_balance) # 电解槽产氢关系：按60kWh/kgH2计算 def hydrogen_prod(model, t): return model.f_h2[t] == 0.05 * model.p_el[t] model.h2_con = Constraint(T, rule=hydrogen_prod) # 合成氨装置产氨下限（大M法表示启停联动） M = 1000 def am_min_output(model, t): return model.f_am[t] >= 10 * model.y_am[t] - (1 - model.y_am[t]) * M model.am_min = Constraint(T, rule=am_min_output) # 储氢动态平衡：1吨氨约需176kg氢 def storage_dynamic(model, t): if t == 0: return model.h2_sto[t] == 0 return model.h2_sto[t] == model.h2_sto[t-1] \ + model.f_h2[t-1] - 176 * model.f_am[t-1] model.sto_dyn = Constraint(T, rule=storage_dynamic) solver = SolverFactory("highs") results = solver.solve(model, tee=False) print("求解状态:", results.solver.termination_condition)

有几个地方要特别提醒。第一，储氢量、产氢速率、耗氢速率三者单位必须一致。上面代码里统一使用kg和kg/h，避免了单位混乱，这也是最省心的做法。第二，合成氨产量f_am单位是吨每小时，乘以176才是对应的耗氢质量流量，这个换算来自化学计量关系。第三，大M法里的M取值不能太小，要大于约束可能出现的最大值，用1000在这个问题里是安全的。

3.3 求解与后处理：结果可视化

求解只是第一步，结果可视化是论文复现工作量的大头。通常要画系统功率分配堆叠图、储氢量与氨产量时序图、成本结构图。我一般优先画前三张：

import matplotlib.pyplot as plt # 取出求解结果 p_el_val = [value(model.p_el[t]) for t in T] p_grid_val = [value(model.p_grid[t]) for t in T] h2_sto_val = [value(model.h2_sto[t]) for t in T] f_am_val = [value(model.f_am[t]) for t in T] fig, axes = plt.subplots(3, 1, figsize=(10, 8), sharex=True) # 功率平衡 axes[0].bar(T, wind_in, label="风电出力", color="steelblue") axes[0].bar(T, pv_in, bottom=wind_in, label="光伏出力", color="gold") axes[0].bar(T, p_grid_val, label="电网购电", color="lightcoral") axes[0].set_ylabel("功率 / MW") axes[0].legend() # 储氢量 axes[1].plot(T, h2_sto_val, marker="o", label="储氢量") axes[1].set_ylabel("储氢量 / kg") axes[1].legend() # 合成氨产量 axes[2].bar(T, f_am_val, label="氨产量", color="mediumseagreen") axes[2].set_ylabel("氨产量 / t") axes[2].set_xlabel("时段 / h") axes[2].legend() plt.tight_layout() plt.savefig("results/scheduling_result.png", dpi=300)

调试心得：当求解结果出现振荡或者异常跳变时，先别急着怀疑目标函数错了，去看功率平衡子图。如果风电和光伏完全消纳、购电没有突变，说明电平衡没问题；再往下看储氢量和氨产量曲线，通常能快速定位到是哪条约束在捣乱。很多问题一眼就能从图上判断出来。

3.4 与论文结果对标：先跑通再谈优化

复现完成的第一标志不是"模型能求解"，而是"结果趋势能和论文基本对上"。我常用的做法是固定论文给出的容量配置，只做运行调度优化，验证起止储氢量、典型日氨产量、弃电率这些关键指标是否匹配。如果匹配，再放开容量变量做全优化。这样做的好处是，把配置层和运行层解耦，出问题的时候能快速定位是容量变量导致的偏差，还是运行约束导致的偏差。

如果结果和论文差异很大，先检查三个方面：数据场景是否一致、单位是否一致、论文目标函数里是否包含你没有注意到的惩罚项。优先级是这样排序的：模型逻辑错误最严重，其次单位错误，再次数据场景差异，最后才是求解器容差。

4. 复现过程中踩过的坑与排查技巧

4.1 常见问题速查表

我复现过程中遇到的典型问题，整理成了下面这个表，基本覆盖了大部分新手会碰到的坑：

最重要的排查方法只有一个：把模型导出成LP格式，用文本编辑器打开，逐条对照数学公式。对，就是最土的方法，但效率极高。

model.write("model.lp", io_options={"symbolic_solver_labels": True})

这条命令会生成一个可读的LP文件，变量和约束都带名字。我用这一招校出过无数低级错误，强烈推荐。

4.2 非线性约束线性化的三种典型操作

不是所有模型都要线性化，但如果你想用开源求解器稳定求解大规模MILP，线性化几乎是绕不开的。三处最常见的线性化处理：

第一种是设备效率曲线的分段线性化。电解槽的电耗和产氢量不是严格线性关系，在论文里经常被写成非线性曲线。复现时可以把它拆成若干段线性函数，引入二进制变量表示当前运行在哪一段。图解法很简单，把横轴划分成区间，每个区间用一条直线近似，然后加一组大M约束激活对应的线段。

第二种是if-then逻辑。比如"合成氨装置如果开启，最低产氨量不得低于10吨每小时"。这类约束直接用大M法处理。上面代码里的am_min_output就是这么写的。y等于1时约束生效，y等于0时约束自动松弛。

第三种是最大化最小值或最小化最大值这类目标。比如"最小化最大购电功率"，引入辅助变量gamma，加约束p_grid[t]小于等于gamma，目标函数里放gamma。这样一来，非线性目标被转换成线性约束加线性目标，模型性质完全变了。

复现论文时如果发现求解器报非线性或不收敛，先看看目标函数和约束里有没有绝对值、最大值、乘积项。这几个关键词出现任何一个，都要考虑做线性化处理。

4.3 求解慢、不收敛怎么办

求解性能优化的策略，按性价比排序主要有五条：

先跑确定性模型，再跑随机或鲁棒模型。场景一多，问题规模指数级上升。我刚接触两阶段随机优化时，直接把100个场景丢进去，求解器跑了几个小时都不出结果，后来先把场景削减到12个才跑通。

给变量设置紧的上下界。特别是储氢罐初始储量、合成氨日产量上限这些，给一个贴近实际的界能极大缩小搜索空间。求解器很多时候会浪费大量时间探索无意义的解空间。

尝试热启动。把上一个求解结果作为初始解传入模型，对MILP的启停变量特别有效。

对设备启停变量，可以用对称性破缺约束减少等价分支。比如多个电解槽如果完全相同，可以加约束让它们的启停顺序固定，避免求解器在等价解之间反复切换。

如果以上都不行，就直接用商业求解器。个人做研究不用纠结，开源求解器够用就用开源的，规模大了再换商业的，效率差距确实明显。

我复现项目里有一次就是储氢约束过多导致求解时间爆炸。后来把全年8760个时段削减成12个典型场景,给每个场景分配权重，求解时间从数小时降到分钟级，优化目标变化不到3%。从此我养成了一个习惯：凡是大规模时序优化，先做场景削减再求解。

5. 扩展思路：从确定性优化走向不确定性优化

5.1 灵敏度分析：模型很快能找到关键瓶颈

模型跑通之后，第一步建议做灵敏度分析。方法很简单，固定其他参数不变，逐个改变关键参数，比如电解槽投资单价、上网电价、氨销售价格、风光资源系数，看优化结果如何变化。这一步能快速判断系统的瓶颈在哪里。

我做灵敏度分析时发现，电解槽投资成本对系统年化成本的影响最大，其次是年风光有效利用小时数。这其实符合直觉：风光制氢合成氨系统中，电解槽是整个能量转换链条的核心，它的投资和运行特性直接决定全系统的经济性。氨合成单元的产能约束在某些场景下会成为瓶颈，这时需要增加储氢容量来平滑氢气供需波动。

灵敏度分析的结果不仅能验证模型合理性，还能为后续论文写作或项目决策提供素材。很多论文的结论部分其实就是在讲参数变化如何影响最优方案。

5.2 两阶段随机规划：论文进阶的常见方向

确定性模型的局限在于，它假设风光出力是已知的。实际工程中，风光出力具有很大的不确定性。论文复现做到一定程度，你大概率会想往随机优化方向扩展。

两阶段随机规划的经典思路是把决策分成两阶段。第一阶段在不确定性实现之前做容量配置决策，确定装多少风机、多少光伏、多大电解槽。第二阶段在不确定性实现后做运行调度决策，针对每个随机场景给出最优调度方案，目标函数里是所有场景期望成本的总和。

# 两阶段随机规划的目标函数示意 # 第一阶段成本 + sum(场景概率 * 第二阶段运行成本) model.obj = Objective( expr=first_stage_cost + sum(prob[s] * second_stage_cost[s] for s in SCENARIOS), sense=minimize )

实现时要注意，第二阶段变量必须为每个场景单独定义，第一阶段变量跨场景共享。这个结构决定了随机规划模型的规模会迅速膨胀，所以场景削减非常关键。常见的场景削减方法有基于概率距离的同步回代消除法、聚类方法等。削减后的场景集要在保持原始概率分布特征的前提下，尽量压缩数量。

5.3 多目标优化的处理思路

如果你面对的论文涉及经济性和可再生能源利用率两个目标，处理方案通常是两种。

第一种是加权求和法，给不同目标乘以权重因子并相加。简单直观但需要反复试权重，且难以处理量纲差异。第二种是ε-约束法，把其中一个目标降为约束，要求它不大于某个阈值，然后优化另一个目标，遍历不同阈值生成帕累托前沿。

我复现过的模型在经济目标和弃光弃风率之间做权衡时，采用ε-约束法效果很好。把弃光弃风率上限依次设为30%、25%、20%、15%，分别求解经济最优，最后能画出一条平滑的帕累托前沿曲线，给决策者提供了非常直观的参考。多目标优化的结果分析和可视化又是一块单独的工作，画出来的图在论文里会很加分。

最后说点实在的

写这篇文章的时候，我想到自己第一次完整跑通风光制氢合成氨优化模型时的情形。盯着功率分配堆叠图看了好久，储氢曲线是平滑的，氨产量稳定在额定区间附近，所有约束都对上了，那一刻确实挺解压的。

总结下来，复现这类优化论文，核心就三件事：把系统物理架构画明白，把数学模型公式写严谨，把求解边界条件调合理。真正的功夫都花在数据预处理和约束调试上，求解器反而是最省心的环节。

最后再分享一个小技巧：如果论文里没有明确的约束细节，可以先用简化模型跑通流程，然后逐步增加约束的复杂程度。每加一组约束，就记录一下结果和目标值的变化。这样做的好处是，出问题的时候你能清楚知道是哪一组新约束引入了不可行性。这个小习惯帮我省掉了大量返工时间。后面如果你也在做类似的复现工作，欢迎一起交流建模和调试中的经验。