关键场景辨别驱动的两阶段鲁棒微网调度与CCG实现

接手园区微网调度项目那阵子，光伏装机占到全园区负荷的将近一半，天气好的时候发电量一路冲到上限，云层一厚又瞬间掉下去，柴油发电机启动需要时间，储能容量又有限。最初用确定性优化模型做日前调度，预测曲线和实际偏差一大，当天实时调整就手忙脚乱，甚至出现过几次切负荷。后来我换成了两阶段鲁棒优化框架，结果又遇到新问题：蒙特卡洛抽了几百个不确定性场景，每次迭代主问题规模巨大，模型跑一个调度周期要二十多分钟，根本无法滚动执行。

最后是“关键场景辨别算法”帮了大忙。它的思路很朴素：与其把几百个场景全部塞进两阶段鲁棒模型，不如先通过聚类和违约度排序，挑出少数几个对调度决策最有威胁的关键场景，再用列约束生成（C&CG）迭代求解。整套流程用Matlab和Yalmip实现，单次调度时间从二十分钟压缩到三分钟以内，而且保守性损失很小。这篇文章把我完整的实现思路、核心代码骨架和踩过的坑写清楚。

如果你也在做微网优化调度、鲁棒优化或者可再生能源消纳相关的研究，这篇内容可以直接作为入门参考，尤其是“关键场景怎么辨别”“C&CG怎么落地”“Matlab代码怎么组织”这三块，都是有实操价值的东西。

1. 为什么确定性调度模型在真实运行中总是差一口气

微网优化调度，基础逻辑是“预测-计划-调整”。给定负荷曲线、光伏出力曲线和电价曲线，优化出柴油机出力和储能充放电计划。但真实运行的困难在于，预测值从来不会刚好等于实际值。

1.1 微网不确定性到底来自哪些维度

微网里的不确定性来源比大电网更密集，因为可调节资源少、惯性小，任何一个环节的波动都可能让计划失效。按我的经验，至少有这么几个维度：

• 光伏出力：受云层、温度、空气质量影响，分钟级波动就能超过额定出力的30%。日前的光照强度预测，哪怕用数值天气预报，误差也经常在15%以上。
• 风电出力（如果有）：风速的随机性更强，尾流效应也让单机出力关系复杂。
• 负荷波动：园区负荷受生产计划、天气、人员行为影响，尤其是空调负荷，夏季下午的短时爬坡很凶。
• 上级电网购电价格：如果参与现货市场或峰谷电价，价格本身就是一个不确定性参数。
• 设备故障和通信延迟：这个容易被忽略，但在实际运行中会导致控制指令无法按时执行。

这些不确定性如果只用“预测值+固定备用”来处理，往往不是过于乐观就是过于保守。乐观了，遇到极端天气就调节不过来；保守了，柴油机持续高功率运行，经济性很差。

1.2 确定性模型和随机期望模型的共同盲区

确定性调度模型是最常见的一类，它的目标函数一般是：

min 总运行成本 subject to 功率平衡、储能动态、机组出力上下限等

这类模型里，光伏出力和负荷都是固定值。随机规划则进了一步，把不确定性建模成多个场景，目标函数变成“期望成本最小”。随机规划能处理概率分布已知的情况，但它本质上是在“平均情况”下优化，对极端场景的容忍度不够。实际项目里，最可怕的不是平均偏差，而是小概率大偏差事件——比如连续阴雨天叠加负荷高峰，这种场景在期望优化里权重很低，模型不会为它专门预留充足的安全裕度。

鲁棒优化干脆换了个思路：我不求在所有场景下最优，但我要保证在任何可能的场景下都可行。用行话说，这是从“期望视角”切换到“最坏情况视角”。两阶段鲁棒优化则是把决策拆成两段，第一阶段做日前预决策（机组启停、储能充放电计划等调整代价高的变量），第二阶段是在不确定性揭晓后的实时再调度（调整出力、切负荷等快速动作）。

1.3 两阶段鲁棒优化的基本结构

两阶段鲁棒优化的一般形式可以写成：

