综合能源系统两阶段随机优化：源荷不确定性场景生成与容量配置实战

直接说结论：这篇论文复现的难度不在“两阶段随机优化”这个数学框架本身，而在“源荷不确定性”如何生成、如何缩减、如何嵌入优化模型而不让求解器直接卡死。我第一次跑通这个模型用了将近三周，中间踩了很多坑，走了不少弯路。这篇文章把我最终落地的一套完整方案拆开来讲，包括场景生成、两阶段建模、求解器配置、线性化处理，以及我实际排过的几个典型故障，希望能帮你少走几个月的弯路。

1. 问题拆解：这个模型到底在优化什么

1.1 容量配置和运行调度为什么必须放在一起

先说清楚一个问题：为什么这个模型要把容量配置和运行调度同时放进一个优化问题里，而不是分两步做——先定容量，再做调度？

我刚开始也这么想，甚至真去试过。先按典型日数据跑一个规划模型，得到风电、光伏、储氢、P2G设备的容量，然后固定这些容量，再跑一个调度模型。表面上逻辑没问题，但实际上结果非常失真。原因在于：容量配置决策会影响运行可行性，运行可行性反过来又决定容量配置的经济性，两者是强耦合的。举个例子，如果P2G的容量配置小了，风光大发时氢产量受限，弃风弃光严重，这时候就算运行模型优化得再好，整体经济性也救不回来。反过来，如果储氢容量配得过大，前期投资成本直接压垮收益，运行阶段再灵活也弥补不了。

联合优化的数学本质是在同一组变量空间里同时求解投资决策变量和运行决策变量。投资决策变量一般是连续的容量值（比如电解槽的额定功率、储氢罐的体积）加一批二元变量（比如某个候选设备是否投建）。运行决策变量则是小时级的功率流、储能SOC、氢流量、气流量。优化目标函数里同时含有投资成本的年值化项和运行阶段的期望成本项。

1.2 源荷不确定性的“源”和“荷”具体指什么

标题里“源荷不确定性”的“源”指的是风电、光伏的出力不确定性，“荷”指的是电、氢、气三类负荷的不确定性。这里有个很容易忽略的细节：综合能源生产单元的输出不是单一产品，而是电、氢、气三种能源产品，其中氢和天然气还可以通过储罐和管网缓冲。因此负荷侧的三种需求曲线本身就具有不同的波动特征和季节特性，不能简单套用同一个概率分布。

我在做场景数据时采用的思路是：用历史统计参数构造风电、光伏的预测误差分布和三类负荷的预测误差分布，然后通过蒙特卡洛采样生成大量短期场景，再用K-means聚类缩减出具有代表性的场景子集。这里每个场景都带有一个概率权重，两阶段期望成本就是对所有场景的概率加权求和。

1.3 两阶段随机优化的核心思想

两阶段随机规划的标准形式是：

目标函数为第一阶段投资成本加上第二阶段期望运行成本的最小化。对应到这里：

• 第一阶段（here-and-now）决策：决定各个设备的容量，比如风电装机、光伏装机、电解槽功率、燃料电池功率、储氢罐容量、电锅炉容量、吸收式制冷机容量等。这些决策必须在不确定性实现之前做出，一旦确定就不能改变。
• 第二阶段（wait-and-see）决策：在给定容量配置和某个具体场景下的风光出力、负荷需求后，决定逐时段的设备出力、储能充放电量、母线功率交换等，使该场景下的运行成本最小。
• 非预期约束（non-anticipativity）：第一阶段决策对所有场景必须是同一个值，这一点在建模时容易被漏掉。我见过不少复现失败的情况就是把容量配置变量写在了每个场景内部，导致求解器每个场景得到一套容量，最后退化成多个独立优化问题，完全失去了两阶段的意义。

1.4 设备集和能流关系的基本框架

综合能源生产单元的设备集合通常包括风力发电机组、光伏阵列、电解水制氢设备（P2G）、储氢罐、氢燃料电池、燃气轮机组、天然气锅炉、电储能、吸收式制冷机（如果涉及冷负荷）、P2G副产热回收等。这个单元的输入是风光资源和天然气，输出是电、氢、气，必要时还有冷/热。

