前阵子帮朋友复现了一项关于并网/离网风光互补制氢合成氨系统的容量-调度优化研究源码是用 Python 写的。模型本身不算特别复杂但当“容量规划”和“小时级调度”放在同一个优化问题里求解时坑比我想象中多很多。折腾了一周多索性把整个复现过程、数学模型搭建思路、代码实现细节和踩过的坑整理成文给同样在做新能源制氢、合成氨方向或者想复现类似论文的同学做个参考。先说清楚这个项目到底在解决什么问题。所谓“并网/离网风光互补制氢合成氨系统”简单说就是风电场和光伏电站同时给电解槽供电电解水制氢后氢气与氮气在一定温度和压力下反应生成氨。并网模式下系统可以和电网双向买卖电离网模式下则是完全孤岛运行所有电力都来自风光。容量-调度优化则是要同时回答两个问题一是长期投资决策风光装机、电解槽功率、储氢罐容量、合成氨装置规模该建多大二是短期运行决策每个小时风电、光伏出力怎么分配给电解槽、储氢罐怎么充放、合成氨装置运行在什么负荷。这两个问题互相耦合不能分开求解这也是整个项目的核心难点。下面我就按“问题拆解 → 数学建模 → Python实现 → 排坑经验 → 扩展思考”的顺序来写尽量把复现过程中那些论文里没说透的细节补全。1. 项目背景与核心需求解析1.1 风光互补制氢合成氨系统到底长什么样很多人一听“制氢合成氨”就觉得是个大型化工项目确实如此。从能源供给侧看有风力发电机组和光伏阵列分别输出变动功率从负荷侧看主要有三块电解水制氢装置、合成氨装置以及必要的辅助设备。中间还有储氢装置作为缓冲因为风光出力波动很大而合成氨装置更希望稳定进气。如果系统还并网电网就成了一个“无限大”的缓冲池可以买电、卖电离网状态下则只能靠储氢、储能来过平衡可靠性压力一下就上来了。用个生活化的类比这套系统就像一个自成一体的微型“食材加工厂”风机和光伏是“菜园”什么季节收什么菜产量不稳定电解槽是“清洗切配间”负责把水变成氢气和氧气储氢罐是“冷库”把暂时用不完的氢气存起来合成氨装置是“中央厨房”按配方氮氢比3比1连续产出氨。并网模式相当于厨房可以随时从外面超市买菜买电也能把多余的菜卖掉卖电换钱离网模式则只能靠自己的菜园和冷库应付。合成氨之所以在这个系统里这么有吸引力是因为氨本身就是大宗化工原料也是理想的氢能载体——液氨储运比高压氢容易能量密度也不低。很多可再生能源丰富但消纳困难的区域都在规划“风光发电—电解水制氢—合成氨”这条产业链。1.2 容量-调度耦合决策一个双向影响的问题新手拿到这个题目时第一反应往往是“先算最优容量再算最优调度”拆成两个步骤。但实际做下来会发现这两层决策是强耦合的。举个例子如果光伏装得多中午时段发电严重过剩离网模式下如果不配置足够大的电解槽或者储氢罐来消纳那这部分电力就只能白白弃掉最终影响项目经济性。反过来电解槽容量选大了在夜间风电低、光伏为零的时候又可能因为电力不足而长期低负荷运行调度的灵活性也不够。所以必须做“容量-调度一体化优化”即把长期投资变量和全年逐小时运行变量放进同一个模型里一次性算出同时满足建设约束和运行约束的最优方案。这样做的好处有两个一是避免分层优化导致的局部最优二是可以直接看到装机容量和运行策略之间的边际关系——比如多装1兆瓦风电能在多大程度上降低电解槽购电成本。这个问题的规模也会让初学者吃惊。假设典型日只有24个时段一年12个典型日每个典型日又有几十个连续变量和整数变量加上建设容量变量模型规模通常会达到几万个决策变量和约束。如果直接用全年8760小时逐时刻建模变量数更多求解时间可能长达数小时甚至数天。因此复现时大多会选择“典型日集合”来近似全年这也是论文里最常见的做法。1.3 为什么这个“复现”值得做以前读文献时经常看到“系统容量优化模型”这一章公式写得漂亮但想照着实现却发现很多细节被省略了。比如电解槽效率到底随负荷怎么变化储氢罐容量怎么和运行约束关联合成氨装置最小负荷是多少都是不写清楚。复现一篇相关研究的价值恰恰在于把这些坑都趟一遍。这份记录适合三类人一是研究新能源制氢或电力系统优化想快速入门建模和求解工具的研究生二是做可再生合成燃料或绿氨可行性评估的工程人员需要一套能落地的计算框架三是用 Python 做能源系统优化、但经常被 Pyomo/Gurobi 折磨的开发者。我会把思路和关键代码片段都放在下面你可以直接改参数、换数据跑自己的场景。2. 数学模型与优化框架拆解2.1 系统拓扑与设备建模先画清楚边界。我复现的模型里主要包含以下设备和变量流风电机组WT输入为风速时间序列输出为有功功率决策变量是装机台数或容量kW。光伏阵列PV输入为水平面或倾斜面辐照度输出为有功功率决策变量是安装容量kW。电解槽EL输入为电力输出为氢气流量kg/h决策变量是额定功率kW运行变量是每时段实际输入功率。储氢罐HST氢气缓冲决策变量是储氢容量kg运行变量是每个时段的充放氢量kg/h和储氢量状态kg。合成氨装置NH3 Plant输入为氢气和氮气输出为液氨t/h决策变量是额定产能kg/h运行变量是每时段负荷比。电网接口Grid仅并网模式存在运行时每时段可买电、卖电购售电价可以不同。可能还有储能电池BESS用于短时功率平衡决策变量是容量kWh和每个时段的充放电功率。下表整理了设备输入输出和关键参数方便后面建模对照设备输入输出关键建模参数风电机组风速序列电功率切入/额定/切出风速、功率曲线、单位容量投资光伏阵列辐照度电功率光伏转换效率、温度折减、单位容量投资电解槽电功率、水氢气额定功率、最低部分负荷率、电耗/kWh/kgH2储氢罐氢气氢气最大储氢量、初始/终末储量、充放速率上限合成氨装置氢气、氮气氨额定产能、最低负荷率、吨氨耗氢量、单位产能投资电网电功率电功率购电价、售电价、联络线功率上限2.2 决策变量与主要约束把决策变量分成两类在代码里也是这么组织的。容量变量(C_{WT}, C_{PV})风电、光伏装机容量kW。(C_{EL})电解槽额定功率kW。(C_{HST})储氢罐容量kg。(C_{NH3})合成氨装置额定产能kg/h。调度变量(P_{WT}(t), P_{PV}(t))风电、光伏实际计入出力kW。(P_{EL}(t), H_{EL}(t))电解槽输入功率和产氢流量。(H_{HST,in}(t), H_{HST,out}(t), S_{HST}(t))储氢罐充氢、放氢和储氢量。(H_{NH3}(t), A_{NH3}(t))合成氨装置用氢流量和产氨流量。(P_{grid,buy}(t), P_{grid,sell}(t))并网购电、售电功率。约束可以分为四组。功率平衡约束并网模式下[ P_{WT}(t) P_{PV}(t) P_{grid,buy}(t) P_{BESS,dis}(t) P_{EL}(t) P_{aux}(t) P_{grid,sell}(t) P_{BESS,ch}(t) ]离网模式下电网功率项为零。这里的辅助负荷 (P_{aux}(t)) 可以包含合成氨装置的用电或者按固定比例估。氢气物料平衡约束[ H_{EL}(t) H_{HST,out}(t) H_{NH3}(t) H_{HST,in}(t) ]再加上储氢罐的动态方程[ S_{HST}(t1) S_{HST}(t) \eta_{HST} H_{HST,in}(t) - \frac{H_{HST,out}(t)}{\eta_{HST}} ]以及容量边界 (0 \le S_{HST}(t) \le C_{HST})。设备运行约束电解槽的输入功率有上下限[ \lambda_{EL}^{min} C_{EL} \le P_{EL}(t) \le C_{EL} ]合成氨装置的用氢量也有负荷区间[ \lambda_{NH3}^{min} C_{NH3} \cdot r_{H2/NH3} \le H_{NH3}(t) \le C_{NH3} \cdot r_{H2/NH3} ]其中 (r_{H2/NH3}) 是吨氨耗氢比。年度产量与储能终值约束合成氨系统通常有年产量目标。可以设置[ \sum_{t} A_{NH3}(t) \ge A_{target} ]或者设定日产量下限。同时为了不让模型在期末把储氢罐“掏空”要加储氢罐终值等于初始值的约束。2.3 目标函数与成本模型目标函数一般是“年化总成本最小”。年化总成本由几部分构成折年投资成本将设备初投资乘以资本回收因子 (CRF) 分摊到每一年。运行维护成本按固定比例或按发电量/运行小时数折算。购电成本并网模式减卖电收益。副产品收益如果氧气也能卖可以加一个负项。写成公式[ \min \quad \sum_{i \in {WT,PV,EL,HST,NH3}} CRF_i \cdot inv_i \cdot C_i \sum_{t} \left( C_{OM} price_{buy} P_{grid,buy}(t) - price_{sell} P_{grid,sell}(t) \right) ]资本回收因子一般取 (CRF \frac{r(1r)^n}{(1r)^n - 1})(r) 是贴现率(n) 是设备寿命。这里有个容易被忽略的细节电解槽寿命往往比风电场短单位投资和运维费率也高合成氨装置的寿命最长。如果把所有设备的寿命统一取成20年得出的年化成本会和实际偏差很大。所以我复现时是分设备取寿命、分开计算折年系数。2.4 求解策略从MINLP到MILP上面这个模型直接写出来由于存在效率随负荷变化的非线性项、储氢量和功率的乘积项它是一个混合整数非线性规划MINLP。如果直接用非线性求解器一方面不一定能保证全局最优另一方面求解速度慢代码调试也很痛苦。主流的做法是线性化把模型等价转成混合整数线性规划MILP。三种最常用的手段分段线性化用于电解槽电耗-产氢曲线、风电功率曲线等。例如把电解槽负荷区间切成5段每段用线性关系近似。Big-M法用于把“设备启停”“是否建设”这类逻辑关系转化成线性不等式。McCormick包络用于线性化两个连续变量的乘积比如储氢量与充放效率的乘积。不过这个模型里我尽量避免直接出现这类乘积而是把效率放进平衡约束的系数里。我最后用的是 Pyomo 建模 Gurobi 求解 MILP。Gurobi 对大规模 MILP 的求解能力非常强如果你没有许可证也可以用开源的 CBC、HiGHS 或 SCIP但求解大案例时速度差异会非常明显。3. Python程序架构与核心实现3.1 工程文件组织复现这种项目代码文件千万别摊成一大坨。我习惯按“数据—模型—求解—后处理”四层组织hydrogen_nh3_optimization/ │ ├── data/ │ ├── wind_speed.csv │ ├── solar_irradiance.csv │ ├── electricity_price.csv │ └── scenarios.pkl │ ├── src/ │ ├── data_process.py │ ├── model.py │ ├── solve.py │ └── plot_results.py │ ├── output/ │ ├── capacity_results.csv │ ├── dispatch_results.csv │ └── figures/*.png │ ├── requirements.txt └── main.pymain.py只是编排入口读取参数后依次调用数据预处理、建模、求解、结果导出。这样别人拿到你的代码不需要翻遍所有文件就知道运行流程。3.2 数据准备构造典型日场景原始数据是逐小时的风速、辐照。全年8760个点直接优化模型太大。所以我先用 k-means 聚类得到若干个典型日并统计每个典型日出现的天数权重。这一步很像视频压缩——不逐帧保存所有信息而是提取关键帧。from sklearn.cluster import KMeans import pandas as pd # df_raw: index为时间列为 [wind_speed, solar_irradiance, price] # 按天重采样成24维向量 days [] day_info [] for day, group in df_raw.groupby(df_raw.index.date): vec group[[wind_speed, solar_irradiance]].values.flatten() days.append(vec) day_info.append(day) days pd.DataFrame(days) # 聚成12个典型日 kmeans KMeans(n_clusters12, random_state42) labels kmeans.fit_predict(days) # 统计每个典型日的权重该簇天数/总天数 weight pd.Series(labels).value_counts(normalizeTrue).sort_index()聚类完成后数据集规模从8760小时缩减为 12×24 288 小时同时每个典型日附带一个权重系数用于目标函数里计算全年成本。3.3 用 Pyomo 构建容量-调度一体化模型Pyomo 建模思路很直白先建一个ConcreteModel然后添加集合Set、参数Param、变量Var再用Constraint添加约束最后用Objective定义目标函数交给求解器。我在模型里定义了三个集合SCEN是典型日TIME是24小时DEVICE是设备名称。import pyomo.environ as pyo m pyo.ConcreteModel() m.SCEN pyo.Set(initializescenario_ids) # 典型日编号 m.TIME pyo.Set(initializerange(24)) # 小时 # 容量变量 m.C_wt pyo.Var(domainpyo.NonNegativeReals, bounds(0, ub_wt)) m.C_pv pyo.Var(domainpyo.NonNegativeReals, bounds(0, ub_pv)) m.C_el pyo.Var(domainpyo.NonNegativeReals, bounds(0, ub_el)) m.C_hst pyo.Var(domainpyo.NonNegativeReals, bounds(0, ub_hst)) m.C_nh3 pyo.Var(domainpyo.NonNegativeReals, bounds(0, ub_nh3)) # 调度变量 m.P_wt pyo.Var(m.SCEN, m.TIME, domainpyo.NonNegativeReals) m.P_pv pyo.Var(m.SCEN, m.TIME, domainpyo.NonNegativeReals) m.P_el pyo.Var(m.SCEN, m.TIME, domainpyo.NonNegativeReals) m.H_el pyo.Var(m.SCEN, m.TIME, domainpyo.NonNegativeReals) m.S_hst pyo.Var(m.SCEN, m.TIME, domainpyo.NonNegativeReals) m.H_in pyo.Var(m.SCEN, m.TIME, domainpyo.NonNegativeReals) m.H_out pyo.Var(m.SCEN, m.TIME, domainpyo.NonNegativeReals) m.H_nh3 pyo.Var(m.SCEN, m.TIME, domainpyo.NonNegativeReals) m.P_grid_buy pyo.Var(m.SCEN, m.TIME, domainpyo.NonNegativeReals) m.P_grid_sell pyo.Var(m.SCEN, m.TIME, domainpyo.NonNegativeReals)然后添加关键约束。功率平衡和氢平衡是核心。def power_balance_rule(m, s, t): return ( m.P_wt[s, t] m.P_pv[s, t] m.P_grid_buy[s, t] m.P_el[s, t] m.P_grid_sell[s, t] p_aux ) m.power_balance pyo.Constraint(m.SCEN, m.TIME, rulepower_balance_rule) def hydrogen_balance_rule(m, s, t): return m.H_el[s, t] m.H_out[s, t] m.H_nh3[s, t] m.H_in[s, t] m.hydrogen_balance pyo.Constraint(m.SCEN, m.TIME, rulehydrogen_balance_rule)电解槽的运行区间约束要同时关联容量变量和调度变量否则模型会把电解槽功率建得很大却很少开机。def el_load_rule(m, s, t): return m.el_min_factor * m.C_el m.P_el[s, t] m.el_load pyo.Constraint(m.SCEN, m.TIME, ruleel_load_rule)合成氨装置的用氢量约束def nh3_load_rule(m, s, t): return m.H_nh3[s, t] m.C_nh3 * h2_per_nh3 m.nh3_max_load pyo.Constraint(m.SCEN, m.TIME, rulenh3_load_rule)目标函数按典型日权重累加运行成本def obj_rule(m): invest_cost ( inv_wt * m.C_wt inv_pv * m.C_pv inv_el * m.C_el inv_hst * m.C_hst inv_nh3 * m.C_nh3 ) * CRF run_cost sum( weight[s] * ( om_cost price_buy[s, t] * m.P_grid_buy[s, t] - price_sell[s, t] * m.P_grid_sell[s, t] ) for s in m.SCEN for t in m.TIME ) return invest_cost run_cost m.total_cost pyo.Objective(ruleobj_rule, sensepyo.minimize)求解时指定求解器和参数比如Gurobi或CBCsolver pyo.SolverFactory(gurobi) solver.options[MIPGap] 0.01 # 1% 最优间隙 solver.options[TimeLimit] 3600 # 秒 results solver.solve(m, teeTrue)值得注意的是上面的氢平衡是简化版没有体现“先满足合成氨再决定充放氢”的逻辑。但因为有储氢罐容量约束和成本目标求解器会自动找到最优的充放策略模型逻辑其实是完备的。3.4 结果后处理与可视化优化结束后先导出容量变量和逐时段调度变量到 DataFrame再画图。最常用的可视化是典型日的功率平衡堆叠面积图和储氢罐储量变化曲线。import matplotlib.pyplot as plt # 画某个典型日的功率平衡 s 0 plt.figure(figsize(10, 5)) plt.stackplot(range(24), [data.P_wt[i] for i in range(24)], [data.P_pv[i] for i in range(24)], [data.P_grid_buy[i] for i in range(24)], labels[WT, PV, Grid Buy]) plt.plot([data.P_el[i] for i in range(24)], labelElectrolyzer, colork) plt.legend() plt.xlabel(Hour) plt.ylabel(Power (kW)) plt.savefig(output/figures/dispatch_scenario0.png, dpi150)除了把图生成出来还要把每个时段的电力去向、氢去向计算成“能量流表格”方便你验证功率平衡约束是否真的被满足。这是排查模型错误的第一手段。4. 复现过程中的坑与解决实录4.1 典型日调度的“场景跳变”问题第一次复现时我图省事从全年数据聚出12个典型日后每个典型日独立求解调度再把结果乘以权重加总。结果算出最优装机容量时发现光伏容量被严重夸大甚至到了不可思议的程度。原因是典型日之间没有关联模型会利用每个典型日的最佳日照来设计容量但实际天气无法每天都是典型日这是一种“数据窥视”偏差。解决办法是让全年或所有典型日进入同一个优化模型所有典型日共享同一套容量变量而调度变量各自独立。这样容量选择必须同时满足所有典型日的运行约束结果才合理。4.2 离网模式下“无解”的经典缘由离网模式跑起来最容易遇到Model is infeasible。一开始我还以为是线性化出错后来逐步加约束检查发现最核心的矛盾在于离网系统的电力完全来自风光而如果连续阴雨天外加静风系统既发不出电也氢库存用尽合成氨装置却要求连续运行模型自然无解。解决方案是在模型里允许一定程度的“失负荷”但不是粗暴达标。具体做法是加一个表示不足电量的变量 (P_{short}(t))目标函数里对缺电量和缺氢量设置很高的惩罚成本。这样模型宁愿选择一个稍微大一点的储氢罐或放弃极端天气下的生产也不至于无解。工程含义也很直接离网系统难以保证100%可靠性通常允许一定概率的供电缺额。4.3 电解槽效率曲线的线性化精度电解槽的电耗和产氢量并不是简单的正比关系低负荷下单位产氢能耗反而更高。论文里常用多项式拟合功率-产氢曲线。我最初用两点线性拟合即恒定效率结果模型给出的最优电解槽容量明显偏小因为系统忽略了低负荷效率下降。改用分段线性化后问题解决了。参考经验把负荷区间从10%到100%分成5到6段每段用线性关系近似分段点设在20%、40%、60%、80%和100%。分段数继续增加对结果影响不超过1%但求解时间却可能慢好几倍。我日常做敏感性分析时用3段最终方案才用6段。4.4 二进制变量过多导致求解器卡死容量变量虽然是连续变量但设备选型通常需要整数比如风机台数。加入整数台数后模型二进制变量会增多求解速度变慢。尤其当容量上界设得太大时变量搜索空间指数增长Gurobi可能在数分钟内都得不到一个可行的整数解。我的做法是先用连续变量求解一次得到近似容量值然后在该值附近设置整数变量的上下界缩小搜索空间。比如风电连续最优解是41.7 MW那就限制风机整数台数范围为[30台, 55台]这样不会影响全局最优又能大幅加速收敛。4.5 结果合理性检验模型输出后不要急着写报告先用物理量纲和能量守恒去验算输入电量 电解槽用电 辅助用电 卖电 弃电量是否成立电解槽产氢总量 合成氨耗氢 储氢净变化是否守恒单位氨耗电总用电/合成氨产量是否在一个合理范围通常是50到60 kWh/kg NH3 左右最优容量再算一次全年调度看看储能初始/终值是否一致弃电率是否过高。这些检查非常有效我曾经发现某个模型结果里电解槽年利用小时数只有2000多小时明显偏离工程预期追查下来是电解槽最小负荷约束写反了。5. 体会与后续扩展5.1 复现之后的一些认知这个项目让我最深的理解是容量优化只是壳里面的调度约束才是灵魂。很多文献会强调算法多先进实际落地时真正影响结果的是电解槽的最小负荷率、储氢罐的初始容量约束、合成氨装置的连续运行要求以及天气数据本身的质量。如果你只是随手给几个常数很容易得到一个看起来“最优”但是在现实里完全无法运行的系统。还有一个体会是Python 生态下 Pyomo Gurobi 已经足够应对这类中小规模的能源优化问题。真正需要花时间的地方是数据前处理和结果后处理而不是模型本身。把数据、场景、约束写清楚任何一个免费求解器也能跑出让人信服的结果。5.2 可以继续扩展的方向如果想把这套模型再往深做我个人比较推荐几个方向一是多目标优化。除了经济成本还可以把单位氨碳排放、弃电率、系统自供电率也放进目标函数用加权求和或帕累托前沿分析找平衡解。二是不确定性优化。风电和光伏都是随机变量可以把容量规划问题扩展成两阶段随机规划或鲁棒优化让容量决策能够应对天气场景的波动而不仅仅是几个典型日的确定性场景。三是与模型预测控制MPC结合。本文模型解决的是长期容量和年度运行策略但实际生产中还需要小时级或分钟级的实时调度。把优化结果送给MPC做滚动修正是一个非常自然且应用价值很高的扩展。四是加入电解槽退化模型。电解槽在高波动功率下运行寿命会加速衰减。把退化成本嵌入目标函数后容量和调度结果都会发生变化尤其是会倾向配置一部分电池储能来平滑功率这在很多研究报告里已经被验证。如果你正在复现类似项目建议先搭一个1个典型日的小规模原型把约束和目标调通再逐步扩展到全年场景。不要一上来就追求完整复现那样调试一个不可行问题时根本分不清是模型错误还是数据问题。这个项目最大的收获不是那几个容量数字而是整个“建模—求解—验证—再修正”的闭环思维希望这份记录也能让你少走几步弯路。