min_x max_u min_y f(x, u, y)

外层是第一阶段决策 x，中间层是不确定参数 u 在最坏情况下的取值，内层是第二阶段决策 y。这个三层结构读起来抽象，但实际意义很清晰：

• 第一阶段：今天根据预测信息，决定明天柴油机开不开机、储能是否允许充电等；
• 中间层：明天实际光伏、负荷朝最不利方向发展；
• 第二阶段：在最不利情况下，用最经济的调节手段维持平衡。

这个结构的好处是不需要精确的概率分布，只需要不确定性集合 U。U 可以是简单的盒式区间，也可以带预算约束的切片箱式集合。相比随机规划动辄几百上千个场景，鲁棒优化对数据的要求低得多，正好符合微网项目实际可获取数据量有限的特点。

2. 关键场景辨别算法的原理与实现流程

两阶段鲁棒优化直接求解是很难的，因为 min-max-min 是嵌套结构，内层 max 和内层 min 互换非常麻烦。主流做法是用 C&CG 把它分解成主问题和子问题迭代求解。但主问题一开始不知道哪些不确定性场景最危险，如果拿一个巨大的连续不确定集直接建模，模型不可解。这时候就需要“关键场景辨别算法”。

2.1 从海量蒙特卡洛场景到关键场景集

场景生成的常用方法是蒙特卡洛采样。比如用 Beta 分布拟合光伏出力的历史规律，用正态分布拟合负荷预测误差，随机生成 N=500 个未来 24 小时的场景。每个场景是一条 24 维（或 96 维，15分钟一个点）的时序曲线。

但把 500 个场景全部代入 C&CG 主问题，主问题就有 500 组约束和 500 组第二阶段变量，模型规模直接膨胀。而且这些场景之间有大量重叠信息，它们对调度决策的影响基本上是重复的。场景辨别算法的目标就是：从 N 个场景里选出 K 个代表性场景（K远小于N），使得这 K 个场景下的调度决策，在其它的场景下大概率依然可行，而且目标值接近完整场景集的结果。

2.2 两种关键场景筛选思路：聚类中心法与违约度排序法

我在项目中测试过两种筛选思路，各有适用场景。

第一种是我最后采用的“聚类中心+边界样本”法。核心是用 K-medoids 或 DBSCAN 把 500 个场景聚类成若干个簇，然后每个簇里取“中心场景”和“边界场景”。中心场景代表这个簇最典型的形态，边界场景代表这个簇里距离中心最远、最有可能触发恶劣运行的形态。把中心场景放进模型，能保证决策对大多数情况不偏不倚；把边界场景放进去，能保证决策对极端情况有足够裕度。

聚类的特征空间需要做一点处理。如果直接用 24 维原始序列做欧氏距离，容易受整体功率水平主导，而忽略了曲线形态差异。我试过的做法是：先对 24 小时场景做 PCA 降维，保留贡献率 95% 以上的主成分，再在降维空间里用欧氏距离聚类。这样曲线形态相近的场景更容易被聚到一起，关键场景的代表性更强。

第二种是“违约度排序”法。对每个采样场景，先固定一个第一阶段基准决策，求解对应场景下的第二阶段经济调度，记录每个场景下约束违反量或者运行成本异常高的情况。把目标值最高、或者约束修改成本最大的 Top-K 场景挑出来，作为关键场景。这种方法的好处是完全面向调度目标，不用考虑聚类几何；缺点是要先有基准决策，一般需要预跑一轮确定性优化。

为了兼顾效率和鲁棒性，我的实际做法是两者结合：先用聚类中心法生成一批关键场景，保证场景覆盖度；然后跑一轮 C&CG 迭代，把每一轮子问题识别出的最恶劣场景再动态加入关键场景集。这正是“辨别”二字的含义——关键场景不是一次性生成的，而是随着迭代不断更新的。

2.3 关键场景辨别算法的完整伪代码

下面这段伪代码是我实际实现时采用的流程，读者可以直接对照翻译成 Matlab 代码。