能流关系的核心是一组耦合矩阵方程，本质上是各母线（电母线、氢母线、气母线、热母线）的功率平衡方程。例如电母线平衡方程为：风电出力 + 光伏出力 + 燃气轮机发电 + 燃料电池放电 + 储能放电 + 外购电 = 电负荷 + 电解槽耗电 + 电锅炉耗电 + 储能充电 + 卖给电网的电。氢母线和气母线类似，需要注意的是气母线里除了天然气网购之外，P2G产氢和甲烷化产气之间还有可能通过氢转气（Power-to-Methane）进一步耦合。

我把这套能流关系在Matlab里维护成一个稀疏矩阵结构，每个设备对应一条或多条能流路径，母线平衡约束通过线性等式实现。这个结构的好处是后期加设备或者改场景数时不需要重写平衡方程，只需增减矩阵的行列。

2. 场景生成与缩减：模型精度的分水岭

2.1 不确定性参数的概率模型设计

不确定性参数建模这一步直接决定整个模型的可信度。我的做法是把风电出力视作“预测值+随机误差”的形式，而不是直接对出力本身做分布拟合。

• 风电出力基准值由Weibull风速分布通过功率曲线转换得到，误差项采用正态分布，标准差取预测值的15%-20%。
• 光伏出力基准值由Beta分布的太阳辐照度配合光伏转换效率计算，误差项同样用正态分布，标准差取预测值的10%-15%。
• 三类负荷（电、氢、气）的基准值采用历史典型日曲线，误差项用多元正态分布建模，并通过相关系数矩阵反映负荷之间的时空相关性。

为什么要用“基准值+误差项”而不是直接采样总出力？原因有两个。第一，预测误差模型更贴近实际调度场景——调度员在安排日前计划时掌握的是预测值，不是真实值。第二，误差项相对基准值来说数值范围小得多，采样后的场景不会出现“风电负数”这种物理上不合理的极端值。

2.2 蒙特卡洛采样+聚类缩减的实现细节

先采样，再缩减，这是目前的主流做法。我当时的操作流程如下：

1. 设定采样场景数Nsample = 2000。
2. 对每个场景，随机抽取24小时的风速误差、辐照度误差、负荷误差序列，叠加到基准曲线上，得到完整24小时场景数据。
3. 对2000个原始场景做K-means聚类，聚类距离采用欧式距离，特征向量为所有不确定量的24维序列拼接。
4. 用肘部法则确定聚类数，最终取K = 20或K = 30。每个聚类中心作为典型场景，场景概率等于该簇样本数占比。

缩减后的场景集要同时具备两个性质：一是概率之和为1，二是能覆盖原始采样空间的主要波动范围。我在实测中发现，场景数从5个增加到20个，优化结果的期望成本变化很大，但从20个增加到50个，变化幅度不超过3%。这说明20个场景对于这个模型已经足够。

场景缩减后还应该做一次“场景平滑”处理。聚类中心经常会比原始样本更平坦，直接使用聚类中心会让模型低估波动性，导致容量配置结果偏小。我的做法是：在聚类中心基础上，随机叠加一个幅值为原标准差20%-30%的小扰动，重新生成多个扰动场景并平均，这样既保留了聚类中心的代表性，又不至于过度平滑。

2.3 场景树的构建与非预期约束

两阶段模型中场景树的结构很简单：第一阶段是一个根节点，第二阶段是K个叶子节点，每个叶子节点对应一个场景和一组运行变量。非预期约束在这个树结构下不需要显式写出来——因为第一阶段变量本身就是全局共享的。但如果你把模型写成了逐场景独立的子问题，那就必须显式添加非预期约束来强制各个场景的第一阶段变量相等。

实现上有两种方式：

• 方式A：所有场景的变量全部放在一个模型里，容量变量下标不含场景维度，天然满足非预期约束。这种方式简单直观，适合中小规模问题。
• 方式B：按场景拆分成K个子问题，每个子问题有自己的容量变量，再额外添加容量变量相等的等式约束。这种方式适合用Benders分解或拉格朗日松弛求解，但直接调用求解器时反而更慢。

