我复现这个“并/离网风光互补制氢合成氨系统容量-调度优化”项目时最大的感触是把CPLEX调通只是刚开始真正花时间的地方在于把风力发电、光伏出力、电解槽制氢、氢储罐缓冲、氨合成装置连续运行这些物理环节统一翻译成一个能让CPLEX高效求解的混合整数线性模型。标题里的“复现”两个字恰恰意味着论文里不会给你写好的Matlab代码你手里只有几张资源曲线、一张设备参数表、十几个约束公式。我在这篇里会把我复现时的建模思路、Matlab调用CPLEX的代码骨架、以及来回踩过的坑都整理出来目标是你拿到自己的风资源数据和光伏数据之后能照着搭出同样一个可运行、可解释的双层容量-调度优化框架。这个项目适合谁读我的判断是三类人最适合第一类是正在复现可再生能源制氢类论文的研究生尤其是遇到“容量和调度耦合”“并网离网对比”这类问题的新手第二类是从事氢能或风电光伏项目前期规划、想做经济性测算的工程师第三类是已经会用Matlab但还没有把CPLEX接到自己的优化模型里、希望看一个完整调用链路的技术人员。我不打算讲太多可再生能源的常识而是直接聚焦到“这个优化是怎么建模的、怎么用CPLEX解的、哪里容易出错”这三件事上。1. 项目整体从风光到氨的物理链条与它背后的优化问题1.1 系统的四个关键环节风、光、氢、氨先把这个系统的物理链条捋清楚。风电场和光伏电站发电之后电能优先送入电解槽电解水得到氢气。氢气可以先存进储氢罐这是一个缓冲环节因为风电和光伏出力在小时级尺度上有很大的随机波动而下游氨合成装置恰恰希望原材料供给尽量稳定。氨合成装置把电解产生的氢气与空分装置提供的氮气在高温高压下合成氨产物是液氨或者气氨。以完整的制氢合成氨流程作为优化背景涉及到的设备一般至少包括风电机组、光伏阵列、电解槽碱性或PEM、氢气储罐、空分装置、氨合成塔。部分论文中还会加入蓄电池或压缩机组但我们就先用最精简的流程来理解后面扩展也容易。这套系统的核心特征在于“跨能流”和“跨时间尺度”。电能是瞬时的发出来那一刻就要被消耗或者被丢弃所以在离网场景下每一小时都必须满足电力平衡。氢气是有存储缓冲的白天光伏出力大时多制氢存起来夜里用储存的氢继续维持氨合成。氨生产也是连续性的从系统运行经济性来讲不希望氨合成装置频繁启停因为频繁启停不仅损失催化剂寿命还会带来额外的能耗。这些特性最终都会变成优化模型里的约束条件。1.2 “容量”和“调度”到底分别优化什么标题里提到的“容量-调度优化”很容易被理解成“先定容量再做调度”但我在实际复现中发现这篇论文的建模重点其实是两者的联立。容量层解决的是规划问题风电场要装多大、光伏要装多大、电解槽额定功率是多少、储氢罐容量多大、氨合成装置规模多大。这些变量的量纲是MW、MWh、kg或t/d。容量变量决定整个项目的固定投资成本所以目标函数里资本成本的量级比较大。调度层解决的是运行问题给定一组容量之后在一天24小时或者一年8760个小时的时间序列上每小时我应该让风电多发还是弃掉一部分、从电网买电还是自己发电、电解槽在哪个功率水平运行、储氢罐是充放还是保持不动、氨合成装置产率是否调节。调度变量通常都是小时级甚至15分钟级的时间序列变量目标里对应的是运行成本和收益。这两种决策天然耦合在一起因为调度变量的可行域会受到容量变量的限制。比如电解槽额定功率越大调度时能够消纳的风电就越多弃风率可能降低但它每小时运行时的单位成本、以及购买电解槽的成本也在提高。如果在仿真的过程中只做一个单层优化会得到一个过于乐观或者过于保守的结果。正确做法是通过一个统一的目标函数来权衡容量投资和长时间运行的效果同时优化两类变量。这个技术细节是整个复现项目的关键点。1.3 并网与离网两种场景对模型约束的影响论文里以“并网”和“离网”两种系统的对比来开展我的理解是这两个场景下模型结构与难度区别很大。并网系统允许与公共电网交换功率。优化时需要考虑分时电价比如夜间电价低就可以优先从电网购电来制氢既降低了制氢成本又避免了白白扩大风电场容量。但并网通常不允许同一根母线上同时买电和卖电这意味着需要引入整数变量区分购电/售电状态模型从线性规划变成混合整数线性规划求解复杂度和前一种情况完全不同。离网系统完全靠自产新能源供给没有大电网作为后盾。这时候电力平衡约束从“柔性可交易”变成“刚性自平衡”每一小时都必须有充足的风光出力和存储缓冲来满足电解槽、空分等负荷的用电。在设备容量和调度计划上它对储氢罐容量的依赖会显著增加。我复现下来还发现离网系统的可行域经常非常狭窄稍有不慎就会在某个时段出现不可行解。所以在构建模型时我倾向于把主模型搭成“可切换模式”的形式加一个模式开关通过修改边界条件来分别模拟并网和离网而不是维护两套代码。2. 容量-调度耦合模型的建模细节与关键参数2.1 目标函数用年化费用统一规划与运行很多刚接触这种优化模型的人都会问一个问题容量变量是长期的调度变量是逐时的两个量纲不一样怎么放到一个目标函数里我的做法和大多数文献一致采用“年化成本”的方式。先把所有一次性投资按寿命年限做等年值化再与一年的运行费用相加。等年值化公式很简单$$ A C_{\text{总}} \cdot \frac{r(1r)^n}{(1r)^n - 1} $$这里的(C_{\text{总}})是初始单位投资(r)是贴现率(n)是设备寿命。乘出来的结果叫资本回收系数它把一次性投资拆到每一年。% 年化系数计算示例 r 0.08; % 年贴现率 n_wind 20; % 风电场寿命 crf_wind r * (1r)^n_wind / ((1r)^n_wind - 1); annual_wind_inv crf_wind * unit_cost_wind * W_cap; % W_cap单位kW这样目标函数可以写成最小化 sum(设备年化投资费用) sum(年运行维护费用) sum(每年外购电费用) sum(每年弃电惩罚如果有的话) - sum(售电或售氨收益如果有的话)在并网场景里外购电费是根据分时电价计算的每小时的购电功率乘电价再累加。离网场景中没有外购电这一项同时为了确保供电可靠性一般会在目标里加入“负荷失供惩罚”或者“弃电惩罚”这样模型就不会随便跳出一些极端但不符合工程直觉的解。2.2 电力平衡约束从“可以买电”到“必须自平衡”对于并网场景调度层每小时要满足P_wind(t) P_pv(t) P_buy(t) P_el(t) P_airsep(t) P_sell(t) P_curtail(t)离网场景则是这个等式的简化版把(P_{\text{buy}})和(P_{\text{sell}})去掉变成P_wind(t) P_pv(t) P_el(t) P_airsep(t) P_curtail(t)如果你希望通过并网时的购售电价差来决定购电还是售电还需要给(P_{\text{buy}})和(P_{\text{sell}})之间加一个互斥条件。常用的方式是引入二元变量(y_{\text{buy}}(t))P_buy(t) M * y_buy(t) P_sell(t) M * (1 - y_buy(t))这里(M)是一个足够大的常数。不过提醒一句(M)别取得过大否则CPLEX的数学规划松弛质量会变差求解时间大幅上升。我通常会让M取该线路或设备最大允许功率的1.1倍这样既保证能放开约束又不至于破坏分支定界的下界。2.3 电解槽与制氢环节的建模电解槽是最核心的可调负荷。它有一个很实用的特性在相当宽负荷区间内电耗量和产氢量近似可以按线性关系近似处理。虽然严格意义上电解槽效率会随电流密度变化但许多容量优化论文里会用一个等效电耗率比如(5.0\sim5.6\text{kWh/Nm}^3)来建模。为了统一量纲我建议把产氢量用质量单位“kg”而不是体积单位“Nm³”放在模型里。因为储氢罐的约束、氨合成装置的氢耗约束都用质量表示更方便。折算大致为1 kg氢气约等于11.2 Nm³标准状态下所以如果设备电耗是5.6 kWh/Nm³那么相当于约62.7 kWh/kg。电解槽的建模约束m_H2_prod(t) P_el(t) / W_el_per_kg % 每小时产氢量kg 0 P_el(t) P_el_cap如果原论文要求电解槽不能长期在低负荷运行还要加一条最小运行负荷约束P_el(t) min_load * P_el_cap * y_el_on(t) P_el(t) P_el_cap * y_el_on(t)这个约束是典型的“二进制变量容量变量”耦合例子。当把容量(P_{\text{el_cap}})也作为优化变量时模型就出现了容量变量乘以二元变量的双线性项。我在复现类似模型时确实被这个卡过很久。解决办法有两个方向一个是用大M法把双线性项线性化另一个更常用的办法是先把容量决策从小时变量中分离出来在多层迭代框架里处理。我后面会专门讲这两种建模策略的取舍。2.4 储氢罐的动态方程与周期接续约束储氢罐的作用相当于电路中的电容它缓冲氢气的产消不平衡。它的动态约束是一个离散时间状态方程S_H2(t 1) S_H2(t) m_H2_prod(t) - m_H2_to_ammonia(t) - m_H2_vent(t) S_min S_H2(t) S_max其中(S_{\min})一般取最大容量的一定比例比如5%作为安全余量(S_{\max})就是储氢罐容量变量。很多第一次复现这类模型的朋友会漏掉周期接续约束也就是一个调度时间域结束后的储氢量要等于初始储氢量或者至少要满足结束时的水平不低于一个给定值。如果不约束这个CPLEX会钻空子把最后一个时段的储氢全部用光或者全部排空造成结果失真。我在模型中一般加这样的约束S_H2(T 1) S_H2(1) S_H2(T 1) S_H2(1)如果使用的是典型日而不是连续一年的数据则更常见的是让每个典型日首末储氢水平相同。不过这样的话必须注意典型日之间的跨日传递问题有时候需要用“多个典型日串联共享同一个初始状态”的做法来避免过度乐观的逐日独立优化结果。2.5 氨合成环节化学计量比与最小连续运行约束合成氨反应的化学计量式是[ N_2 3H_2 \rightleftharpoons 2NH_3 ]按质量计算每生产1吨氨大约需要0.176吨氢气。这个比例几乎在任何复用Haber-Bosch工艺的系统里都不能随便偏离因为它是化学反应的严格计量约束。调度层里需要约束m_H2_to_ammonia(t) lambda_NH3 * m_NH3(t) 0 m_NH3(t) R_NH3_cap如果追求更真实的运行约束还应对氨合成装置设置最小生产负荷和爬坡速率限制m_NH3(t) L_min * R_NH3_cap * y_NH3_on(t) m_NH3(t) R_NH3_cap * y_NH3_on(t) |m_NH3(t) - m_NH3(t-1)| ramp_rate * R_NH3_cap空分装置的耗电也不容忽视。通常氮气制备会有少量电力消耗估算时可简化为每吨氨耗电0.20.3 MWh具体与空分规模有关。我没有把它当成主体负荷但它在离网场景中可能占电解槽耗电的5%左右在精确的电力平衡计算里不能随手省略。2.6 可再生能源出力数据从气象时序到模型系数风光资源数据是模型的外部输入通常给的是逐小时的风速和辐照度、温度数据。在实际编码中我会先把它们转成“单位容量出力序列”也就是假设风电场额定容量是1 MW时每个小时能发多少电。这样容量变量只是把这段单位时序乘以实际容量不会引入非线性% 单位容量出力序列: pu_w(t), pu_pv(t) % 实际出力 W_cap * pu_w(t) 和 PV_cap * pu_pv(t)这一步骤非常重要。因为很多新手直接把风速曲线乘一个功率曲线常数后作为固定发电量代入那就无法把容量作为决策变量了。我一再强调资源数据必须处理成“每单位装机容量”的出力序列才能让容量变量成为模型的直接决策变量。3. 用Matlab调用CPLEX的代码实现细节3.1 为什么推荐YALMIP来搭桥Matlab本身没有内置CPLEX接口需要单独配置有几种方法。第一种是直接用IBM提供的Cplex类比如Cplex(filename, model.lp)或cplexmilp函数第二种是使用第三方建模工具箱YALMIP先写好符号优化变量和约束再在求解时指定solvercplex由YALMIP自动转换成CPLEX的输入格式。第二种方法对于科研复现来说更友好因为模型可读性强后期改约束很快不容易弄错维度。我实际使用的是YALMIP加CPLEX的组合。YALMIP可以用sdpvar定义连续变量用binvar或integer定义整数变量写约束和写普通数学公式一样阅读起来非常直观。% 决策变量容量变量 W_cap sdpvar(1,1); PV_cap sdpvar(1,1); El_cap sdpvar(1,1); H2_cap sdpvar(1,1); NH3_cap sdpvar(1,1); % 调度变量时间序列 P_buy sdpvar(T,1); P_sell sdpvar(T,1); P_el sdpvar(T,1); S_H2 sdpvar(T1,1);这里sdpvar的维度和规模都直接对应小时数T。如果T是8760小时变量数量会很大如果采用典型日比如12个典型日、每个典型日24小时那么T实际是288规模就很容易控制在CPLEX能秒解的范围内。3.2 典型日选取与权重设计在复现这类容量优化时全年8760小时数据直接代入会让MILP规模超出合理范围尤其当还有大量二元变量时。所以我通常采用两部走先用kmeans算法从全年气象和历史负荷序列中聚类出若干个典型日再把每个典型日的目标函数乘以该类别在全年的天数权重。% 假设 RawData 是 N*48矩阵N为全年天数前24列为风速单位出力后24列为光伏单位出力 [idx, C, sumd, D] kmeans(RawData, K); % K12 或16 weights histcounts(idx, 1:K1);聚类后模型中的目标函数就变成sum(weights(k) * sum_24(运行成本(k)))同时所有时间序列变量都只在一个典型日内部流动典型日之间通过权重在目标函数中汇总。这个处理在容量规划里已经被广泛接受因为设备年利用小时数的整体趋势是靠气象聚类覆盖的而不是用某一个代表日硬扛全年。3.3 单层MILP实现容量-调度联立优化的思路标题里的“容量-调度优化”有很多论文会选择把容量作为决策变量和调度变量放在同一个混合整数规划中一次性求解。这个思路的可实现性很高只要满足两个条件第一所有约束关于容量变量都是线性的第二目标函数中不存在容量变量乘以调度变量的非线性项。我特别想说清楚一个坑。比如电解槽约束要求“如果运行功率要大于最小技术出力”写作P_el(t) min_ratio * El_cap * y_el_on(t)这个约束看起来很简单但实际上同时包含连续变量(El_cap)和二元变量(y_el_on(t))二者相乘是一个双线性项不是标准LP/MILP约束。如果直接扔给YALMIP它会调内置非线性求解器或者报错而不会走CPLEX的MILP分支。解决办法是引入辅助连续变量(z(t))代表“电解槽在t时刻使用的额定容量比例”z(t) y_el_on(t) * El_cap然后再通过大M法线性化。但这样的线性化会增加变量和约束数量而且如果容量变量已经有多个设备、多个时段会非常麻烦。所以我在复现中更常见的做法是将这类约束改写为P_el(t) min_ratio * El_cap - M * (1 - y_el_on(t))这实际上是同一个约束的“松约束”形式含义是当(y_{on}1)时(P_{el}(t))必须高于最小负荷当(y_{on}0)时右边减去一个很大的M约束就不起作用。这样写之后CPLEX可以当作标准MILP处理。3.4 约束编写与优化求解的完整骨架为了让你能直接上手我把一个简化版的Matlab/YALMIP骨架整理出来这个骨架包含并网场景的大多数核心约束。实际使用时你可能还需要扩展设备数量、典型日权重、收益项等。%% 参数 T 24; % 每小时调度一个典型日24点 pu_w 随机或实际生成...; % 1×24, 单位容量出力 pu_pv 随机或实际生成...; %% 容量变量 W_cap sdpvar(1,1); PV_cap sdpvar(1,1); El_cap sdpvar(1,1); H2_cap sdpvar(1,1); NH3_cap sdpvar(1,1); %% 调度变量 P_el sdpvar(T,1); P_wind sdpvar(T,1); P_pv sdpvar(T,1); P_buy sdpvar(T,1); P_sell sdpvar(T,1); P_cur sdpvar(T,1); % 弃电功率 S_h2 sdpvar(T1,1); H2_prod sdpvar(T,1); % 每小时产氢 kg H2_to_nh3 sdpvar(T,1); m_nh3 sdpvar(T,1); % 每小时产氨 kg y_switch binvar(T,1); % 1 表示购电0表示售电 %% 约束 C []; % 风光出力上限 C [C, P_wind W_cap .* pu_w]; % 实际风容量*单位出力 C [C, P_pv PV_cap .* pu_pv]; C [C, P_wind P_pv P_buy P_el P_airsep_const P_sell P_cur]; % 购售互斥 M 5000; % 根据系统最大功率取 C [C, P_buy M .* y_switch]; C [C, P_sell M .* (1 - y_switch)]; % 电解槽 C [C, H2_prod P_el / wh_per_kg]; % wh_per_kg约62 C [C, P_el El_cap]; % 储氢 C [C, S_h2(2:T1) S_h2(1:T) H2_prod(1:T) - H2_to_nh3(1:T)]; C [C, S_h2 0.05 * H2_cap]; C [C, S_h2 H2_cap]; C [C, S_h2(1) S_h2(T1)]; % 氨合成 C [C, H2_to_nh3 0.176 * m_nh3]; C [C, m_nh3 NH3_cap]; %% 目标年化投资运行成本 obj crf_w * cost_w * W_cap crf_pv * cost_pv * PV_cap ... crf_el * cost_el * El_cap crf_h2 * cost_h2 * H2_cap ... crf_nh3 * cost_nh3 * NH3_cap; obj obj 365 * sum(P_buy .* elec_price); obj obj 365 * 1e-3 * sum(P_cur); % 弃电惩罚 %% 求解 ops sdpsettings(solver, cplex, verbose, 2); ops.cplex.mip.tolerances.mipgap 1e-4; ops.cplex.timelimit 600; sol optimize(C, obj, ops);这里的P_airsep_const是我为了简化处理设置的空分固定负荷真实建模时也可以把空分功率和氨产量挂钩变成变量。由于是示意代码很多单位换算和边界条件还需要按你实际数据修正。3.5 CPLEX求解器参数的使用心得CPLEX作为商业求解器它的核心优势在于大规模MILP的分支定界速度和稳定的数值处理能力。实际使用中我一般最关心mipgap、timelimit、nodefile和threads这几个参数。ops.cplex.mip.tolerances.mipgap 5e-3; % 0.5%的gap足够工程使用 ops.cplex.threads 8; % 并行线程数 ops.cplex.mip.strategy.file 2; % 节点文件压缩策略内存紧张时开如果发现在一个小时内没有达到理想gap不用太焦虑。先看CPLEX日志里的MIP gap是百分之几。如果gap已经低于1%但模型规模导致时间很长此时强行追求gap到0.01%意义并不大因为容量规划的最优解在1%的gap范围内已经完全能支撑工程决策。如果模型报“infeasible”而你又很难排查是哪条约束导致不可行我个人建议在目标函数中引入一个松弛变量比如失电负荷或者先固定容量变量成一组合理值再做纯调度优化这样可以把可行性问题定位到是容量不足还是调度约束有矛盾。4. 复现并网/离网系统时的实操经验与常见问题4.1 结果怎么看容量结果与调度结果互相校验跑出结果后不要急着截图写报告。至少要做两层校验。第一层校验是总量平衡。以100 MW风电场、50 MW电解槽为例假设风电场年利用小时数为2200小时年产电量就是2.2亿kWh。如果电解槽制氢单位电耗为5 kWh/Nm³即约56 kWh/kgH2那么一年最多可产氢约3928吨。如果按每吨氨消耗0.176吨氢折算这些氢大约对应年产氨2.23万吨。如果你的优化结果给出年产氨3万吨那就说明风电利用小时数或者电解槽容量设置出了偏差需要回头检查输入条件。第二层校验是趋势合理性。离网场景中储氢罐的储量曲线应该呈现昼夜波动并在连续运行中逐步回补不可能出现连续几天放空后又在几个小时内充满的情况。并网场景的购电行为应与电价波动相对应夜间电价低谷时应该有更多购电制氢白天高电价时段则会减少购电甚至售电。如果完全反着来大概率是目标函数中价格的符号写反了或者单位不匹配。下面这张表是我复现时常用的一个快速自检表能帮你快速定位问题方向。检查项常见异常可能原因购电量与电价序列高电价时段反而大量购电目标函数价格符号错误或功率单位弄错储氢罐水平长期处于下边界或末时刻被清空缺少周期接续约束或储氢容量边界取错弃风/弃光量离网系统弃电量为0且无失负荷可能遗漏了切负荷变量系统被过度约束产氨量高于化学计量上限氢耗系数取错或单位没统一CPLEX求解时间几小时不收敛整数变量过多或M值太大不可行问题solver直接返回infeasible容量变量太小或缺少松弛变量4.2 离网系统里的“假可行”问题离网系统的运行目标强调了供电可靠性。如果在目标函数里既没有失电惩罚项又没有设置可靠性约束优化器很可能会选择极小的容量配置让系统在全年绝大部分时段都无法供电但目标函数还是最优。这显然不符合期望。所以在离网模型中我建议用一个“失电负荷”或“氢产量不足”的连续变量来兜底对它的单位惩罚量定到非常高的水平让优化器只在极个别极端的时段选择牺牲可靠性。对于容量规划还可以进一步限制系统的“自供电率”或“氢供应不足率”不超过某一阈值例如利用年氢缺量比例不超过2%写成累积约束sum(m_H2_shortage(t)) 0.02 * sum(m_H2_demand(t))这类约束能有效抑制“过度牺牲可靠性换经济性”的倾向也更符合实际情况。4.3 并网波动的整数变量爆炸与应对并网模式下如果每小时的购售状态都用二元变量表示全年的二元变量数就是8760个。虽然对CPLEX而言不算极端但一旦加上每个电解槽的启停变量、氨合成装置的最小运行状态变量变量规模会迅速膨胀。再加上典型日方法会引入周期性耦合模型一般会在数分钟内出现解然而如果想分出较优的MIP解常需要较长时间。为了控制模型复杂度我的实践经验有三个方向一是把时间分辨率从1小时拉长到2小时或3小时在某些阶段先粗仿一遍找趋势再选取关键区间加密分辨率二是合并同型号电解槽用聚合机组模型替代逐台模型这样启停变量从N台减少为机群级别三是用滚动时域的思路先固定容量结果再做更细时间粒度的调度校验不必把全年8760个小时都塞进同一个MILP。4.4 关于“复现”这件事代码不是唯一重点最后回到“复现”这个词。我在复现多篇同类型论文后有个感受代码层面的复现和论文结论层面的复现并不是一回事。跑通代码只能说明你把公式写对了但论文里的结果往往基于特定的风资源数据、光伏数据、设备价格和贴现率而论文并不会把所有数据都放到附录里。当你换了一套数据以后最优容量比例、氨生产成本很可能发生较大变化。因此复现的价值不在于得到一个和论文一模一样的电价曲线或容量结果而在于你能解释为什么论文在某个电价水平下选择“多装风电解槽少装氢罐”这才是容量-调度优化问题的核心规律。实际项目中让我印象最深的一次调整是同一套风力数据和光伏数据只把贴现率从10%降到6%最优的光伏容量占比就明显上升因为高贴现率会惩罚初期投资大的设备而光伏单位投资虽然低但利用小时数低在高贴现率下并不占优。这也再次说明容量与调度联立优化的必要性如果不做全日历逐时模拟单靠经验公式很可能会得到错误的装机比例。从工程落地的角度看这类优化结果通常还要再做一次极端气象年的稳健性校验。资源曲线如果换成枯风年或少光年系统是否仍然可靠、成本是否会大幅上升这些都是容量规划中必须回答的问题。而这些扩展大多不需要改变模型结构只要把输入数据替换后再求解一遍即可这也是用YALMIP加CPLEX这种可读性较强的方式建模的额外好处。如果你正在复现类似的“风光互补制氢合成氨”容量-调度优化我个人建议先从离网场景入手因为它约束清晰、可行域窄方便你快速发现逻辑错误等离网模型稳定性、自平衡和储氢缓冲的逻辑都通过验证后再切换到并网场景加入购售电变量。这样的递进式复现路线会比一上来就同时处理两种模式省去大量检查时间。