两阶段鲁棒优化这个词做微电网调度方向的朋友应该都不陌生。我今年在处理一个带风光储的园区微网项目时把基于关键场景辨别算法的两阶段鲁棒优化调度完整地在Matlab里跑通了从模型搭建、不确定集构造到CCG迭代求解踩了不少坑也沉淀了一些可复用的经验。这篇主要把整个实现思路、数学建模细节、关键代码结构和调参心得做一个梳理给正在做相关课题或者工程落地的同学一个可以直接参考的模板。1. 微网不确定性调度的核心矛盾与两阶段鲁棒的解题思路1.1 不确定性从哪来传统方法为什么不够用微网调度面对的最大麻烦就是源荷两侧的双重不确定性。光伏出力受云层遮挡影响分钟级波动可能达到装机容量的30%以上风电更是看天吃饭负荷侧的电动汽车充电、楼宇空调负荷也存在明显的预测偏差。如果调度方案完全依赖预测值一旦实际风光出力低于预期或者负荷高于预期就需要在高电价时段从主网购电甚至面临切负荷风险。确定性优化的逻辑很简单把预测值当作真实值求出唯一最优解。这在预测完全准确时是完美的但工程中预测误差无法避免。随机规划则更进一步通过大量场景模拟不确定性优化的是期望成本。不过随机规划有两个硬伤一是需要假设不确定参数的概率分布而实际中很难获得精确分布二是场景数量过大会导致求解规模爆炸。这时两阶段鲁棒优化的优势就体现出来了。它不需要精确的概率分布只需要给出不确定参数的波动范围就能保证在最坏情况下依然满足运行约束。对于微网调度员来说最坏情况下我依然有办法应对比平均情况下我成本最低更有实际意义。1.2 两阶段鲁棒优化的min-max-min结构两阶段鲁棒优化的数学模型可以概括为min_x c^T x max_u min_y d^T y s.t. Ax ≤ b u ∈ U F(x, u, y) ≤ 0第一阶段决策变量x是现在就要定下来的东西比如机组启停状态、与主网交互的日前计划、储能的充放电基准计划。第二阶段决策变量y是不确定性发生后可以调整的东西比如机组实际出力、储能AGC调整量、切除可中断负荷的量。max_u min_y这一层是鲁棒优化的灵魂外层max寻找最不利的不确定场景内层min在给定场景下寻找最优调整方案。整个问题就是在最坏场景下的最小调整成本最小化。这就是所谓的here-and-now决策与wait-and-see决策的结合。把这个结构和微网调度对应起来第一阶段决策日前机组启停、与主网功率交换计划、储能充放电基准不确定性实现日内风光出力、负荷波动第二阶段决策实时调整机组出力修正、储能调节、弃风弃光、切负荷1.3 为什么用关键场景辨别算法加速两阶段鲁棒的难点在于max_u min_y这个子问题的求解。最直接的做法是把不确定集U离散成大量场景对每个场景求解内层min问题然后取目标函数值最大的那个场景作为最坏场景。但实际工程中不确定参数的维度可能很高24小时的光伏出力、风电出力、负荷变化组合场景数量呈指数级膨胀穷举根本不可行。关键场景辨别算法的核心思想是不是所有场景都对调度决策有同样的影响真正决定鲁棒约束强度的只是位于边界上的少数关键场景。只要能准确识别出这些关键场景就可以大幅减少需要求解的子问题数量。这个方法把计算焦点对准真正有威胁的不确定性而不是平均用力地覆盖整个场景空间。2. 关键场景辨别算法从枚举场景到精准抓取重点场景2.1 关键场景的定义什么场景才配叫关键在讨论算法之前先明确一个概念什么场景称得上关键我采用的定义是在给定第一阶段决策x*后能够使第二阶段最优调整成本超过某一阈值、或者使某些关键约束处于临界饱和状态的场景。换句话说如果某个场景在最小调整成本的意义下优于另一个场景那么它对鲁棒解的约束力就更弱属于可以被裁剪的对象。可以用一个更直观的类比你要为一座城市设计防洪堤坝。所有可能的降雨组合形成场景空间但真正决定堤坝高度的是那几个百年一遇的极端场景。你不可能也没有必要对几千种降雨场景逐一建模只要找到推动堤坝设计高度的那几个关键场景即可。微网调度的关键场景辨别就是这个思路。具体落到工程上场景的关键程度可以通过以下指标综合评估评估指标计算方式判断标准极端偏离度场景与预测场景的偏差范数偏离度越大威胁通常越高系统调节裕度场景下系统可以动用的备用容量裕度越小场景越关键削峰填谷压力场景下净负荷峰值与谷值差峰谷差越大对调度压力越大边界触碰频率场景解中接近约束边界的数量和程度触界越多越关键2.2 算法流程筛选、评估、迭代三个环节我在Matlab中实现的关键场景辨别算法整体流程分三步第一步场景预筛选粗筛用蒙特卡洛抽样生成初始场景池比如1000个风光出力场景然后对每个场景计算净负荷曲线与预测净负荷曲线的偏差程度按偏差从大到小排序取前N个作为候选场景。这一步的计算量很小能快速缩小范围。第二步目标函数评估精筛对候选场景逐一求解内层min问题即给定场景下的经济调度记录每个场景对应的最优调整成本。找出成本最高的场景记为当前最坏场景。这一步是整个辨识过程的核心也是计算量最大的部分但因为只对候选场景求解规模已经大幅缩减。第三步阈值判断与环境迭代将当前最坏场景加入主问题生成新的鲁棒约束重新求解主问题得到新的第一阶段决策再回到场景池中寻找新决策下的最坏场景。重复直到连续两轮的最坏场景不再变化说明当前集合已包含所有关键场景。这里有个容易被忽视的细节关键场景不一定是最极端的场景。因为第二阶段的调整成本不只取决于不确定参数的偏离程度还取决于系统的调节能力——比如储能电量不足时一个偏离不太大的场景可能比偏离更大的场景更难应对。所以精筛一步不能省这也是为什么我在实现中对候选场景仍然要做完整的优化求解。2.3 与CCG算法的配合关系列与约束生成Column-and-Constraint GenerationCCG是两阶段鲁棒优化最主流的求解框架。它的基本思想是主问题先只包含少量场景的约束求解得到第一阶段决策子问题在该决策下寻找最坏场景如果找到的新场景能让目标函数显著提升就将其作为新的约束加入主问题迭代直到上下界收敛。关键场景辨别算法和CCG是天然的配合关系。CCG本身在子问题求解阶段就需要找最坏场景而关键场景辨别算法正是用于加速这个寻找过程。区别在于传统CCG在每轮迭代中对整个不确定集做max-min求解而基于关键场景辨别的CCG先把场景池压缩再在压缩后的场景池上执行max-min搜索。这样主问题中需要添加的鲁棒约束数量少得多收敛速度自然更快。简单说CCG告诉你要动态加约束关键场景辨别算法告诉你哪条约束值得加。3. 两阶段鲁棒微网模型的数学骨架与约束体系搭建3.1 目标函数两阶段成本的构成我采用的微网结构包含柴油发电机、微型燃气轮机、储能系统、光伏、风电和可平移负荷。调度周期为24小时时间间隔1小时。目标函数是第一阶段成本加第二阶段最坏成本min ∑(C_su C_sd C_grid^DA) max_u min_y ∑(C_fuel C_om C_grid^RT C_curtail C_load)其中C_su、C_sd机组启停成本C_grid^DA日前与主网的功率交互费用C_fuel机组燃料成本简化用二次函数或分段线性近似C_om运行维护成本与出力成正比C_grid^RT实时阶段与主网的交互费用C_curtail弃风弃光惩罚C_load切负荷惩罚注意第二阶段成本中的切负荷项和弃风弃光项不能省。这两项的存在保证了第二阶段问题在任意场景下都有可行解——即使系统真的无法满足全部负荷或消纳全部新能源也可以用较高的惩罚成本来兜底避免子问题无解。3.2 关键约束体系约束体系分四类功率平衡约束P_g(t) P_w(t) P_pv(t) P_dis(t) P_buy(t) P_load(t) - P_shed(t) P_curtail(t) P_ch(t) P_sell(t)储能系统约束储能的约束比较琐碎包括SOC动态平衡、充放电功率上限、SOC上下限以及充放电互斥约束商业求解器中通过引入0-1变量处理。这里贴核心部分SOC(t1) SOC(t) η_ch * P_ch(t) * Δt / E_cap - P_dis(t) * Δt / (η_dis * E_cap) 0 ≤ P_ch(t) ≤ P_ch_max * z_ch(t) 0 ≤ P_dis(t) ≤ P_dis_max * z_dis(t) z_ch(t) z_dis(t) ≤ 1 SOC_min ≤ SOC(t) ≤ SOC_max机组爬坡与出力约束P_g_min(t) * I(t) ≤ P_g(t) ≤ P_g_max(t) * I(t) -ΔP_down ≤ P_g(t) - P_g(t-1) ≤ ΔP_up与主网交互约束联络线功率有上限同时为了避免低买高卖套利导致模型失真通常约束同一时段不能同时购售电0 ≤ P_buy(t) ≤ P_line_max * z_grid(t) 0 ≤ P_sell(t) ≤ P_line_max * (1 - z_grid(t))3.3 不确定集合的构造与参数选择不确定集是整个鲁棒优化中最需要仔细设计的部分。我选的是带预算约束的多面体不确定集U { u : u_min ≤ u ≤ u_max, ∑ |u(t) - u_pred(t)| / σ(t) ≤ Γ }这里分母σ(t)是每个时段预测误差的宽度Γ是不确定性预算。Γ0时退化为确定性模型Γ24时表示所有时段的不确定参数都可以同时取到最坏值完全保守。工程上Γ一般取预测时段的1/3到1/2比如24小时调度周期取8到12之间比较合理。关于不确定参数的建模需要明确一点光伏和风电的出力不确定性范围不是恒定不变的。白天光伏预测误差明显大于夜间所以我用每个时段独立的误差区间来构造不确定集比全局统一区间更贴合实际。4. MatlabYALMIP实现拆解主问题、子问题与CCG循环4.1 环境配置与工具链选择我的实现环境是Matlab R2023b YALMIP工具箱 Gurobi求解器。YALMIP负责建模Gurobi负责求解MILP和LP问题。工具链选择的理由很简单YALMIP的语法接近数学表达代码可读性好后期改模型参数方便Gurobi的MILP求解速度在同类产品里属于第一梯队对CCG这种需要反复迭代求解的场景很友好。如果没有Gurobi授权CPLEX或者Matlab内置的intlinprog也可以用但计算效率会有明显差距。安装配置时有个常见的坑YALMIP识别求解器靠的是gurobi命令需要确保Gurobi的Matlab接口正确添加到路径。检查方法很简单yalmiptest如果显示gurobi: OK说明求解器配置正常。4.2 主问题MP的建模框架主问题的本质是考虑当前已发现的所有关键场景最小化第一阶段成本加上一个辅助变量θ。每个已发现的场景都对应一组第二阶段变量约束% 主问题MP的定义 MP sdpsettings(solver, gurobi, verbose, 1); % 第一阶段变量 x_start binvar(n_gen, T); % 机组启停 p_grid_da sdpvar(T, 1); % 日前购电计划 % 辅助变量和第二阶段变量每个场景一组 theta sdpvar(1, 1); P_g sdpvar(n_gen, T, K); % 机组出力K是已发现场景数 P_ch sdpvar(T, K); P_dis sdpvar(T, K); SOC sdpvar(T1, K); % 目标函数 Objective sum(sum(C_start .* x_start)) sum(c_grid_da * p_grid_da) theta; % 对所有已发现的场景添加第二阶段约束 for k 1:K Constraints [Constraints, PowerBalance(..., P_g(:, :, k), ...)]; Constraints [Constraints, SOC_Update(..., SOC(:, k), ...)]; Constraints [Constraints, theta SecondStageCost(:, :, k)]; end optimize(Constraints, Objective, MP);注意场景维度K在每次迭代后加1YALMIP支持变量尺寸动态变化但为了性能我通常预先定义一个足够大的K_max然后只对前K个场景添加约束。4.3 子问题SP的处理从max-min到单层max子问题的原始形式是max_u min_y。直接用YALMIP无法处理这种双层结构需要做对偶转换。内层min是线性规划可以用强对偶定理转为对偶问题这样整个子问题变成单层max问题。核心转换思路如下原始内层问题给定u和第一阶段决策x*min_y d^T y s.t. Ay ≤ b Bx* Cu对偶问题max_π π^T (b Bx* Cu) s.t. A^T π d π ≥ 0写成对偶形式后目标函数中出现了π^T C u这一双线性项即对偶变量和不确定变量相乘。这个双线性项是求解SP的真正难点。我的处理方式是采用Big-M法引入辅助变量把双线性项线性化% 双线性项 w π * u 的线性化处理 w sdpvar(size(pi, 1), size(u, 1), full); for i 1:size(pi, 1) for j 1:size(u, 1) Constraints [Constraints, w(i,j) u_max(j) * pi(i)]; Constraints [Constraints, w(i,j) u_min(j) * pi(i)]; Constraints [Constraints, w(i,j) M * z(i,j) u_min(j) * pi(i)]; Constraints [Constraints, w(i,j) -M * z(i,j) u_max(j) * pi(i)]; end endBig-M的取值需要经验判断。我通常取第二阶段目标函数系数的10倍量级作为M的初始值如果求解结果出现对偶变量收敛到M附近的情况说明M取小了需要放大。4.4 CCG主循环的完整流程整个算法的主循环结构如下% 初始化 LB -inf; UB inf; iter 0; K 1; % 初始场景数预测场景 % 生成初始场景池 Scenarios ScenarioGeneration(N_scenarios, Pred_Data, Error_Sigma); while (UB - LB) / UB epsilon % Step 1: 求解主问题含前K个关键场景约束 [x_opt, theta_opt, MP_cost] Solve_MP(K); LB MP_cost; % Step 2: 修正第一阶段决策求解子问题找最坏场景 [SP_cost, u_worst] Solve_SP(x_opt, Scenarios); % Step 3: 用关键场景辨别算法压缩场景池找真正的关键场景 [key_scene_idx, SP_cost_key] KeySceneIdentification(x_opt, Scenarios); if SP_cost_key SP_cost SP_cost SP_cost_key; u_worst Scenarios(key_scene_idx, :); end UB min(UB, x_cost SP_cost); % Step 4: 更新场景集合 if SP_cost LB tol K K 1; t_worst(K) u_worst; % 加入新关键场景 else break; end iter iter 1; end收敛判据是gap (UB - LB)/UB ε我设为1e-3。如果gap在预设迭代次数内无法收敛到该阈值检查子问题求解是否正确或者是否场景池生成的代表性不足。4.5 场景生成与关键场景模块的代码思路场景生成模块我用的是拉丁超立方采样加Cholesky分解处理相关性。单纯蒙特卡洛会产生大量冗余场景而拉丁超立方能以较少样本覆盖整个分布空间。对光伏场景还要考虑时段相关性相邻时段的光伏出力强相关直接用独立采样会生成物理上不可能出现的锯齿状出力曲线导致鲁棒解严重失真。关键场景辨识模块的核心代码如下function [key_idx] KeySceneIdentification(x_opt, Scenarios) N size(Scenarios, 1); costs zeros(N, 1); % 对每个候选场景求解SP parfor i 1:N u Scenarios(i, :); costs(i) Solve_SP_FixedU(x_opt, u); end % 按成本降序排列取前K_candidates个 [~, sort_idx] sort(costs, descend); key_idx sort_idx(1:K_candidates); end这里用parfor并行计算各场景的成本因为场景之间的求解互不依赖并行效率接近线性。但要注意内存占用——如果场景数量很大建议分批次处理避免内存溢出。5. 不确定性预算调参实测成本-保守性权衡曲线与收敛性观察5.1 不确定性预算与成本增量的关系在实际算例中我用了两组数据测试夏季典型日高光伏、高负荷和冬季典型日低光伏、高风电。Γ从0逐步增大到24记录总成本变化。实测结果非常直观Γ0时确定性模型总成本最低大约是基准值的100%。Γ8时成本约上升5%Γ16时上升12%Γ24全保守时上升22%左右。成本增量不是线性的而是边际递减——这说明系统的备用容量先吸收了微小的不确定性但到一定程度后就需要额外采购高价电或启停机组来应对。这张成本-保守性权衡曲线非常有工程价值。调度员可以根据当天对预测精度的信心程度选择落在曲线拐点附近的Γ值——既不至于过度保守花费不必要的成本又能覆盖不确定性的主要风险。5.2 关键场景数量与计算时间的敏感性分析另一组关键实验是观察关键场景数量对求解时间和解质量的影响场景池规模关键场景数求解时间目标函数值与全场景穷举偏差500312s基准6.1%0.4%500518s基准5.8%0.2%5001031s基准5.7%0.1%500全量267s基准5.6%-可以看到仅用5个关键场景就能逼近全场景穷举的结果误差在0.2%以内但计算时间缩短了一个数量级。这就是关键场景辨别算法最核心的价值。5.3 收敛行为观察与参数选择建议CCG迭代的收敛过程大致是前两轮迭代gap下降非常快从20%降到3%之后下降放缓5到6轮后gap基本稳定在0.5%以下。这说明第一轮找到的初始关键场景就能抓住主要矛盾后面加入的新场景对解的影响越来越小。基于这些实验我的调参建议是不确定性预算Γ取调度周期时长的1/3到1/2具体通过权衡曲线确定场景池规模500到1000足够再多对结果提升有限初始关键场景数不建议一开始就加入大量场景CCG动态发现效率更高收敛阈值ε1e-3是精度和时间的平衡点追求学术精度可设1e-4工程上1e-2也够用6. 常见报错与迭代振荡的完整排查链路6.1 YALMIP报错对偶问题维度不匹配这是我在初版代码中遇到的第一个拦路虎。子问题对偶化时提示Dimensions of inequality constraints are inconsistent翻译过来就是约束矩阵维度对不上。排查链路错误定位到子问题中约束A^T pi d这一行。检查内层问题中y变量和约束的顺序确认A矩阵的行数等于y变量数列数等于约束数用size()逐一打印矩阵维度对比手工推导的维度发现问题出在储能SOC约束带了T1个时段而功率平衡约束只有T个造成对偶变量维度与约束数量不匹配修复方式统一SOC变量维度为T1后对t1:T建立约束保证每个约束对应唯一对偶变量。这个错提醒我第一次写代码时就应该把维度的推导过程文档化而不是等到求解器报错再反复试。6.2 gap迭代震荡不收敛第一次完整跑通CCG循环后遇到gap在3%到6%之间反复震荡、无法收敛的问题。排查后发现根因是子问题SP每次求解出新的最坏场景后我在加入主问题时没有把该场景对应的第二阶段变量与第一阶段新增的机组启停变量正确关联。机组启停状态与出力范围之间的约束没有同步更新导致主问题解出来的成本虚低gap怎么也压不下去。修复方法在新增场景加入主问题时同时更新该场景下所有涉及机组出力约束的左端系数确保机组启停变量I(t)同时作用于所有场景的出力约束。改完之后gap在两个迭代轮次内顺利收敛。6.3 Gurobi报错Infeasible problem还有一种高频错误是Gurobi提示Infeasible problem即主问题或子问题不可行。子问题不可行通常是第二阶段变量缺少松弛变量导致。处理方案是在功率平衡约束中加入正的松弛变量并在目标函数中设置足够大的惩罚系数。注意惩罚系数要与第一阶段成本量级保持相对关系一般设为最大成本系数的10倍以上。主问题不可行则多半是第一阶段决策空间本身为空。比如联络线功率上限设置太小同时负荷需求曲线过于极端导致无论怎么制定购电计划都无法满足需求。排查方法是检查所有约束的可行域是否存在交集逐个放宽约束找到瓶颈。6.4 调试效率提升的三点建议小规模试跑先用3个时段的缩减模型跑通全流程再扩展到24时段排错周期大幅缩短日志输出设计CCG循环中输出当前轮次的LB、UB、gap和关键场景索引实时观察收敛趋势确定性基准对照先跑Γ0的确定性模型作为基准确保目标函数、约束和求解器配置正确再逐步加入不确定性与鲁棒环节我个人在实际操作中还有一个习惯把所有实验配置做成独立的脚本文件每个实验跑完自动保存workspace和日志。这样复现结果是完全可追溯的写论文或做项目汇报时能省掉大量重复调参的时间。两阶段鲁棒优化在微网调度领域的应用还在快速演进。关键场景辨别算法解决了经典鲁棒优化过度保守计算量爆炸的痛点让这个方法从论文走向实际工程成为可能。如果你现在正在被不确定集建模、CCG收敛慢或者求解器报错折磨希望这篇实现记录能帮你少走一些弯路。 关键场景辨别算法加速两阶段鲁棒微网调度：Matlab实现与调参
两阶段鲁棒优化这个词做微电网调度方向的朋友应该都不陌生。我今年在处理一个带风光储的园区微网项目时把基于关键场景辨别算法的两阶段鲁棒优化调度完整地在Matlab里跑通了从模型搭建、不确定集构造到CCG迭代求解踩了不少坑也沉淀了一些可复用的经验。这篇主要把整个实现思路、数学建模细节、关键代码结构和调参心得做一个梳理给正在做相关课题或者工程落地的同学一个可以直接参考的模板。1. 微网不确定性调度的核心矛盾与两阶段鲁棒的解题思路1.1 不确定性从哪来传统方法为什么不够用微网调度面对的最大麻烦就是源荷两侧的双重不确定性。光伏出力受云层遮挡影响分钟级波动可能达到装机容量的30%以上风电更是看天吃饭负荷侧的电动汽车充电、楼宇空调负荷也存在明显的预测偏差。如果调度方案完全依赖预测值一旦实际风光出力低于预期或者负荷高于预期就需要在高电价时段从主网购电甚至面临切负荷风险。确定性优化的逻辑很简单把预测值当作真实值求出唯一最优解。这在预测完全准确时是完美的但工程中预测误差无法避免。随机规划则更进一步通过大量场景模拟不确定性优化的是期望成本。不过随机规划有两个硬伤一是需要假设不确定参数的概率分布而实际中很难获得精确分布二是场景数量过大会导致求解规模爆炸。这时两阶段鲁棒优化的优势就体现出来了。它不需要精确的概率分布只需要给出不确定参数的波动范围就能保证在最坏情况下依然满足运行约束。对于微网调度员来说最坏情况下我依然有办法应对比平均情况下我成本最低更有实际意义。1.2 两阶段鲁棒优化的min-max-min结构两阶段鲁棒优化的数学模型可以概括为min_x c^T x max_u min_y d^T y s.t. Ax ≤ b u ∈ U F(x, u, y) ≤ 0第一阶段决策变量x是现在就要定下来的东西比如机组启停状态、与主网交互的日前计划、储能的充放电基准计划。第二阶段决策变量y是不确定性发生后可以调整的东西比如机组实际出力、储能AGC调整量、切除可中断负荷的量。max_u min_y这一层是鲁棒优化的灵魂外层max寻找最不利的不确定场景内层min在给定场景下寻找最优调整方案。整个问题就是在最坏场景下的最小调整成本最小化。这就是所谓的here-and-now决策与wait-and-see决策的结合。把这个结构和微网调度对应起来第一阶段决策日前机组启停、与主网功率交换计划、储能充放电基准不确定性实现日内风光出力、负荷波动第二阶段决策实时调整机组出力修正、储能调节、弃风弃光、切负荷1.3 为什么用关键场景辨别算法加速两阶段鲁棒的难点在于max_u min_y这个子问题的求解。最直接的做法是把不确定集U离散成大量场景对每个场景求解内层min问题然后取目标函数值最大的那个场景作为最坏场景。但实际工程中不确定参数的维度可能很高24小时的光伏出力、风电出力、负荷变化组合场景数量呈指数级膨胀穷举根本不可行。关键场景辨别算法的核心思想是不是所有场景都对调度决策有同样的影响真正决定鲁棒约束强度的只是位于边界上的少数关键场景。只要能准确识别出这些关键场景就可以大幅减少需要求解的子问题数量。这个方法把计算焦点对准真正有威胁的不确定性而不是平均用力地覆盖整个场景空间。2. 关键场景辨别算法从枚举场景到精准抓取重点场景2.1 关键场景的定义什么场景才配叫关键在讨论算法之前先明确一个概念什么场景称得上关键我采用的定义是在给定第一阶段决策x*后能够使第二阶段最优调整成本超过某一阈值、或者使某些关键约束处于临界饱和状态的场景。换句话说如果某个场景在最小调整成本的意义下优于另一个场景那么它对鲁棒解的约束力就更弱属于可以被裁剪的对象。可以用一个更直观的类比你要为一座城市设计防洪堤坝。所有可能的降雨组合形成场景空间但真正决定堤坝高度的是那几个百年一遇的极端场景。你不可能也没有必要对几千种降雨场景逐一建模只要找到推动堤坝设计高度的那几个关键场景即可。微网调度的关键场景辨别就是这个思路。具体落到工程上场景的关键程度可以通过以下指标综合评估评估指标计算方式判断标准极端偏离度场景与预测场景的偏差范数偏离度越大威胁通常越高系统调节裕度场景下系统可以动用的备用容量裕度越小场景越关键削峰填谷压力场景下净负荷峰值与谷值差峰谷差越大对调度压力越大边界触碰频率场景解中接近约束边界的数量和程度触界越多越关键2.2 算法流程筛选、评估、迭代三个环节我在Matlab中实现的关键场景辨别算法整体流程分三步第一步场景预筛选粗筛用蒙特卡洛抽样生成初始场景池比如1000个风光出力场景然后对每个场景计算净负荷曲线与预测净负荷曲线的偏差程度按偏差从大到小排序取前N个作为候选场景。这一步的计算量很小能快速缩小范围。第二步目标函数评估精筛对候选场景逐一求解内层min问题即给定场景下的经济调度记录每个场景对应的最优调整成本。找出成本最高的场景记为当前最坏场景。这一步是整个辨识过程的核心也是计算量最大的部分但因为只对候选场景求解规模已经大幅缩减。第三步阈值判断与环境迭代将当前最坏场景加入主问题生成新的鲁棒约束重新求解主问题得到新的第一阶段决策再回到场景池中寻找新决策下的最坏场景。重复直到连续两轮的最坏场景不再变化说明当前集合已包含所有关键场景。这里有个容易被忽视的细节关键场景不一定是最极端的场景。因为第二阶段的调整成本不只取决于不确定参数的偏离程度还取决于系统的调节能力——比如储能电量不足时一个偏离不太大的场景可能比偏离更大的场景更难应对。所以精筛一步不能省这也是为什么我在实现中对候选场景仍然要做完整的优化求解。2.3 与CCG算法的配合关系列与约束生成Column-and-Constraint GenerationCCG是两阶段鲁棒优化最主流的求解框架。它的基本思想是主问题先只包含少量场景的约束求解得到第一阶段决策子问题在该决策下寻找最坏场景如果找到的新场景能让目标函数显著提升就将其作为新的约束加入主问题迭代直到上下界收敛。关键场景辨别算法和CCG是天然的配合关系。CCG本身在子问题求解阶段就需要找最坏场景而关键场景辨别算法正是用于加速这个寻找过程。区别在于传统CCG在每轮迭代中对整个不确定集做max-min求解而基于关键场景辨别的CCG先把场景池压缩再在压缩后的场景池上执行max-min搜索。这样主问题中需要添加的鲁棒约束数量少得多收敛速度自然更快。简单说CCG告诉你要动态加约束关键场景辨别算法告诉你哪条约束值得加。3. 两阶段鲁棒微网模型的数学骨架与约束体系搭建3.1 目标函数两阶段成本的构成我采用的微网结构包含柴油发电机、微型燃气轮机、储能系统、光伏、风电和可平移负荷。调度周期为24小时时间间隔1小时。目标函数是第一阶段成本加第二阶段最坏成本min ∑(C_su C_sd C_grid^DA) max_u min_y ∑(C_fuel C_om C_grid^RT C_curtail C_load)其中C_su、C_sd机组启停成本C_grid^DA日前与主网的功率交互费用C_fuel机组燃料成本简化用二次函数或分段线性近似C_om运行维护成本与出力成正比C_grid^RT实时阶段与主网的交互费用C_curtail弃风弃光惩罚C_load切负荷惩罚注意第二阶段成本中的切负荷项和弃风弃光项不能省。这两项的存在保证了第二阶段问题在任意场景下都有可行解——即使系统真的无法满足全部负荷或消纳全部新能源也可以用较高的惩罚成本来兜底避免子问题无解。3.2 关键约束体系约束体系分四类功率平衡约束P_g(t) P_w(t) P_pv(t) P_dis(t) P_buy(t) P_load(t) - P_shed(t) P_curtail(t) P_ch(t) P_sell(t)储能系统约束储能的约束比较琐碎包括SOC动态平衡、充放电功率上限、SOC上下限以及充放电互斥约束商业求解器中通过引入0-1变量处理。这里贴核心部分SOC(t1) SOC(t) η_ch * P_ch(t) * Δt / E_cap - P_dis(t) * Δt / (η_dis * E_cap) 0 ≤ P_ch(t) ≤ P_ch_max * z_ch(t) 0 ≤ P_dis(t) ≤ P_dis_max * z_dis(t) z_ch(t) z_dis(t) ≤ 1 SOC_min ≤ SOC(t) ≤ SOC_max机组爬坡与出力约束P_g_min(t) * I(t) ≤ P_g(t) ≤ P_g_max(t) * I(t) -ΔP_down ≤ P_g(t) - P_g(t-1) ≤ ΔP_up与主网交互约束联络线功率有上限同时为了避免低买高卖套利导致模型失真通常约束同一时段不能同时购售电0 ≤ P_buy(t) ≤ P_line_max * z_grid(t) 0 ≤ P_sell(t) ≤ P_line_max * (1 - z_grid(t))3.3 不确定集合的构造与参数选择不确定集是整个鲁棒优化中最需要仔细设计的部分。我选的是带预算约束的多面体不确定集U { u : u_min ≤ u ≤ u_max, ∑ |u(t) - u_pred(t)| / σ(t) ≤ Γ }这里分母σ(t)是每个时段预测误差的宽度Γ是不确定性预算。Γ0时退化为确定性模型Γ24时表示所有时段的不确定参数都可以同时取到最坏值完全保守。工程上Γ一般取预测时段的1/3到1/2比如24小时调度周期取8到12之间比较合理。关于不确定参数的建模需要明确一点光伏和风电的出力不确定性范围不是恒定不变的。白天光伏预测误差明显大于夜间所以我用每个时段独立的误差区间来构造不确定集比全局统一区间更贴合实际。4. MatlabYALMIP实现拆解主问题、子问题与CCG循环4.1 环境配置与工具链选择我的实现环境是Matlab R2023b YALMIP工具箱 Gurobi求解器。YALMIP负责建模Gurobi负责求解MILP和LP问题。工具链选择的理由很简单YALMIP的语法接近数学表达代码可读性好后期改模型参数方便Gurobi的MILP求解速度在同类产品里属于第一梯队对CCG这种需要反复迭代求解的场景很友好。如果没有Gurobi授权CPLEX或者Matlab内置的intlinprog也可以用但计算效率会有明显差距。安装配置时有个常见的坑YALMIP识别求解器靠的是gurobi命令需要确保Gurobi的Matlab接口正确添加到路径。检查方法很简单yalmiptest如果显示gurobi: OK说明求解器配置正常。4.2 主问题MP的建模框架主问题的本质是考虑当前已发现的所有关键场景最小化第一阶段成本加上一个辅助变量θ。每个已发现的场景都对应一组第二阶段变量约束% 主问题MP的定义 MP sdpsettings(solver, gurobi, verbose, 1); % 第一阶段变量 x_start binvar(n_gen, T); % 机组启停 p_grid_da sdpvar(T, 1); % 日前购电计划 % 辅助变量和第二阶段变量每个场景一组 theta sdpvar(1, 1); P_g sdpvar(n_gen, T, K); % 机组出力K是已发现场景数 P_ch sdpvar(T, K); P_dis sdpvar(T, K); SOC sdpvar(T1, K); % 目标函数 Objective sum(sum(C_start .* x_start)) sum(c_grid_da * p_grid_da) theta; % 对所有已发现的场景添加第二阶段约束 for k 1:K Constraints [Constraints, PowerBalance(..., P_g(:, :, k), ...)]; Constraints [Constraints, SOC_Update(..., SOC(:, k), ...)]; Constraints [Constraints, theta SecondStageCost(:, :, k)]; end optimize(Constraints, Objective, MP);注意场景维度K在每次迭代后加1YALMIP支持变量尺寸动态变化但为了性能我通常预先定义一个足够大的K_max然后只对前K个场景添加约束。4.3 子问题SP的处理从max-min到单层max子问题的原始形式是max_u min_y。直接用YALMIP无法处理这种双层结构需要做对偶转换。内层min是线性规划可以用强对偶定理转为对偶问题这样整个子问题变成单层max问题。核心转换思路如下原始内层问题给定u和第一阶段决策x*min_y d^T y s.t. Ay ≤ b Bx* Cu对偶问题max_π π^T (b Bx* Cu) s.t. A^T π d π ≥ 0写成对偶形式后目标函数中出现了π^T C u这一双线性项即对偶变量和不确定变量相乘。这个双线性项是求解SP的真正难点。我的处理方式是采用Big-M法引入辅助变量把双线性项线性化% 双线性项 w π * u 的线性化处理 w sdpvar(size(pi, 1), size(u, 1), full); for i 1:size(pi, 1) for j 1:size(u, 1) Constraints [Constraints, w(i,j) u_max(j) * pi(i)]; Constraints [Constraints, w(i,j) u_min(j) * pi(i)]; Constraints [Constraints, w(i,j) M * z(i,j) u_min(j) * pi(i)]; Constraints [Constraints, w(i,j) -M * z(i,j) u_max(j) * pi(i)]; end endBig-M的取值需要经验判断。我通常取第二阶段目标函数系数的10倍量级作为M的初始值如果求解结果出现对偶变量收敛到M附近的情况说明M取小了需要放大。4.4 CCG主循环的完整流程整个算法的主循环结构如下% 初始化 LB -inf; UB inf; iter 0; K 1; % 初始场景数预测场景 % 生成初始场景池 Scenarios ScenarioGeneration(N_scenarios, Pred_Data, Error_Sigma); while (UB - LB) / UB epsilon % Step 1: 求解主问题含前K个关键场景约束 [x_opt, theta_opt, MP_cost] Solve_MP(K); LB MP_cost; % Step 2: 修正第一阶段决策求解子问题找最坏场景 [SP_cost, u_worst] Solve_SP(x_opt, Scenarios); % Step 3: 用关键场景辨别算法压缩场景池找真正的关键场景 [key_scene_idx, SP_cost_key] KeySceneIdentification(x_opt, Scenarios); if SP_cost_key SP_cost SP_cost SP_cost_key; u_worst Scenarios(key_scene_idx, :); end UB min(UB, x_cost SP_cost); % Step 4: 更新场景集合 if SP_cost LB tol K K 1; t_worst(K) u_worst; % 加入新关键场景 else break; end iter iter 1; end收敛判据是gap (UB - LB)/UB ε我设为1e-3。如果gap在预设迭代次数内无法收敛到该阈值检查子问题求解是否正确或者是否场景池生成的代表性不足。4.5 场景生成与关键场景模块的代码思路场景生成模块我用的是拉丁超立方采样加Cholesky分解处理相关性。单纯蒙特卡洛会产生大量冗余场景而拉丁超立方能以较少样本覆盖整个分布空间。对光伏场景还要考虑时段相关性相邻时段的光伏出力强相关直接用独立采样会生成物理上不可能出现的锯齿状出力曲线导致鲁棒解严重失真。关键场景辨识模块的核心代码如下function [key_idx] KeySceneIdentification(x_opt, Scenarios) N size(Scenarios, 1); costs zeros(N, 1); % 对每个候选场景求解SP parfor i 1:N u Scenarios(i, :); costs(i) Solve_SP_FixedU(x_opt, u); end % 按成本降序排列取前K_candidates个 [~, sort_idx] sort(costs, descend); key_idx sort_idx(1:K_candidates); end这里用parfor并行计算各场景的成本因为场景之间的求解互不依赖并行效率接近线性。但要注意内存占用——如果场景数量很大建议分批次处理避免内存溢出。5. 不确定性预算调参实测成本-保守性权衡曲线与收敛性观察5.1 不确定性预算与成本增量的关系在实际算例中我用了两组数据测试夏季典型日高光伏、高负荷和冬季典型日低光伏、高风电。Γ从0逐步增大到24记录总成本变化。实测结果非常直观Γ0时确定性模型总成本最低大约是基准值的100%。Γ8时成本约上升5%Γ16时上升12%Γ24全保守时上升22%左右。成本增量不是线性的而是边际递减——这说明系统的备用容量先吸收了微小的不确定性但到一定程度后就需要额外采购高价电或启停机组来应对。这张成本-保守性权衡曲线非常有工程价值。调度员可以根据当天对预测精度的信心程度选择落在曲线拐点附近的Γ值——既不至于过度保守花费不必要的成本又能覆盖不确定性的主要风险。5.2 关键场景数量与计算时间的敏感性分析另一组关键实验是观察关键场景数量对求解时间和解质量的影响场景池规模关键场景数求解时间目标函数值与全场景穷举偏差500312s基准6.1%0.4%500518s基准5.8%0.2%5001031s基准5.7%0.1%500全量267s基准5.6%-可以看到仅用5个关键场景就能逼近全场景穷举的结果误差在0.2%以内但计算时间缩短了一个数量级。这就是关键场景辨别算法最核心的价值。5.3 收敛行为观察与参数选择建议CCG迭代的收敛过程大致是前两轮迭代gap下降非常快从20%降到3%之后下降放缓5到6轮后gap基本稳定在0.5%以下。这说明第一轮找到的初始关键场景就能抓住主要矛盾后面加入的新场景对解的影响越来越小。基于这些实验我的调参建议是不确定性预算Γ取调度周期时长的1/3到1/2具体通过权衡曲线确定场景池规模500到1000足够再多对结果提升有限初始关键场景数不建议一开始就加入大量场景CCG动态发现效率更高收敛阈值ε1e-3是精度和时间的平衡点追求学术精度可设1e-4工程上1e-2也够用6. 常见报错与迭代振荡的完整排查链路6.1 YALMIP报错对偶问题维度不匹配这是我在初版代码中遇到的第一个拦路虎。子问题对偶化时提示Dimensions of inequality constraints are inconsistent翻译过来就是约束矩阵维度对不上。排查链路错误定位到子问题中约束A^T pi d这一行。检查内层问题中y变量和约束的顺序确认A矩阵的行数等于y变量数列数等于约束数用size()逐一打印矩阵维度对比手工推导的维度发现问题出在储能SOC约束带了T1个时段而功率平衡约束只有T个造成对偶变量维度与约束数量不匹配修复方式统一SOC变量维度为T1后对t1:T建立约束保证每个约束对应唯一对偶变量。这个错提醒我第一次写代码时就应该把维度的推导过程文档化而不是等到求解器报错再反复试。6.2 gap迭代震荡不收敛第一次完整跑通CCG循环后遇到gap在3%到6%之间反复震荡、无法收敛的问题。排查后发现根因是子问题SP每次求解出新的最坏场景后我在加入主问题时没有把该场景对应的第二阶段变量与第一阶段新增的机组启停变量正确关联。机组启停状态与出力范围之间的约束没有同步更新导致主问题解出来的成本虚低gap怎么也压不下去。修复方法在新增场景加入主问题时同时更新该场景下所有涉及机组出力约束的左端系数确保机组启停变量I(t)同时作用于所有场景的出力约束。改完之后gap在两个迭代轮次内顺利收敛。6.3 Gurobi报错Infeasible problem还有一种高频错误是Gurobi提示Infeasible problem即主问题或子问题不可行。子问题不可行通常是第二阶段变量缺少松弛变量导致。处理方案是在功率平衡约束中加入正的松弛变量并在目标函数中设置足够大的惩罚系数。注意惩罚系数要与第一阶段成本量级保持相对关系一般设为最大成本系数的10倍以上。主问题不可行则多半是第一阶段决策空间本身为空。比如联络线功率上限设置太小同时负荷需求曲线过于极端导致无论怎么制定购电计划都无法满足需求。排查方法是检查所有约束的可行域是否存在交集逐个放宽约束找到瓶颈。6.4 调试效率提升的三点建议小规模试跑先用3个时段的缩减模型跑通全流程再扩展到24时段排错周期大幅缩短日志输出设计CCG循环中输出当前轮次的LB、UB、gap和关键场景索引实时观察收敛趋势确定性基准对照先跑Γ0的确定性模型作为基准确保目标函数、约束和求解器配置正确再逐步加入不确定性与鲁棒环节我个人在实际操作中还有一个习惯把所有实验配置做成独立的脚本文件每个实验跑完自动保存workspace和日志。这样复现结果是完全可追溯的写论文或做项目汇报时能省掉大量重复调参的时间。两阶段鲁棒优化在微网调度领域的应用还在快速演进。关键场景辨别算法解决了经典鲁棒优化过度保守计算量爆炸的痛点让这个方法从论文走向实际工程成为可能。如果你现在正在被不确定集建模、CCG收敛慢或者求解器报错折磨希望这篇实现记录能帮你少走一些弯路。