我推荐用方式A。对于本文的模型规模——比如候选设备8个、时段24、场景20个——变量总数大约是20×24×(20~30个运行变量)+20×(8个容量变量)+若干二元变量，规模大概在1万到2万之间，用YALMIP建模后交给CPLEX或Gurobi求解完全没问题。

3. 两阶段随机优化模型的数学框架

3.1 目标函数的具体形式

目标函数我建议写成年度总成本最小化，包含四项：设备投资年值成本、运行维护成本、从外部电网购气和购电成本、弃风弃光惩罚成本（可选），然后减去可能的售电收益。

投资成本年值化处理公式是：等年值系数r×(1+r)^n / ((1+r)^n - 1)，其中r是贴现率，n是设备寿命。风电、光伏、电解槽、燃料电池的寿命不同，要分开计算。我之前见到过有人对所有设备统一用20年寿命和8%贴现率，结果容量配置严重偏向寿命长的设备。正确做法是每个设备单独查寿命：光伏25年、风电20年、电解槽10~15年、燃料电池10年、储氢罐20年、电储能10年。

运行阶段的成本项需要对每个场景逐时段求和，再按场景概率加权。这里所有成本都是线性的，所以整个目标函数是线性表达式，不会引入求解困难。

3.2 约束体系：母线平衡、设备模型和储能的SOC逻辑

约束体系分三层来写，清晰且不容易漏。

第一层是母线平衡约束。前面已经提过，电、氢、气三条母线的逐时段功率平衡。这个约束在每个场景的每个时段都必须严格成立。

第二层是设备模型约束。比如电解槽的输入电功率和产氢速率之间的线性转换关系，燃料电池的产电功率和耗氢速率之间的线性转换关系，燃气轮机的产电和耗气关系。需要注意的是转换效率取常数还是分段线性函数。实测下来，对于容量配置阶段的研究，取常数效率足够——因为模型的时间分辨率是小时级，设备的动态特性在小时尺度上体现得不明显。当然，如果你想做得更精细，可以把电解槽效率处理成关于负载率的分段线性函数，但这会大幅增加模型的非线性和求解时间，需要权衡。

第三层是储能约束。储氢罐的SOC（荷电状态）更新方程：SOC_t = SOC_{t-1} + 充入量 × 充入效率 - 放出量 / 放出效率。储能容量上下限约束、充放速率约束、始末SOC相等约束（循环约束）。对于电储能，还要额外加一组充放电状态互斥的二元变量，防止同时充放电。储氢罐通常不做互斥要求——因为氢气系统本身有缓冲能力，同时充放虽然不常见但并不物理违背，加互斥反而会让模型更紧、容易无解。

3.3 不确定性条件下的约束处理方式

对含不确定性的约束，比如“电网交互功率不超过联络线容量”，在随机规划框架内有三种处理方式：

• 方式1：对所有场景都强制满足。这是最保守的，也是最常见的做法。意味着联络线容量必须能应对所有场景下的极端潮流，容量配置结果会偏大。
• 方式2：引入机会约束，允许违反概率不超过某个阈值（比如5%）。这需要额外引入二进制变量来标记哪些场景允许越限，问题变成混合整数规划。这个改法在工程上有意义，但会增加求解复杂度。
• 方式3：把约束从确定性约束改为期望值约束，允许个别场景越限，只惩罚期望越限量。

我最终采用的是方式1加弃风弃光惩罚。因为两阶段随机优化的目的本来就是用场景集来覆盖不确定性，如果再把约束做成机会约束，整个模型的复杂度会上升一个量级。对于复现论文的场景，方式1足够。

3.4 二元变量与线性化处理

模型里会出现哪些二元变量？我总结下来有三个来源：

• 设备投建二元变量：某个候选设备是否建设。
• 电储能充放电状态变量：防止同一时段同时充放电。
• 分段线性函数的分段选择变量：如果设备效率或成本是分段线性的。

前两类必需要，第三类可以作为简化选择的备选项。如果采用常数效率，第三类就不需要。