输入: uncertaintyModel # 光伏、负荷不确定性参数分布 N # 初始采样场景数，本文取500 K # 关键场景最大数量，本文取20 epsGap # C&CG收敛间隙，本文取1e-3 输出: x_opt # 第一阶段最优决策，例如机组启停、储能计划 1. 生成初始场景集 scenes = sampleScenes(uncertaintyModel, N) 2. 场景预处理: PCA降维 + 归一化 feat = normalize(pca(scenes, 0.95)) 3. 聚类，得到K个簇 clusters = kmedoids(feat, K) 4. 从每个簇中提取中心样本和边界样本 keyScenes = [] for each cluster: center = medoid(cluster) boundary = sample_farthest_from_medoid(cluster) keyScenes.add(center) keyScenes.add(boundary) 若 keyScenes 数量少于K，补充簇内离散度高的样本 若 keyScenes 数量超过K，按目标值敏感度排序保留前K个 5. 初始化C&CG LB = -inf, UB = inf, gap = inf, iter = 0 6. while gap > epsGap: iter = iter + 1 求解主问题MP: min 第一阶段成本 + theta s.t. 第一阶段约束 对每个 keyScene in keyScenes: 第二阶段运行约束 + theta下界约束 获得 x_cur 和 LB 求解子问题SP: 给定 x_cur，在不确定集U中寻找最恶劣场景worstScenario 并计算对应第二阶段最小运行成本 f_sp 若 f_sp < UB: UB = f_sp 计算 gap = |UB - LB| / |UB| 若 worstScenario 不属于当前 keyScenes: keyScenes.add(worstScenario) # 动态辨别关键场景 若 keyScenes.size() > K_max: 根据场景出现频率和对目标影响度裁剪 7. 返回 x_opt

这段流程的关键是第 4 步和第 6 步的最后一行。前者保证初始场景集有代表性，后者保证 C&CG 迭代中真正恶劣的场景不会被丢到模型外面。

3. 列约束生成（C&CG）算法：两阶段鲁棒问题的主循环

关键场景集准备好之后，剩下的就是两阶段鲁棒优化的核心求解器。目前工程中最常用的算法是列约束生成，比传统的 Benders 分解在微网这类问题上收敛更快、数值表现更稳定。

3.1 min-max-min 结构如何拆分成主问题和子问题

把刚才的三层结构拆开：

• 主问题（MP）处理第一阶段决策 x 和辅助变量 theta。它假设不确定参数只能取当前关键场景集里的若干个离散值，所以是一个普通的混合整数线性规划（MILP）。
• 子问题（SP）给定第一阶段决策 x，让不确定参数 u 在连续不确定集 U 中找最坏取值，然后计算第二阶段的调节成本。

主问题求解得到的是原问题的一个下界，因为在不确定性上我们只考虑了有限个关键场景，相当于“放窄”了不确定性的范围，最优成本肯定不高于真实最坏情况下的最优成本。子问题求解得到的是原问题的一个上界，因为它是针对某一个 x 的最坏情况，不一定是最优 x 下的最坏情况，所以成本不低于真正的最优值。

C&CG 的逻辑就是反复用主问题产生的 x 去刺激子问题，找出更恶劣的场景加进主问题，直到上下界间隙收敛到要求精度。

3.2 子问题中的 max-min 转化与不确定性集合构建

子问题内部还有一个 max-min 嵌套。幸好，对固定 x 和固定 u 来说，第二阶段是一个线性规划（LP），满足强对偶条件，可以把内层 min 对偶成 max，于是整个子问题变成一个单层的 max 问题：

SP(x*) = max_{u, λ} constant(x*, u) + λ^T * RHS(u)

这个单层问题目标函数里包含 u 和 λ 的乘积项，如果 u 是连续变量，就出现双线性项，需要用大M法逐段线性化，或者用迭代逼近技巧处理。我在微网算例中采用的做法是：把不确定集 U 限定为预算约束下的盒式集合，即每个时刻的 u 最多偏离预测值一定比例，同时所有时刻的总偏离量受到预算参数 Gamma 限制。类似问题在电力系统文献中有多种处理方式，本文实现中为控制难度，对子问题双线性项采用了 big-M 离散化，效果稳定。

