这个题目我第一次看到是在一个新能源项目的交流群里有人发了一段代码截图配了一句话“风光制氢合成氨容量和调度一起优化Cplex跑起来特别酸爽。”我当时第一反应是这标题也太长了但仔细一看这其实是一个特别典型的“新能源电力系统耦合化工生产”的优化问题而且用Matlab配合Cplex来做容量规划与日前调度的联合优化确实是目前学术论文里最常见的复现方向。我后来花了一周时间把这个模型从零搭起来把代码跑通也踩了不少坑。这篇就把我整个复现过程中最核心的建模思路、求解器配置、代码结构、调试经验和结果分析套路都写出来给想复现同类工作的朋友做个参考。1. 项目到底在做什么风光制氢合成氨的容量优化与调度优化很多人一看到“容量-调度优化”就发懵觉得这是两个问题。实际上容量优化是决定“系统里每种设备装多大”调度优化是决定“在每一时刻设备怎么运行”。这两个问题相互耦合不能分开做。比如电解槽如果装大了风光充裕时能多制氢但投资成本高闲置时就是浪费储能罐装大了能平滑氢气的供需波动但同样带来成本压力。而调度策略又决定了在给定容量下系统能不能把风光的弃电率压下来、能不能满足合成氨的连续供氢需求。所以必须把容量变量和调度变量放在同一个优化模型里同时算。1.1 风光互补制氢合成氨的完整链条先把这个系统的物理链条捋清楚不然看代码容易看晕。整个系统的能量流和物质流大概是这样的风力发电机和光伏板发出直流/交流电经过整流后供给电解槽。电解槽消耗水和电产出氢气和氧气氢气经过纯化后压入储氢罐。储氢罐缓冲氢气供应避免风光波动直接冲击合成氨工段。合成氨装置需要氢气和氮气氮气通常来自空分装置在高温高压下合成氨。氨可以储存在液氨储罐中作为最终产品输出。在这个链条里电能、氢能、热能合成氨反应放热、电解槽需要热能但通常是电加热都耦合在一起。做优化的时候核心是电力平衡和氢气平衡两条主线其他如水的消耗、空分能耗可以作为常数处理或简化约束。并网和离网的区别在于系统是否与外部电网存在功率交换。并网模式下系统可以在风光出力不足时从电网购电补充电解槽和合成氨的用电也可以把多余的电卖给电网离网模式下系统完全不依赖电网只能靠风光发电、储能如果有电储能和负荷侧的调节来维持平衡。离网模式对容量配置的要求更高因为必须保证任何时刻都不能缺电缺氢所以往往需要更大的储能容量和更保守的调度策略。1.2 并网与离网模式下的优化逻辑差异先说并网模式。并网时电网相当于一个无限大的“缓冲池”缺电就买富余就卖。这个模式下的优化目标通常是整个系统在全生命周期内的总成本最小包括设备投资成本、运维成本、购电费用减去售电收益和氨产品销售收入。容量优化的重点在于风光装机、电解槽功率、储氢容量之间存在经济平衡点。比如如果当地上网电价高、购电电价低那么系统会倾向于少装电解槽多用电网电来制氢如果风光资源特别好可能更多装电解槽去吃弃电。离网模式则完全不同。离网系统没有电网兜底本质上是个孤岛。这时候“可靠性”变成硬约束任何时刻电力都要平衡氢气也不能断供。为了满足极端天气连续阴天、无风期下的负荷需求要么把风光装机做得很大要么配置更多的储能还要在调度策略上允许适当削减负荷或停产。容量优化结果往往比并网模式“看起来更浪费”因为多出来的容量是为了极端工况准备的虽然利用率低但能保证系统不停摆。1.3 为什么容量和调度要联合优化而不是分步做有些论文为了简化会先优化容量再在给定容量下优化调度。这种“两步法”的问题是容量方案和调度策略是互相影响的。如果你先用一个粗糙的调度规则去优化容量得到的容量配置可能在一个更聪明的调度策略下并不是最优的。比如电解槽如果在低价电时段集中满负荷运行、在缺电时段完全停机那么需要的储氢容量就比较小如果要求电解槽连续平稳运行那储氢罐就要大得多。反过来调度策略又受容量限制——储氢罐就那么大的容积你不可能让它在某个时段超额存储。所以联合优化的本质是在同一套数学模型中把设备规模通常是整数或连续变量和逐时段运行功率连续变量一起当作决策变量来求解这样才能得到真正全局最优的设计方案。2. 系统建模先把物理过程翻译成数学约束建模是整个复现工作的核心也是最容易出现“改了一晚上约束还是不可行”的地方。建模的思路其实不复杂把所有设备都抽象成“输入-输出”的关系用变量表示每个时段的状态把物理规律写成等式或不等式约束最后把所有成本或收益加起来写成一个目标函数。2.1 风光出力数据与时间序列处理风光出力是整个模型的外部输入。通常的做法是用典型日数据或全年8760小时数据。复现的时候如果手头没有实测数据一般会用NASA或气象数据库的典型年数据或者直接生成几组典型日场景比如春夏秋冬四季、晴阴雨天、大风无风天。我复现时用的是论文里常见的“典型日法”挑出春、夏、秋、冬四个典型日每个典型日按24小时离散成24个时段。风电出力用威布尔分布拟合风速然后转换成功率光伏出力用辐照度转换。注意这里有一个很关键的细节风光出力的数值要归一化后再乘以装机容量这样在优化容量变量时比较方便。举个例子如果风电归一化出力曲线在某时段是0.45风光装机容量变量是120 MW那么该时段风电实际出力就是54 MW。这样容量变量和出力曲线可以解耦建模更灵活。2.2 电解槽、储氢罐与合成氨单元的建模电解槽的建模一般用效率法。电解槽消耗电功率 (P_{el}(t))产氢速率 (Q_{el}(t)) 可以用公式 (Q_{el}(t) \eta_{el} \cdot P_{el}(t)) 来表示其中 (\eta_{el}) 是电解槽的电氢转换效率单位是 (Nm^3/(kWh)) 或者 (kg/(MWh))取决于你要以质量还是体积计。注意电解槽有最小负载率限制比如质子交换膜电解槽的负载范围通常是30%~100%碱性电解槽可能是20%~110%甚至短时过载。这个约束必须写成 (P_{el,min} \cdot u_{el}(t) \leq P_{el}(t) \leq P_{el,max} \cdot u_{el}(t))其中 (u_{el}(t)) 是二进制启停变量否则调度结果会出现电解槽在极低功率下运行的工程上不合理的状态。储氢罐的建模就是动态平衡方程(S_{st}(t1) S_{st}(t) Q_{in}(t) - Q_{out}(t) - Q_{loss}(t))其中 (Q_{loss}) 可以是泄漏损失或者为了简化设置为零。储氢罐容量变量 (V_{st}) 会进入容量优化的决策变量集合而每个时段的存量 (S_{st}(t)) 必须小于等于 (V_{st})。这里要注意如果你把 (V_{st}) 作为连续变量那么储氢成本就是线性的如果储氢罐规格是离散的比如只能选50、100、200 Nm³就得引入整数变量。合成氨单元相对简单一点因为它通常要求平稳运行。合成氨反应是连续过程氢氮比固定3:1反应条件是高温高压所以短期内不宜频繁调节负荷。在做日前调度时简单模型把合成氨负荷设为一个固定值 (Q_{NH3,const})或者允许在额定负荷的70%到100%之间调节但必须加上“爬坡率约束”比如相邻时段氢消耗量的变化不能超过一定比例。否则优化器会让合成氨负荷剧烈波动来适应风光变化这在工程上是不现实的。2.3 并离网切换与能量平衡约束并网模式和离网模式的主要区别体现在电功率平衡方程中。并网模式[ P_{wind}(t) P_{pv}(t) P_{grid,buy}(t) P_{el}(t) P_{NH3}(t) P_{grid,sell}(t) ]其中 (P_{grid,buy}(t)) 和 (P_{grid,sell}(t)) 不能同时为正需要引入二进制变量避免同时购售电的荒谬结果。离网模式直接去掉这两项[ P_{wind}(t) P_{pv}(t) P_{el}(t) P_{NH3}(t) ]如果还有电储能蓄电池则电能平衡还要加入充放电项。不过在制氢合成氨系统里电储能经常被氢储能替代因为氢本身就是能量载体储氢罐既缓冲氢气又相当于能量存储。离网系统如果发现功率不平衡要么弃风弃光要么切负荷这两者在优化目标里都会有惩罚项。氢气平衡约束则是所有模式下都必须满足的[ Q_{el}(t) Q_{st,out}(t) Q_{NH3,in}(t) Q_{st,in}(t) ]这里每个变量都是非负的储氢罐不能同时进氢和出氢同样需要二进制变量或者通过线性化的方式约束。很多初学者在这里会漏掉“同一时刻不能同时充放”的约束导致优化结果出现储氢罐一边充一边放的自损行为。3. 优化模型与Cplex求解从非线性到MILP的关键转换很多人在这一步卡住因为一开始搭出来的模型可能是非线性的Cplex根本解不了。Cplex是一个求解线性规划、混合整数线性规划和二次规划的商业求解器它不认识非线性函数。所以建模的重头戏就是“线性化”。3.1 目标函数怎么设计目标函数通常分成两部分年度化投资成本设备容量乘以单位投资成本再用资本回收因子折算成每年成本。比如风力发电机单位投资成本是8000元/kW寿命20年折现率8%那么资本回收因子 (CRF \frac{r(1r)^N}{(1r)^N - 1} \approx 0.10185)即每kW风机的年化投资成本约815元。运行成本与收益包括购电费、售电收益、制氢/制氨销售收入、弃电惩罚等。有些论文还会把碳排放作为一个约束或目标项这里先不提复杂化以后再说。总之目标函数是线性项的组合这样才能被Cplex高效求解。我在复现时用的目标函数是[ \min \quad C_{inv} \sum_{t} \left( c_{buy}P_{buy}(t) - c_{sell}P_{sell}(t) c_{curtail}P_{curtail}(t) \right) - R_{NH3} ]其中 (R_{NH3}) 是卖氨的收益按全年产量计算(C_{inv}) 是所有设备的年化投资成本。如果你只是求调度优化容量变量固定那投资成本是常数可以直接去掉只保留运行调度项。3.2 线性化处理二进制变量和大M法的使用模型里最典型的非线性来源有三个电功率平衡中的“购电、卖电不能同时发生”通过二进制变量 (u_{grid}(t)) 和大M法处理(P_{buy}(t) \leq M u_{grid}(t))(P_{sell}(t) \leq M (1-u_{grid}(t)))。电解槽的“开停机”状态和“最低负载率”通过二进制变量 (u_{el}(t)) 来处理本质上也是大M法。储氢罐的“不能同时充放”可以用两种方式一是用二进制变量强行分开二是通过状态方程自然避免。其实如果目标函数没有故意奖励同时充放由于能量守恒同时充放只会白白增加损耗优化器通常不会选所以有些模型可以省略这个约束但为了稳妥我建议还是加上。还有一个经常被忽略的非线性是“容量变量乘以负荷系数”。比如储氢罐容量 (V_{st}) 和储氢状态 (S_{st}(t)) 之间的约束 (S_{st}(t) \leq V_{st})这其实是线性约束因为两个都是变量但一个随时间变化一个是全局容量。只要不出现“容量 × 出力系数”这种形式就都还是线性的。Cplex可以直接处理二进制和整数变量所以最后得到的模型是MILP这是Cplex最拿手的问题类型。如果模型里有平方项比如风机出力与风速的二次关系那就该考虑用分段线性化近似或者直接查表用离散点表示。3.3 Cplex在Matlab中的配置与调用这一步是很多复现者的噩梦。Cplex是商业软件需要安装IBM ILOG CPLEX Optimization Studio并且版本最好和Matlab兼容。安装时需要勾选“MATLAB”接口支持安装完成后在Matlab里运行“addpath(‘C:...\cplex\matlab\x64_win64’)”把接口路径加入环境变量。配置好以后有两种调用方式。一种是用Cplex的Matlab工具箱函数比如Cplex、cplexmilp等另一种是先用Yalmip建模再把模型交给Cplex求解。我个人更推荐Yalmip因为Yalmip写约束非常直观特别是带二进制变量的MILP模型Yalmip的代码可读性远高于直接用Cplex API。比如%% 定义变量 P_el sdpvar(24,1); % 电解槽功率 u_el binvar(24,1); % 电解槽启停 S_st sdpvar(24,1); % 储氢罐储量 V_st sdpvar(1,1); % 储氢罐容量 P_buy sdpvar(24,1); P_sell sdpvar(24,1); u_grid binvar(24,1);然后写上约束最后的求解命令是optimize(Constraints, Objective, sdpsettings(solver,cplex))这个方案的好处是你不需要记忆Cplex原生API的那些复杂的矩阵拼装过程Yalmip会自动帮你把模型转成Cplex能接受的稀疏矩阵。当然如果模型规模特别大原生API可能效率更高但Yalmip对于复现论文级别的模型完全够用。3.4 为什么用Cplex而不是fmincon或遗传算法这是复现者经常被问到的问题。制氢合成氨优化模型带有大量二进制变量天然是个混合整数问题。fmincon只能处理连续非线性规划即使把二进制变量当成0到1的连续变量松弛得到的解往往不是整数可行的而且非线性求解器对初值敏感常常收敛到局部最优。遗传算法倒是可以处理整数变量但它是启发式算法不保证全局最优而且对于带复杂约束的问题遗传算法的可行域探索非常麻烦经常给出违反约束的坏解。Cplex作为商业求解器采用分支定界法和割平面法能在可接受的时间内找到全局最优解或证明最优解。对于24个时段、五六个设备变量的模型Cplex通常在几十秒内就能搞定。所以选择一个合适的求解器本身就是此类优化问题的正确打开方式。4. Matlab代码实现模块拆解与关键代码片段复现这类代码最重要的不是把模型一股脑塞进一个脚本里而是拆成几个清晰的功能模块。我复现时把代码分成数据准备、模型构建、求解配置、结果输出四个部分每个部分单独调试最后组合。4.1 数据准备与典型日场景生成数据准备阶段要生成用于约束构建的所有常数。这些常数包括24小时的风电归一化出力系数、光伏归一化出力系数、分时电价、设备单位投资成本、技术参数电解槽效率、储氢罐初始存量、合成氨额定氢需量等。我习惯用结构体保存所有参数% 技术参数 param.eta_el 0.65; % 电解槽效率kWh/kgH2的倒数 param.P_el_max 100; % 电解槽最大功率单位MW这里可调 param.P_el_min_ratio 0.2; % 最小运行比例 param.H2_demand 5.5; % 合成氨每时段氢耗kg/h param.S0 100; % 储氢罐初始存量kg这里特别提醒单位必须统一。我一开始就因为把电解槽产氢速率单位混用了“kg/h”和“Nm³/h”导致数值差了一个数量级结果模型怎么调都不对。建议在代码里注释清楚每个变量的单位别嫌麻烦。典型日的选取直接决定了结果质量。如果你只用夏季典型日可能过度乐观只用冬季典型日配置会过于保守。复现时可以用k-means聚类从全年8760小时数据中选出几个代表日也可以直接手动挑四组数据每组赋权重形成“多场景加权”的模型。权重就是该典型日在全年中的占比这样容量配置能兼顾不同季节。4.2 模型搭建Yalmip约束构建示例下面给出一段我实际用的核心约束片段简化版帮助你理解Yalmip的写法%% 电功率平衡并网模式 Constraints []; Constraints [Constraints, P_wind P_pv P_buy P_el P_NH3 P_sell]; %% 购电和卖电不能同时进行 Constraints [Constraints, P_buy M * u_grid]; Constraints [Constraints, P_sell M * (1 - u_grid)]; %% 电解槽功率上下限与启停 Constraints [Constraints, P_el P_el_min_ratio * P_el_max * u_el]; Constraints [Constraints, P_el P_el_max * u_el]; %% 储氢罐动态平衡 Constraints [Constraints, S_st(2:end) S_st(1:end-1) Q_el(1:end-1) - Q_NH3(1:end-1)]; Constraints [Constraints, S_st V_st]; Constraints [Constraints, S_st(1) param.S0];注意Q_el在代码里是用P_el乘以效率算出来的我图省事直接写了表达式实际代码中可以定义一个辅助变量或者直接用表达式代入。4.3 求解与结果后处理求解前一定要确认Yalmip的sdpsettings里指定了求解器Cplex并设置合适的参数比如MIP间隙容忍度options sdpsettings(solver,cplex, cplex.mip.tolerances.mipgap, 0.01);0.01的间隙意味着求解器找到的解与最优解的差距在1%以内就停止。对于工程问题这个精度足够而且能明显缩短求解时间。如果你的模型特别大、一时半会儿跑不完可以适当放宽到0.02或0.05。求解完后用value()取出所有变量的值然后画图。我一般会画个三图合一上面是电力平衡图风电、光伏、电解槽、电网购售电中间是氢气储量图下面是电解槽启停状态和功率图。这三张图基本能把调度策略的合理性展示清楚。5. 复现过程中踩过的坑错误诊断与性能调优这部分是我最想讲的因为代码跑不通或者结果不合理时新手往往会盲目调参数结果越调越乱。下面这几个坑是我实际遇到过的几乎每个都能让人卡上一天。5.1 模型不可行先查“同时充放”和“同时购售”模型不可行最常见的元凶就是缺了“同一设备不能同时双向运行”的约束。Yalmip报错Infeasible problem后可以用check(Constraints)检查每条约束的残差也可以把约束拆开逐条排查。我建议养成一个习惯先用一小段数据比如只取6个时段构建简化模型跑通后再扩展成24小时。6小时模型求解快、容易定位错误等约束逻辑没问题了再用24小时。另一种不可行原因是“储氢罐初始存量和末端存量”的约束过于生硬。有些模型要求一个调度周期结束时储氢罐回到初始状态这可能导致无解。你可以加上“末端存量不小于某下界”的软约束或者把初始存量也作为决策变量允许优化器自由选择周期性调度。5.2 求解时间过长砍时段、加约束、调参数24时段MILP通常几十秒内能解完。如果你遇到求解器像死机一样长时间跑不完先考虑模型里是否出现了不必要的大M值。M值如果取得太大比如1e7会导致CPLEX的数值稳定性变差分支定界过程中出现大量虚假可行分支。正确做法是M取实际物理量级的1.2倍左右。比如购电功率最大不会超过风光装机加负荷总和那么M取这个总和的1.5倍就足够了。另外为每个时段引入的二进制变量数量也会影响求解效率。如果你有几个设备都需要启停变量变量规模呈线性增长问题不大但如果你用二进制变量去表示“多档位运行”或“分段线性化”二进制变量数量会爆炸式增长。这时候可以考虑减少分段数或者用特殊顺序集SOS1/SOS2来替代二进制变量Cplex对SOS约束的求解效率往往更高。5.3 数值尺度问题你可能被单位坑了这是复现里最阴间的问题。电解槽功率动辄上百MW储氢罐容量可能是几千kg而目标函数里的购电成本可能只有每MWh几百元这些数字放在一起时Cplex内部计算的矩阵条件数会非常差导致求解器返回“精度不足”或者结果异常。解决办法是把所有变量的量纲统一到同一个基准。我实际使用的做法是以“kg”为氢气单位、以“MWh”为能量单位、以“万元”为货币单位并把所有功率变量的数值控制在0.1~1000的量级内。如果发现某个变量的数值超过1e4就考虑调整单位或除以一个缩放因子。Cplex官方文档也建议用户把数据归一化到差不多的数量级这样数值稳定性会好很多。5.4 别忘了Cplex的许可证配置很多人在配置环节就卡住了。Cplex免费版Community Edition最多只能求解1000个变量或1000个约束对于稍大模型是不够用的。如果是学生可以使用IBM Academic Initiative申请免费full license这需要注册一个学术邮箱。拿到许可证后在Matlab里运行cplex.setup()或者在Yalmip里指定路径确保Cplex真的被调用而不是后台用了默认的内点法求解器。一个很坑的点是如果你安装了多个版本的CplexMatlab可能会误调用旧版本导致接口报错。建议在系统环境变量里明确设置路径并在每次运行前打印出cplex.getVersion()确认版本。6. 结果分析的维度和扩展方向模型跑通了结果也出来了别急着写“符合预期”就完事。你需要从结果里读出两层信息第一层是容量配置是否合理第二层是调度策略是否体现了系统特性。6.1 容量配置结果怎么读容量优化结果会给出风电机组装机、光伏装机、电解槽功率、储氢罐容量等最优值。你可以做几组对照去掉并网模式、只做离网看容量如何变化改变电解槽的单位投资成本看装机如何移动。这些敏感性分析是论文和工程报告中最重要的部分。我复现时发现并网模式下最优容量组合倾向于“风光大、电解槽小”因为风光零边际成本多发电可以卖钱而电解槽贵不如从电网买电制氢。离网模式下则相反电解槽会相对大一档储氢罐容量也会明显增加因为系统必须在无风无光期间靠储氢来保障连续供氨。这个规律可以作为检验模型是否合理的参考。6.2 调度策略的典型模式观察调度曲线你会发现几种明显的运行模式。风光充裕时电解槽满负荷运行多余电力卖电网并网或弃掉离网夜间光伏为零风电可能较小时电解槽降载到最低运行比例或者干脆停机氢气由储氢罐供出。这种“充放互补”的模式如果和你的工程经验匹配说明调度约束写得对。另一点要看储氢罐的利用率。假如最优结果里储氢罐基本上始终满着或始终空着说明储能容量没起到缓冲作用那要么是储氢罐成本太低导致随便装大要么是氢气供需曲线本来就匹配不需要缓冲。可以通过调整储氢罐单位投资成本看看调度曲线是否会发生变化。6.3 从确定性优化走向随机优化与鲁棒优化这个模型最自然的扩展方向是把风光的不确定性显式建模。当前模型用的是典型日确定性数据隐含假设风光预测完全准确。实际工程中风光预测误差很大所以很多新论文会引入多场景随机优化或分布鲁棒优化。场景法把典型日扩成几十个场景每个场景有概率权重目标函数变成所有场景期望成本最小约束在某些场景下可松弛。Cplex依然能有效地求解这类大规模MILP只是计算时间会成倍增加。另一个扩展方向是加入碳捕集或绿氨认证约束或者把氨作为终端产品直接卖给下游建立更完整的产品价值链条。如果对电化学模型更感兴趣还可以把电解槽的过载特性、热管理约束引入调度让模型更贴近实际设备。我在重复跑完这些扩展场景之后最大的感受是这个模型本质上是一个“多能互补系统的通用骨架”你换掉一组数据、改掉几个约束就能迁移到电制甲醇、电制甲烷、氢燃料电池热电联供等其他领域。所以花时间把这个骨架吃透后续做任何综合能源系统的优化项目都会事半功倍。最后再说一个小技巧做敏感性分析时批量跑场景记得用Matlab的并行计算工具箱parfor一次把几十组投资成本参数全部算完别一个接一个地等。