在线性化方面有一个很关键的技巧：投资容量和设备投建二元变量之间的耦合约束。例如电解槽的额定功率P_P2G = 投建变量 × 候选容量，这个等式在数学上是非线性的（二元变量乘连续变量）。处理方法是用大M约束来线性化：

P_P2G ≤ M×z，P_P2G ≤ 候选容量，P_P2G ≥ 候选容量 - M×(1-z)。

这里还有一个容易被忽视的细节：容量配置变量如果取连续值，那么线性化之后，未投建设备的容量会被逼到0。但如果你给投建变量加了“最小投建容量”约束，比如要么不建，要么至少建10MW，那处理起来稍微复杂一点。常见做法是引入分段常数候选容量集合，用整数变量从离散集合中选一个容量档位。我的建议是：尽量让建设容量保持连续，不要加最小建设约束，这样模型最干净，求解最快。

4. 从零搭建到跑通：完整实现路径

4.1 环境选择和求解器配置参数

我在复现时用的是Matlab R2022a + YALMIP + CPLEX 12.10。也可以换成Gurobi，区别不大。有一个值得说的经验：YALMIP在建模大规模问题时很方便，但有一个坑——默认情况下它对二阶锥约束的处理不如Gurobi直接。这个模型里基本不会出现二阶锥约束（因为储能模型的SOC方程是线性的），所以问题规模还算友好。

求解器参数方面，我做了如下配置（以CPLEX为例）：

• MIP相对间隙（MIPGap）：设置为0.01，也就是让求解器在1%最优性间隙内停止。如果设到0.0001，模型可能要跑几个小时，而结果改善微乎其微。
• 线程数（Threads）：设为本机物理核心数。
• 时间上限（TimeLimit）：设置为7200秒，超过就接受当前最优整数解。
• 节点文件（NodeFile）：如果机器内存不够，设置节点文件路径以便在内存不足时降级到磁盘。

这套配置下，我跑的模型规模约为1.5万个变量、2万条约束、200个二元变量，求解时间大约10-20分钟。如果超过1小时还没解完，大概率是模型结构有问题，不是求解器不够快。

4.2 YALMIP建模的关键代码逻辑与结构

我给出一个简化的建模骨架（不写完整表达式，只展示结构）：

% 决策变量 X_cap = sdpvar(n_device, 1); % 容量配置变量 Z_build = binvar(n_device, 1); % 设备投建0-1变量 for s = 1:n_scenario % 逐场景运行变量 P_grid{s} = sdpvar(T, 1); % 电网交互功率 P_p2g{s} = sdpvar(T, 1); % 电解槽输入电功率 Q_h2{s} = sdpvar(T, 1); % 产氢速率 SOC_h2{s} = sdpvar(T+1, 1); % 储氢罐SOC % ... 其他运行变量 % 逐场景约束（母线平衡、设备模型、储能SOC） Constraints = [Constraints, 母线平衡约束]; Constraints = [Constraints, 设备模型约束]; Constraints = [Constraints, SOC更新和边界约束]; end % 目标函数 Objective = 投资年值成本 + sum(prob(s) * 运行成本(s)); % 求解 optimize(Constraints, Objective, sdpsettings('solver', 'cplex', ... 'cplex.mip.tolerances.mipgap', 0.01));

这里最需要注意的一个工程问题是：场景之间的约束和变量要独立，不要交叉引用。如果场景s里误用了场景s'的变量，两阶段模型的场景独立性就被破坏，优化结果会严重失真且难以发现。我在代码里专门写了断言逻辑，在求解前检查每个场景的约束是否只引用了本场景变量。

另一个工程技巧是目标函数的构建速度。初学者经常用循环拼接大段表达式，导致YALMIP建模时间很长。我的做法是：把目标函数的系数向量化，用向量点积构建，避免在循环里反复增加表达式项。实测下来，建模时间从几分钟降到十几秒。

4.3 求解结果的后处理与正确性验证

跑完模型之后不能直接看结果，至少要做三个方面的验证：

第一，可行性验证。检查所有关键约束的残差是否在允许范围内。尤其是母线平衡约束和SOC终值约束。SOC终值如果不等于初值，说明循环约束没有生效，需要检查模型是否真的包含了最后一个时段到第一个时段的衔接方程。