不确定集的构造看起来是细节，其实直接决定优化结果的保守度。Gamma 越大，允许的偏差总和越大，结果自然越保守。我在实际算例里对 Gamma 做了灵敏度扫描，从 4 到 12（24小时尺度）都跑过，最后选了一个既能覆盖 90% 以上历史极端场景、又不会让柴油机时刻高负荷运转的值。

3.3 收敛判断与加速技巧

C&CG 的迭代终止判据很多教程写成上下界间隙小于 1e-3 或 1e-4，但实际中千万别一上来就求高精度。微网调度问题中，第一阶段变量包含大量0-1整型变量，主问题本身就是MILP，每轮迭代求解时间都不短。如果间隙设成 1e-4，可能要多跑十几轮，而这些额外轮次对最终调度方案的改善非常有限。

我的经验是：

• 间隙设成 5e-3 到 1e-2 就够用了，对应误差大约几十到一两百块钱，对运行调度来说完全可接受；
• 给子问题设置合理的求解时间上限，Gurobi 里用timelimit参数；
• 主问题求解时把前一轮的整数解作为 warm start 传给求解器，能明显加快后续迭代；
• 如果关键场景数量接近上限，优先用新增最恶劣场景替换掉历史迭代中从未“起作用”的旧场景，保持主问题规模可控。

4. Matlab环境下基于Yalmip/Gurobi的实现细节

这一章写给想在 Matlab 里复现代码的读者。我不打算贴完整源码，因为每个微网拓扑和参数不同，完整代码不具备通用性，但核心骨架和关键细节必须交代清楚。

4.1 环境配置与求解器选型

我用的是 Matlab R2022b + Yalmip（最新版）+ Gurobi 10.0。你如果用的是 Cplex 或 Mosek，也完全可以，只是下面的求解器参数写法略有不同。

有一个特别重要的配置：Gurobi 安装后必须在 Matlab 里调用gurobi_setup完成路径配置，并在环境变量里设置GRB_LICENSE_FILE指向你的 license 文件。这一步卡住过很多人，运行solve时如果报license error，优先检查环境变量，而不是 Yalmip 安装。

Yalmip 的安装很简单，把整个文件夹放进 Matlab 路径即可。但注意版本兼容性：太老的 Yalmip 对 Gurobi 10 的支持不完整，会提示无法识别求解器。把 Yalmip 更新到 2023 年之后的版本可以避免很多兼容性问题。

4.2 核心代码骨架：关键场景生成与C&CG主循环

我把代码分成三个部分。第一部分是场景生成和聚类，第二部分是主问题建模，第三部分是子问题建模和迭代。

场景生成与聚类部分：

% 加载不确定参数模型 [T, N] = deal(24, 500); scenes = sampleScenes(N, T); % 自定义函数，返回 N*T 矩阵 % PCA降维 + K-medoids聚类 [coeff, ~, ~] = pca(scenes, 'NumComponents', ceil(0.95*min(size(scenes)))); feat = normalize(scenes * coeff, 'range'); [idx, C] = kmedoids(feat, 20, 'Distance', 'sqeuclidean'); % 提取中心场景与边界场景 keyScenes = extractKeyScenes(scenes, idx, C);

注意kmedoids返回的质心点 C 已经是降维空间的坐标，需要映射回原始空间找对应场景，不能直接拿去建约束。

主问题建模部分，用 Yalmip 定义变量和约束：