第二，场景一致性验证。对每个场景逐一检查：给定该场景的容量配置（应为全场景同一值），固定第一阶段变量，单独对该场景求解运行子问题，目标函数值应该等于原模型里该场景的期望项除以概率。如果两者不相等，说明场景之间发生了耦合，或者目标函数里概率权重写错。

我实际遇到过一次这个情况：场景概率在for循环里写成了1/n_scenario + 0.05，结果每个场景的权重都不是真实权重，目标函数严重偏斜。这种错位用上面对比法很容易发现。

第三，经济指标汇总。输出各设备容量、总投资成本、年均运行成本、弃风弃光率、最优场景下的设备利用小时数。这些指标有助于判断结果的工程合理性。比如如果你发现燃料电池年利用小时数只有200小时，但投资成本很高，那大概率是容量配置不合理，需要检查模型是否遗漏了某些约束，导致燃料电池成为“可有可无”的设备而恰好又被选中。

5. 求解中的常见故障：三条典型排错路径

5.1 模型无解或不可行：从松弛分析找根因

无解是复现这类模型最常见的坑。我的排查路径是这样的：先把所有整数变量固定到初始可行解（比如固定不建设新设备），看纯线性问题是否可行。如果线性问题都不可行，那就是连续变量约束之间产生了矛盾。

最常见的矛盾来源是储能SOC方程的时间衔接。我曾经在SOC递推方程里把充入效率写在了放出项上，导致一个时段内能量不守恒，24小时循环后SOC偏差巨大，模型直接报不可行。排查办法是：先只保留单场景，把时段数从24缩减到2，手工推演一遍SOC轨迹，看是否满足能量守恒。

5.2 求解时间爆炸：多数原因指向大M取值

模型求解了3小时还不收敛，先检查所有大M约束中的M值。我经常看到有人把M设成1e6，以为“足够大”，但实际上过大的M会给线性松弛引入极差的数值条件，CPLEX和Gurobi会在分支定界时产生灾难性的数值误差，导致节点数爆炸。

正确做法是：对每个约束单独设置最紧的大M值。比如设备功率上限约束中的M，取该设备候选容量的最大值即可；联络线容量约束中的M，取联络线最大允许功率即可。M宁可多算一步，也不要偷懒用一个全局大数。

5.3 结果反直觉：先怀疑场景问题，再怀疑模型问题

如果模型能求解，但结果不符合预期——比如光伏容量配得异常大，或者储氢罐容量配为零——不要先改模型结构，先把场景数据的分布特征画出来看。

我的实测经验是：50%以上的“模型不合理”问题，根因在场景生成阶段。比如K-means聚类时用了未经归一化的数据，导致量纲大的特征（电价负荷）主导了聚类距离，量纲小的特征（光伏辐照）被忽略。结果是场景集虽然看起来多样，但光伏场景的波动完全没被代表，优化结果自然不肯配置光伏。

正确做法是在聚类前对所有特征做Z-score标准化。

5.4 目标函数不符：概率权重与成本口径

这个坑比较隐蔽。两阶段模型的期望成本需要在第二阶段对场景概率加权求和，但如果你把第一阶段投资成本也写进了场景循环里，就会产生“每个场景都算一遍投资成本”的错误，期望成本变成单场景投资成本的N倍。

我的检查方法是：单独固定第二阶段运行变量为零，看目标函数是否等于第一阶段投资成本。如果不等于，说明投资成本和场景耦合了。这个测试在建模完成后立刻做一遍，几十秒钟就能验证。

6. 容量配置结果的分析方法与调参方向

6.1 从结果中提取投资组合的敏感性规律

模型跑通后，我习惯先做一组“基准场景”试验：把随机场景替换为单一确定性场景（即所有误差设为零），对比两阶段方案和确定性方案的结果。这个对比非常有价值，它能直观量化不确定性带来的成本代价——两阶段方案的投资成本通常会略高，但运行成本更低，总期望成本低于确定性方案。

如果两阶段方案的总成本反而高于确定性方案，通常说明场景数过少或者场景误差设置过小，不确定性没有被充分体现。

6.2 参数敏感性分析实验设计