x = binvar(nG, 1, 'full'); % 机组启停，nG=柴油机数量 s = sdpvar(nS, 1, 'full'); % 储能充电功率等第一阶段变量 theta = sdpvar(1, 1); % 辅助变量 constraints_MP = []; obj_MP = c_first' * [x; s] + theta; for k = 1:size(keyScenes, 1) y_k = sdpvar(nY, 1, 'full'); % 当前关键场景下的第二阶段变量 u_k = keyScenes(k, :)'; constraints_MP = [constraints_MP, power_balance_mp(x, s, y_k, u_k), operation_limits_mp(x, s, y_k), theta >= c_second' * y_k]; end

子问题建模部分要稍微说明一下。假设在 C&CG 主循环里我们已经得到第一阶段的变量数值x_val和s_val，然后构建第二阶段 LP 并求对偶。直接调用dual函数就能得到对偶变量，再把对偶问题里涉及 u 的双线性项做线性化。

C&CG 主循环的骨架如下：

LB = -1e6; UB = 1e6; gap = 1; iter = 0; keyScenes = initialKeyScenes; while gap > 1e-3 && iter < 50 iter = iter + 1; % 求解主问题 optimize(constraints_MP, obj_MP, sdpsettings('solver','gurobi')); LB = max(LB, value(obj_MP)); x_val = value(x); s_val = value(s); % 求解子问题，得到最恶劣场景和对应成本 [worstCost, worstScene] = solveSP(x_val, s_val, uncertaintySet); UB = min(UB, worstCost); gap = abs(UB - LB) / max(abs(UB), eps); fprintf('iter=%2d, LB=%.4e, UB=%.4e, gap=%.4f\n', iter, value(LB), value(UB), gap); % 关键场景动态加入 if distance(worstScene, keyScenes) > 1e-6 keyScenes = [keyScenes; worstScene(:)']; if size(keyScenes, 1) > 30 keyScenes = pruneScenes(keyScenes); end % 重新构建主问题约束 rebuild_MP(); end end

4.3 我在实际调试中踩过的坑

这段是花了最多时间绕过的坑，写出来给各位省一点头发。

第一个坑：第二阶段变量里不小心混入了整数变量。子问题如果想通过对偶转换成单层 max，必须满足第二阶段是连续 LP。一旦里面出现 0-1 变量，强对偶条件不成立，整个子问题解法就会失效。我一开始为了建模方便，把储能充放电状态也放在第二阶段，导致子问题成为 MILP，C&CG 循环直接卡死。处理办法是把所有离散决策统一放到第一阶段，第二阶段只保留连续量。

第二个坑：不确定参数的维度方向要和约束矩阵对齐。光伏出力 u 一般是 T 维列向量，而主问题里的节点注入矩阵是按时刻排列的。如果用的是行向量和列向量混合，Yalmip 会默默帮你扩展成矩阵，但矩阵形状和你预期完全不同。这种错误很难发现，因为模型能正常求解，但结果明显不合理。

第三个坑：Gurobi 求解 MILP 主问题时，输出的 MIP 间隙默认是 1e-4，这在 C&CG 外循环里非常费时间。我在实际迭代中把主问题的 MIP Gap 单独放宽到 1e-2，外循环收敛间隙设成 1e-2，整体求解时间能下降 60%，而最终目标值差异只有百分之零点几。

第四个坑：场景距离度量。一开始我用的是 24 维时域曲线直接算欧氏距离，结果聚出来的“边界场景”经常是整体功率特别高的场景，导致关键场景集里全是高光伏、高负荷的大功率场景，对紧急情况没有辨识度。后来改成 PCA 降维后聚类，问题才解决。

5. 算例结果：关键场景数量、鲁棒性与经济性的权衡

场景辨别算法到底值不值得用？我用一个小型微网算例做了完整测试。

5.1 测试微网结构与参数

微网采用单母线结构，包含一台 200kW 柴油发电机，一台 300kW 光伏，一套 100kW/200kWh 储能，和一个峰值为 250kW 的负荷。调度周期为 24 小时，分辨率 1 小时。柴油机发电成本设为 0.8 元/kWh，光伏成本忽略，储能充放电效率 95%，购电价格按峰谷分时设置。

不确定性描述如下：光伏预测值取历史晴天典型曲线，实际出力允许在预测值的 70%~120% 之间波动；负荷允许在预测值的 90%~110% 之间波动；预算参数 Gamma 取 8，表示 24 小时内最多累计偏差不超过一定总量。

5.2 不同关键场景数量下的优化结果对比

我做了四组对照实验，关键场景集数量分别是 5、10、20、50，同时保留一个直接使用全部 500 个场景做两阶段随机规划的对照组作为计算基准参考。

这里有一个很关键的观察：从 5 个关键场景增加到 20 个时，成本和计算时间都在快速上升，但成本偏差从 +2.7% 缩小到 +0.4%。超过 20 个场景后，成本改善已经非常有限，但计算时间还在线性增长。所以在算例中，K=20 是一个“性价比转折点”。

为什么 5 个场景的偏差反而有 2.7%？因为 5 个关键场景太少，只能覆盖少数典型形状，无法体现不确定性集合边界的多样性。最恶劣场景一旦超出关键场景覆盖范围，鲁棒优化可能找不到安全解，于是只能通过更保守的调度来凑合，成本自然偏高。这正好说明了“关键场景辨别”不是越少越好，而是要达到覆盖度与规模的平衡。

5.3 C&CG收敛曲线与运行时间

以 20 个初始关键场景为例，C&CG 的迭代过程如下：

• 第 1 轮：主问题下界 LB=5810 元，子问题上界 UB=7050 元，间隙 17.6%。原因是第一阶段决策还比较乐观，子问题找出的恶劣场景给了一个很残酷的成本上界。
• 第 2 轮：恶劣场景加入主问题，LB 上升到 6280，UB 下降到 6820，间隙缩小到 7.9%。
• 第 3 轮：又有两个新场景被主问题吸收，LB 提高到 6600，UB 稳定在 6690，间隙 1.35%。
• 第 4 轮：间隙低于 1%，算法终止。

值得说明的是，后加入的“恶劣场景”并不是极端高负荷场景，而是“傍晚负荷高峰但光伏出力下滑、储能又已放空”的场景。这种场景在单纯蒙特卡洛采样的前 20 个关键场景里不容易出现，因为概率低；但 C&CG 反复迭代后，它能通过子问题识别出来并被动态加入。这正是两阶段鲁棒相比一阶段鲁棒的优势。

6. 这套方法可以继续往哪个方向扩展

写完这套实现之后，我又把目光投向了几个相关方向，个人觉得拓展空间很大。

第一个是配电网级的两阶段鲁棒重构。微网还可以作为一个“节点”进入配电网优化，不确定性除了分布式电源还有馈线故障概率。关键场景辨别算法在配电网场景数更多、维度更高的情况下，优势会更明显。前提是主问题的网络重构部分要建好线性化潮流模型。

第二个是把场景生成从纯统计模型换成深度学习模型。比如用 VAE 或 GAN 生成更逼真的光伏和负荷场景，再用本文的关键场景辨别算法做筛选。我在另一个小项目里试过用生成对抗网络模拟极端天气下的光伏出力，效果不错，但训练数据要有较长时间尺度的历史记录。

第三个是滚动时域在线调度。两阶段鲁棒优化原本偏日前，但如果在每个滚动窗口内都用关键场景识别算法快速刷新场景集，就能把它变成准在线算法。配合 MPC 思想，调度模型每 15 分钟更新一次，第一次阶段决策只执行第一个时段，后面重新优化。我在算例里测过，K=10 时单次滚动优化 3 秒以内就能完成，具备在线应用潜力。

第四个是扩展目标函数，把碳排放、设备磨损、电池寿命都放进去。两阶段鲁棒优化的内层如果是多目标，需要做一些妥协处理，但外层框架不用大改。

最后再分享一个小技巧：如果你在调试 C&CG 时发现上下界一直不收敛，先别怀疑算法，去检查子问题的可行域是不是被第一阶段决策“卡死”了。比如储能 SOC 初始值设得不好，可能导致第二阶段无论怎么调都无法满足功率平衡。这种情况本质是第一阶段没有保证第二阶段可行，需要额外加入“可行割”约束。我在第一次实现时就踩了这个坑，花了整整一天才发现是储能 SOC 初值的问题。把它修正后，迭代收敛就顺滑多了。