建议做三组敏感性分析：

• 场景数敏感性：固定其他参数，将场景数从5逐步增加到50，记录期望成本和求解时间的变化。
• 误差幅度敏感性：将风电/光伏预测误差标准差从10%逐步增加到30%，观察容量配置的变化方向。
• 贴现率敏感性：从6%到12%变化，观察投资组合中不同寿命设备的取舍变化。

每组敏感性分析结束后，把结果整理成表格：场景数vs期望成本、误差水平vs风电装机容量等。这些表格是论文的“证据面”，也帮助你判断模型的稳健性。

这里有个经验：如果误差幅度增加20%，但最优容量配置几乎不变，说明模型对不确定性不敏感，这时候要回去检查场景生成是否真的包含了波动信息。我试过一种情况：在生成场景时把误差序列的相关系数矩阵设成了对角阵，导致各个时段的误差完全独立，整体波动性被滤波效应大幅削弱，模型自然表现不出对不确定性的反应。

6.3 设备耦合关系对容量决策的影响

综合能源生产单元的容量配置结果最值得分析的是设备之间的耦合。以P2G和储氢罐为例：如果电解槽的候选容量大、成本低，储氢罐候选成本高，那么最优解往往是“大电解槽+小储氢罐”——电解槽跟随风电出力灵活制氢，氢气直接送入管网，储氢罐只起少量缓冲作用。如果储氢成本低而电解槽成本高，则倾向于“中电解槽+大储氢罐”——电解槽尽量满负荷运行，多余氢气存入储罐等待高需求时段释放。

这种设备间的替代关系反映在求解结果上，就是容量配置比例的变化。分析这个替代关系时，我习惯做CapEx（单位容量投资成本）和OpEx（单位电量运行成本）的对比热力图，横轴是电解槽成本，纵轴是储氢罐成本，颜色是配置比例。

7. 复现过程中的额外经验与建议

7.1 和论文原始结果对不上怎么办

复现论文时最常见的问题是结果和原文数值对不上。我的经验是：不要追求逐位一致，先定位“量级一致”和“趋势一致”。论文里的具体数字高度依赖场景生成随机数种子、历史数据范围、边界条件设定，甚至求解器版本都会导致微小差异。你需要对比的是：容量配置的排序是否一致、成本构成的量级是否一致、设备组合的取舍趋势是否一致。

如果趋势都不一致，优先检查隐含假设。比如论文可能假设氢负荷可以在管网的缓冲下不满足短时波动，而你建模时对氢负荷做了逐时段刚性平衡——这两种假设会导致完全不同的储氢罐容量结果。

7.2 模型扩展的方向：从两阶段到多阶段

跑通两阶段模型后，下一步大概率是往多阶段延伸。这里给一个提示：两阶段到多阶段的跨越，难点不在数学形式，而在场景树的管理。多阶段需要每阶段都做一次分支，场景数阶乘增长，需要提前设计缩减策略。

我倾向于用滚动时域（rolling horizon）方案替代完整的多阶段随机规划：把模型在日尺度上保持两阶段结构，但在跨日边界用滚动更新的方式传递状态量。这种方式虽然没有严格的多阶段随机规划最优性，但工程实用性强，求解时间可控，而且和实际调度流程天然对接。

7.3 数据准备的清单管理

最后整理一份我在复现中使用的数据准备清单，你可以直接对照检查：

数据质量决定模型质量，这句话在综合能源系统优化里体现得尤为明显。头几次跑出来的结果不合理，超过一半都在源数据的细节上——比如电价曲线是夏季的还是冬季的、负荷曲线是工作日还是节假日，这些标签必须和数据值一起管理，否则场景生成环节会把完全不同的季节特征混合在一起，产生物理上不可能出现的场景组合。

这次复现最大的收获不是熟练了两阶段随机规划的求解流程，而是意识到“不确定性建模”这件事在综合能源领域本质上是一个概率分布工程——你需要花至少一半的时间在场景数据的生成、缩减和验证上，模型求解反而是相对确定性的工作。如果你正在复现这个模型并且卡在某一步，按照上面的排查路径走一遍，大概率能定位到问题所在。