两阶段鲁棒微网调度这套东西最近在圈子里面是真的火。标题里“基于关键场景辨别算法”这几个字我猜不少人是一眼扫过去就开始找代码了但真正自己动手把两阶段鲁棒和场景辨识揉到一起跑通一个Matlab算例才发现坑全在后头。这篇就围绕我复现这类模型时的完整思路和实际踩坑记录来写不讲空话只讲怎么从零搭起这套两阶段鲁棒微网优化调度框架。1. 整体设计与思路拆解为什么是“两阶段”加“鲁棒”1.1 微网调度为什么绕不开不确定性微网里的分布式光伏和风电出力天生就是波动的负荷也有随机性。传统确定性调度只取一个预测值来算结果往往是“看着最优实际一跑就崩”。我遇到过最典型的情况光伏预测出力拉满调度方案里外购电很少结果第二天早上云层一厚实际出力只有预测的60%缺的功率只能临时高价外购甚至切负荷。这就是确定性模型的脆弱性。所以在微网优化调度里“不确定性”不是理论上的概念而是每天都要面对的现实。处理不确定性的主流思路无非三种随机优化、鲁棒优化、分布鲁棒优化。随机优化需要知道不确定量的精确概率分布现实里很难拿到分布鲁棒虽然好但建模和求解复杂度高。对于大多数工程场景两阶段鲁棒优化其实是性价比最高的选择——它不依赖精确概率分布只需要知道不确定量的波动区间就能给出一个“在最坏情况下依然可行且经济性可接受”的调度方案。1.2 两阶段鲁棒的核心逻辑先决策后调整两阶段鲁棒优化的“两阶段”本质上是把决策变量分成两组。第一阶段是“现在必须定下来”的变量比如机组启停状态、与主网的购售电协议量、储能是否充电的粗计划第二阶段是“等不确定量揭晓后再调整”的变量比如实际出力下的储能充放电功率、机组实际出力、切负荷量等。用大白话讲第一阶段是“先定一个能保底的方案”第二阶段是“在这个保底方案下针对最恶劣的场景做最小代价的调整”。这就形成了一个三层结构——外层最小化费用、中层最大化不确定量找最恶劣场景、内层最小化调整费用。其中中层和内层嵌套着作为第一阶段方案的“最坏情况评估器”也就是所谓的“寻优-对抗”结构。我当初第一次看这个结构的时候绕了很久才想明白一件事鲁棒优化不是把所有场景都列出来而是只针对“让系统最难受的那个场景”做优化。只要能在这个场景下保证可行且经济那其他场景自然也能应付。这个“最恶劣场景”就是标题里说的“关键场景”的雏形。1.3 为什么需要“关键场景辨别”而不是枚举场景理论上如果光伏出力有10个可能取值区间负荷有10个可能取值区间那组合起来就是100个场景。如果再加上风电、电价不确定性场景数量会爆炸式增长。而且两阶段鲁棒里第二阶段是个优化问题每个场景都要算一遍枚举所有场景的算力成本完全不可接受。关键场景辨别算法的思路就是不枚举全部场景而是通过某种规则、指标或迭代策略只挑出那些对调度结果影响最大的几个场景来代表整个不确定集合。这跟机器学习里的“样本筛选”逻辑有点像——不是所有样本都有训练价值有些样本反而会误导模型挑出有代表性的关键样本才是重点。这里我补充说明一下标题里说的“关键场景辨别算法”在实际文献里有几种不同实现路子。一种是基于“风险指标排序”计算每个场景下系统的失负荷量或调整成本排序后取最大的几个另一种是在CCG列与约束生成迭代过程中自然产生的“最恶劣场景”——每次迭代解出来的那个场景就是当前最关键的场景。还有一种是基于场景相似度聚类的思路把大量场景聚成几类每类取一个代表场景。我们这次复现主要用的是CCG框架下“每次迭代主动辨识最恶劣场景”的思路这也是目前两阶段鲁棒里最主流、最稳定的做法。2. 核心细节解析不确定性集合模型与关键场景辨识2.1 不确定性集合怎么建盒式、椭球式还是多面体式两阶段鲁棒的建模第一步是定义不确定性集合。不同集合形状直接决定求解难度和保守程度。最常用的是盒式不确定集合也就是给每个不确定量定义一个上下界。比如光伏出力P_pv的实际值落在[预测值-偏差预测值偏差]区间内。盒式集合简单直观但有个问题它假设所有不确定量同时取到最恶劣值这在实际中几乎不可能发生所以结果偏保守。比盒式集合稍微精细一点的是带预算约束的盒式集合也就是“1-范数无穷范数”那种带预算的约束。这种集合限制“最多只能有Γ个不确定量同时达到极端值”能够有效控制保守度。我们在Matlab实现中采用的是带预算的盒式集合。这个选择的理由很实在椭球式集合需要转二阶锥约束处理起来烦多面体式集合虽然灵活但参数标定麻烦带预算的盒式集合只加了两个整数参数Γ和Γ_inf效果立竿见影而且能和CCG算法无缝配合。预算值Γ的物理意义要解释清楚Γ越大允许同时偏离预测值的不确定量就越多方案越保守Γ0时就退化成确定性模型。实际调试时我习惯把Γ从0开始一点点增大观察总成本的变化曲线——成本增幅突然变大的那个拐点往往就是“性价比最高”的预算值。2.2 关键场景辨别的实现路径CCG迭代中的子问题贡献CCG算法的核心思想是“主问题-子问题”交替迭代而关键场景辨识就藏在子问题求解的过程中。主问题MP是在当前已知的关键场景集合下求解第一阶段的决策变量和对应的最小成本。子问题SP是在给定第一阶段方案后寻找让系统调整成本最大的不确定量取值——说白了就是当前方案最大的“软肋”在哪里。如果找出来的这个最大成本大于主问题给出的成本就把这个场景对应的约束添加到主问题里继续迭代如果两者相等说明已经收敛。这个过程中“找最大成本对应场景”的这一步就是关键场景辨识。不需要枚举场景只需要在每一轮迭代中通过求解优化问题直接定位到“最恶劣”的那个不确定量组合。随着迭代进行关键场景会被逐个“挖”出来——通常迭代3到8轮就能收敛而枚举场景可能需要上万次计算。具体到子问题的求解两阶段鲁棒里最常用的处理办法是对第二阶段线性规划取其强对偶把内层的min转化为max这样中层和外层就合并成了一个单层max问题。如果是混合整数第二阶段问题处理起来会更复杂一些需要线性化处理或引入大M法。我们这里的模型第二阶段全是连续变量储能充放电、机组出力所以可以直接用强对偶转化这也是这类模型最常见的假设条件。2.3 场景辨识的数学表达与注意点子问题转化后的对偶问题会引入对偶变量乘子同时会出现双线性项不确定量与对偶变量的乘积。这个双线性项一般用大M法引入辅助变量来处理把不确定量的连续取值离散化到若干个取值点。这里有一个实操中很容易踩的坑大M值的选取。M太小可能把真正的最恶劣场景排除在外M太大会造成数值病态求解器容易报“numerical issues”。我这边调试下来比较稳妥的方式是根据不确定量的实际物理边界比如光伏出力不可能超过装机容量×效率把M设定为该边界值的2到3倍然后逐步微调。子问题解出来的不确定量取值就是本轮“辨别”出的关键场景。把它传给主问题在主问题中加入对应的第二阶段变量和约束然后继续迭代。这整个过程在Matlab里的实现并不复杂关键是逻辑要清晰主问题越“壮”约束越多鲁棒性越强子问题找的“茬”越准收敛越快。3. 实操过程与核心环节实现Matlab代码框架3.1 环境配置Yalmip Cplex/Gurobi Matlab这套代码不是纯Matlab就能跑的必须要装优化求解器。我的建议配置是Matlab 2020a及以上版本R2019b也能跑但部分语法会有小坑Yalmip工具箱求解器用Cplex或者Gurobi。装Yalmip很简单去官网下载压缩包解压后添加到Matlab路径即可。Cplex或Gurobi相对麻烦一些需要注册学术账号下载安装后要把对应路径添加到Matlab环境变量里。我建议装完先跑一个小测试代码% 测试求解器是否可用 x sdpvar(1,1); optimize([x 0, x 1], x, sdpsettings(solver,cplex));如果这能正常求解说明环境配置没问题。如果报错说找不到求解器多半是路径没加对或者Matlab版本和求解器版本不兼容。另外我强烈建议把求解器详细的输出关掉不然Cplex那满屏的迭代日志看得人脑壳疼。在sdpsettings里设置verbose, 0只保留我们自己打印的迭代信息既清爽又能看清楚CCG的收敛过程。3.2 微网系统结构与参数设置这里以一个典型的交流微网为例来说明参数设置思路。系统中包含一台燃气轮机MT、一台柴油机DE、一组储能BESS、光伏电站PV、本地负荷以及与配网的联络线。关键的参数我建议按下表来设置实际数值可以根据自己的算例调整参数值说明MT容量1.5 MW燃气轮机最大出力DE容量1.0 MW柴油机最大出力BESS容量2.0 MWh储能容量BESS功率0.5 MW最大充/放电功率PV装机1.2 MW光伏峰值联络线功率1.0 MW与主网交换功率上限负荷峰值2.0 MW典型日负荷峰值调度周期24h时间分辨率1h不确定量参数设置如下光伏出力和负荷预测误差分别取±15%和±10%预算值Γ_pv和Γ_load都设为4意思是全天24个时段里最多有4个时段允许光伏和负荷同时取到预测区间边界值。成本参数方面MT的燃料成本系数设为0.6DE的设为0.75单位是元/kWh。储能充放电成本设为0.1元/kWh用于计及损耗向主网购电价格采用分时电价峰时1.2元/kWh、平时0.8元/kWh、谷时0.4元/kWh。注意购电价是离散分档的这是个比较实用的细节跟实际市场机制能对上。3.3 CCG主循环的Matlab实现下面贴一段主循环的核心伪代码具体细节做了一些精简但整体框架是能跑的% 主程序两阶段鲁棒CCG迭代 %% 初始化 MP_cost inf; % 主问题最优值 SP_cost 0; % 子问题最优值 UB inf; % 上界主问题提供 LB -inf; % 下界子问题提供 iter 0; max_iter 20; tol 1e-4; % 初始场景取预测值场景作为第一个关键场景 U_scenarios nominal_scenario; %% CCG主循环 while true iter iter 1; % ---------- 主问题求解 ---------- [MP_cost, x1] solve_MP(U_scenarios); % 主问题给的解必然可行是上界 UB min(UB, MP_cost); % ---------- 子问题求解 ---------- [SP_cost, u_new] solve_SP(x1, uncertainty_set); % 子问题是对抗性的在当前方案下找最恶劣场景 % 如果最恶劣场景下的调整成本大于当前UB说明方案还不够鲁棒 if SP_cost UB tol break; % 收敛了关键场景已经全部被考虑 else % 把这个新场景加入场景集合 U_scenarios [U_scenarios, u_new]; % 下界更新 LB max(LB, SP_cost); end if iter max_iter warning(达到最大迭代次数可能未完全收敛); break; end end % 输出结果 fprintf(CCG迭代次数: %d\n, iter); fprintf(最优成本: %.4f 元\n, MP_cost);注意第七行初始场景我取的是预测值场景也就是确定性场景。这个选择不是随便定的——从预测值场景出发CCG能够保证至少第一轮就能产生一个有意义的方案后续迭代在这个基础上逐步增强鲁棒性整体收敛会比较平稳。3.4 子问题求解强对偶与线性化上面代码里solve_SP这个函数是核心中的核心。它的任务是给定第一阶段的机组启停和储能计划求在不确定集合内让“运行调整成本”最大的场景。这一步难在如何处理第二阶段优化问题的内层min。第二阶段变量的形式如下第一阶段给出的燃气轮机启停状态二进制决定了第二阶段燃气轮机的出力可行域。在第二阶段我们要最小化的是“偏离计划调度后的惩罚成本”——包括储能调整成本、切负荷惩罚、弃光惩罚、与主网的实时交易不平衡惩罚等。子问题写成数学形式后把所有第二阶段约束取对偶得到一个以不确定量为自变量的最大化问题。由于第二阶段的目标函数和约束都是线性的对偶问题依旧线性。对偶问题的目标里会出现不确定量乘对偶变量的双线性项这个地方我专门说明一下处理方式。function [obj, u_out] solve_SP(x1, unc) % 输入第一阶段决策 x1不确定集合 unc % 输出子问题目标值 obj最恶劣场景 u_out % 定义不确定变量 u_pv sdpvar(24, 1); u_load sdpvar(24, 1); % 不确定集合约束带预算 Constraints []; Constraints [Constraints, -abs_dev_pv u_pv abs_dev_pv]; Constraints [Constraints, -abs_dev_load u_load abs_dev_load]; Constraints [Constraints, sum(abs(u_pv)) Gamma_pv]; Constraints [Constraints, sum(abs(u_load)) Gamma_load]; % 第二阶段变量在u给定时的调整量 y sdpvar(...); % 储能充放电、MT调整量、购电调整量等 % 对偶问题的目标函数会包含双线性项 % 使用大M法线性化 M 10; % 需要根据实际参数调整 % ... % 求解这个 max 问题 options sdpsettings(solver, cplex, verbose, 0); optimize(Constraints, -obj_expr, options); % Yalmip默认min取负号转为max obj value(obj_expr); u_out [value(u_pv), value(u_load)]; end这里要重点提醒一个细节Yalmip默认求解的是min问题要把max问题转成min需要在目标函数前面加负号。很多人第一次写这里会把符号搞反结果子问题一直在找“最不恶劣”的场景整个迭代完全跑偏。这类双线性项的线性化是代码里最烦人的部分。不过如果用的是Cplex或Gurobi的新版本支持直接在目标里处理部分二次项可以把双线性项直接写成u * lambda的形式交给求解器两个大厂求解器都能直接求解非凸二次问题Cplex的QP和Gurobi的bilinear。当然为了兼容性和稳定性我还是推荐自己动手做大M线性化尤其当不确定集合里含有绝对值约束时。3.5 主问题求解主问题的形式相对简单function [cost, x1] solve_MP(U_scenarios) % U_scenarios 是已经积累的关键场景集合矩阵 % 每一列是一个场景或者一个场景一个时序向量 % 第一阶段变量 u_MT binvar(24, 1); % MT启停 u_DE binvar(24, 1); % DE启停 p_MT sdpvar(24, 1); % MT出力 p_DE sdpvar(24, 1); % DE出力 p_BESS_ch sdpvar(24, 1); % 储能充电 p_BESS_dis sdpvar(24, 1); % 储能放电 p_grid sdpvar(24, 1); % 购售电正购负售 Constraints []; % 对所有已知场景添加约束 for k 1:size(U_scenarios, 2) u_pv U_scenarios(1:24, k); % 第k个场景的光伏出力 u_load U_scenarios(25:48, k); % 第k个场景的负荷 % 功率平衡约束每个场景都必须满足 Constraints [Constraints, ... p_MT p_DE p_BESS_dis - p_BESS_ch p_grid u_pv ... u_load]; % 储能动态约束每个场景都要满足 Constraints [Constraints, ... SOC(:, k1) SOC(:, k) eta_ch * p_BESS_ch - p_BESS_dis / eta_dis]; % ... 其他约束按需添加 end % 目标函数包含第一阶段成本 所有场景下的第二阶段成本期望/最坏情况成本 obj sum(price_fix .* u_MT c_MT .* p_MT ... price_fix_DE .* u_DE c_DE .* p_DE ... c_grid .* p_grid c_BESS .* (p_BESS_ch p_BESS_dis)); % 求解 options sdpsettings(solver, cplex, verbose, 0); optimize(Constraints, obj, options); cost value(obj); x1 value([u_MT, u_DE, p_MT, p_DE, p_BESS_ch, p_BESS_dis, p_grid]); end这段代码里隐藏了一个很关键的设计决策主问题中所有第一阶段变量对所有场景共享但每个场景有自己独立的第二阶段变量SOC等。这意味着储能系统在每个场景下会有不同的运行轨迹但在第一阶段决策启停状态、容量配置上保持一致。这在CCG中是标准的处理方式也是保证收敛正确性的前提。如果你在复现时发现主问题和子问题来回震荡不收敛建议优先检查主问题中是否为每个场景都正确添加了独立的储能SOC变量和功率平衡约束——很多复现代码的错误就是这里把不同场景的SOC混在了一个约束里。4. 常见问题与排查技巧实录4.1 不收敛或收敛极慢这是两阶段鲁棒复现中最常遇到的问题。我之前调试一个类似模型时CCG迭代了30多次还在振荡最后排查出来是因为子问题里的预算约束写错了绝对值求和忘了加abs()导致不确定集合实际上无界子问题每次都往无穷大跑。还有一个常见原因是主问题目标函数没有把子问题的“最坏情况成本”包含进去。正确的做法是主问题目标应当等于第一阶段成本加上所有已辨识场景的第二阶段成本的线性组合每个场景一个权重一般是等权或者取最坏值。如果只加了第一阶段成本主问题就会疯狂压低第一阶段成本然后子问题每次都找一个巨大的调整成本导致上下界差距一直很大。排查技巧把每次迭代的主问题成本UB和子问题最大成本SP分别打印出来。正常情况下UB应该是单减的主问题约束越加越多可行域越来越小成本不会上升而SP应该是大概率渐增或波动的。如果UB在增加说明主问题约束写错了可能是把不同场景的约束耦合错了。4.2 求解器报数值问题或“Infeasible”表现为Cplex报infeasible或者Gurobi报Numerical trouble。check这几处常见原因排查要点解决办法大M值不匹配M设置过小导致本质可行解被错误剪掉逐步增大M观察解是否稳定量纲不一致功率单位kW、MW混用成本单位元/万元混用统一量纲建议全用p.u.或全用kW、元储能SOC上下界与功率约束矛盾SOC最小上限和最大下限之间没有可行区间检查是否满足SOC_min ch_energy SOC_max 之类的耦合约束对偶转化中遗漏约束第二阶段的某些约束没有取对偶导致子问题不可行仔细核对强对偶条件补充所有约束的对偶数值问题的经验Cplex对数字精度极其敏感我建议在sdpsettings里加上cplex.barrier.tol, 1e-7这类参数或者直接用Gurobi的gurobi.NumericFocus, 1来增强数值鲁棒性。但治本的办法是统一量纲——不要一股脑用MW也别一股脑用kW混合使用最容易出事。4.3 结果比确定性模型贵太多保守度过高如果算出来的总成本比确定性模型高了30%以上大概率是保守度设置得不对或者不确定性集合定义得过宽。我的调试经验是先用Γ1跑一遍看看成本和方案跟确定性相比变化多少然后逐渐增大Γ观察成本和储能充放电策略的变化。如果某个Γ值下成本突变剧烈很可能是方案从“平时充电、峰时放电”变成了“为了防止最坏情况而全天保留裕度”——这种方案虽然鲁棒但经济性很差实际运营中往往不可取。保守度、预算值的经济学解读很有意思预算值Γ其实就是“运维人员的风险偏好”。你愿意为最坏情况多准备多少预算本质上是在为“风险规避”定价。这个认识在实际工程中非常有用——拿同一套代码跑几个不同Γ的方案交给运营人员选比直接给定一个方案要专业得多。4.4 Yalmip函数维度不匹配Matlab复现这类代码时sdpvar的维度声明是个高频错误来源。尤其是repmat、kron这两个函数在构建多时段约束时如果维度不匹配Yalmip有时不会立刻报错而是会悄悄生成错误的约束矩阵最终导致结果完全不符合物理逻辑。我的建议是每定义一个关键约束后立即用size(Constraints)检查维度。如果约束矩阵的行数不是期望的24假设24时段马上回溯是哪一步出了问题。不要等整个模型构建完再调试那时候排查的复杂度会成倍增加。5. 从复现到改进这个模型还能怎么延展5.1 分布鲁棒优化DRO扩展两阶段鲁棒最被人诟病的一点是结果偏保守。如果你算出来的方案在实际运营中“太浪费”可以考虑升级为分布鲁棒优化。核心思路是不需要精确概率分布只需要知道不确定量的一阶矩和二阶矩信息均值和协方差然后构建一个包含所有可能分布的模糊集在最坏分布下做优化。分布鲁棒的Matlab实现并不比两阶段鲁棒复杂太多关键在于模糊集的建模方式矩模糊集或Wasserstein球模糊集第二阶段的处理方式跟CCG类似也是强对偶线性化。这个方向作为后续扩展相当顺滑。5.2 多微网互联与需求响应联动现在的微网很少单打独斗多微网互联、微网与配网的互动调度是热门方向。两阶段鲁棒框架天生适合处理“多个微网各自有不确定性但共享联络线容量”的场景——第一阶段决定各自机组启停第二阶段在互联约束下做联合调整关键场景辨别会从“单微网最恶劣场景”变成“多微网组合最恶劣场景”计算复杂度会显著上升但算法框架不变。5.3 考虑储能寿命衰减的调度策略储能电池的充放电循环会带来寿命衰减这在两阶段鲁棒里是一个很容易被忽略、但实际中非常重要的问题。我见过不少复现代码把储能当成“永动机”来调度——每一轮都满充满放算出来的方案看着很美好但实际储能两年就报废了。比较实用的改进方案是在目标函数中加入储能SOC偏移惩罚项让储能尽量避免长期保持在满充或全放的状态或者在第二阶段约束中加入循环次数上限。这两种方式都能在一定程度上反映寿命因素且不破坏两阶段鲁棒的模型结构。5.4 并行计算加速多场景评估如果后续把模型扩展到多微网或者更长时间尺度比如168小时的周调度CCG中每轮子问题求解的时间会显著增加。一个可行的加速方案是把子问题按场景拆分对每个时段或每个子系统的子问题用parfor并行求解。Matlab的并行计算工具箱在这里非常好用我试过在8核机器上跑能把子问题求解时间压到串行的四分之一左右。实操中我自己的一点体会最后说一个个人经验层面的东西。两阶段鲁棒微网调度这套模型论文里看着高深真正落地时难的不是数学推导而是对物理模型细节的把握和对求解器特性的理解。我建议第一次复现的人先跑一个确定性模型作为基准把储能SOC变化、机组出力等曲线画出来确认每个物理过程都合理再逐步引入不确定性集合和CCG迭代。这样即使后面迭代出现问题也能快速定位到是哪一层约束导致的不一致。另外一个小技巧调试CCG时可以把每一轮辨识出的关键场景画出来看看这些场景的时序特征。通常第一轮是“光伏最低负荷最高”的极端场景第二轮可能是“光伏最高负荷最低”的逆调峰场景后续几轮则是介于两者之间的过渡场景。看到这些场景的形状基本就能判断迭代是否正确——如果关键场景看起来全都是同一种形态多半是约束写错了导致搜索空间被错误压缩了。我自己的习惯是在代码里加上一段把每轮迭代的关键场景写入mat文件的逻辑跑完之后用load把这些场景导出来画图。这个对排查问题的帮助比打印一堆数字大得多。 基于关键场景辨别的两阶段鲁棒微网调度Matlab实现
两阶段鲁棒微网调度这套东西最近在圈子里面是真的火。标题里“基于关键场景辨别算法”这几个字我猜不少人是一眼扫过去就开始找代码了但真正自己动手把两阶段鲁棒和场景辨识揉到一起跑通一个Matlab算例才发现坑全在后头。这篇就围绕我复现这类模型时的完整思路和实际踩坑记录来写不讲空话只讲怎么从零搭起这套两阶段鲁棒微网优化调度框架。1. 整体设计与思路拆解为什么是“两阶段”加“鲁棒”1.1 微网调度为什么绕不开不确定性微网里的分布式光伏和风电出力天生就是波动的负荷也有随机性。传统确定性调度只取一个预测值来算结果往往是“看着最优实际一跑就崩”。我遇到过最典型的情况光伏预测出力拉满调度方案里外购电很少结果第二天早上云层一厚实际出力只有预测的60%缺的功率只能临时高价外购甚至切负荷。这就是确定性模型的脆弱性。所以在微网优化调度里“不确定性”不是理论上的概念而是每天都要面对的现实。处理不确定性的主流思路无非三种随机优化、鲁棒优化、分布鲁棒优化。随机优化需要知道不确定量的精确概率分布现实里很难拿到分布鲁棒虽然好但建模和求解复杂度高。对于大多数工程场景两阶段鲁棒优化其实是性价比最高的选择——它不依赖精确概率分布只需要知道不确定量的波动区间就能给出一个“在最坏情况下依然可行且经济性可接受”的调度方案。1.2 两阶段鲁棒的核心逻辑先决策后调整两阶段鲁棒优化的“两阶段”本质上是把决策变量分成两组。第一阶段是“现在必须定下来”的变量比如机组启停状态、与主网的购售电协议量、储能是否充电的粗计划第二阶段是“等不确定量揭晓后再调整”的变量比如实际出力下的储能充放电功率、机组实际出力、切负荷量等。用大白话讲第一阶段是“先定一个能保底的方案”第二阶段是“在这个保底方案下针对最恶劣的场景做最小代价的调整”。这就形成了一个三层结构——外层最小化费用、中层最大化不确定量找最恶劣场景、内层最小化调整费用。其中中层和内层嵌套着作为第一阶段方案的“最坏情况评估器”也就是所谓的“寻优-对抗”结构。我当初第一次看这个结构的时候绕了很久才想明白一件事鲁棒优化不是把所有场景都列出来而是只针对“让系统最难受的那个场景”做优化。只要能在这个场景下保证可行且经济那其他场景自然也能应付。这个“最恶劣场景”就是标题里说的“关键场景”的雏形。1.3 为什么需要“关键场景辨别”而不是枚举场景理论上如果光伏出力有10个可能取值区间负荷有10个可能取值区间那组合起来就是100个场景。如果再加上风电、电价不确定性场景数量会爆炸式增长。而且两阶段鲁棒里第二阶段是个优化问题每个场景都要算一遍枚举所有场景的算力成本完全不可接受。关键场景辨别算法的思路就是不枚举全部场景而是通过某种规则、指标或迭代策略只挑出那些对调度结果影响最大的几个场景来代表整个不确定集合。这跟机器学习里的“样本筛选”逻辑有点像——不是所有样本都有训练价值有些样本反而会误导模型挑出有代表性的关键样本才是重点。这里我补充说明一下标题里说的“关键场景辨别算法”在实际文献里有几种不同实现路子。一种是基于“风险指标排序”计算每个场景下系统的失负荷量或调整成本排序后取最大的几个另一种是在CCG列与约束生成迭代过程中自然产生的“最恶劣场景”——每次迭代解出来的那个场景就是当前最关键的场景。还有一种是基于场景相似度聚类的思路把大量场景聚成几类每类取一个代表场景。我们这次复现主要用的是CCG框架下“每次迭代主动辨识最恶劣场景”的思路这也是目前两阶段鲁棒里最主流、最稳定的做法。2. 核心细节解析不确定性集合模型与关键场景辨识2.1 不确定性集合怎么建盒式、椭球式还是多面体式两阶段鲁棒的建模第一步是定义不确定性集合。不同集合形状直接决定求解难度和保守程度。最常用的是盒式不确定集合也就是给每个不确定量定义一个上下界。比如光伏出力P_pv的实际值落在[预测值-偏差预测值偏差]区间内。盒式集合简单直观但有个问题它假设所有不确定量同时取到最恶劣值这在实际中几乎不可能发生所以结果偏保守。比盒式集合稍微精细一点的是带预算约束的盒式集合也就是“1-范数无穷范数”那种带预算的约束。这种集合限制“最多只能有Γ个不确定量同时达到极端值”能够有效控制保守度。我们在Matlab实现中采用的是带预算的盒式集合。这个选择的理由很实在椭球式集合需要转二阶锥约束处理起来烦多面体式集合虽然灵活但参数标定麻烦带预算的盒式集合只加了两个整数参数Γ和Γ_inf效果立竿见影而且能和CCG算法无缝配合。预算值Γ的物理意义要解释清楚Γ越大允许同时偏离预测值的不确定量就越多方案越保守Γ0时就退化成确定性模型。实际调试时我习惯把Γ从0开始一点点增大观察总成本的变化曲线——成本增幅突然变大的那个拐点往往就是“性价比最高”的预算值。2.2 关键场景辨别的实现路径CCG迭代中的子问题贡献CCG算法的核心思想是“主问题-子问题”交替迭代而关键场景辨识就藏在子问题求解的过程中。主问题MP是在当前已知的关键场景集合下求解第一阶段的决策变量和对应的最小成本。子问题SP是在给定第一阶段方案后寻找让系统调整成本最大的不确定量取值——说白了就是当前方案最大的“软肋”在哪里。如果找出来的这个最大成本大于主问题给出的成本就把这个场景对应的约束添加到主问题里继续迭代如果两者相等说明已经收敛。这个过程中“找最大成本对应场景”的这一步就是关键场景辨识。不需要枚举场景只需要在每一轮迭代中通过求解优化问题直接定位到“最恶劣”的那个不确定量组合。随着迭代进行关键场景会被逐个“挖”出来——通常迭代3到8轮就能收敛而枚举场景可能需要上万次计算。具体到子问题的求解两阶段鲁棒里最常用的处理办法是对第二阶段线性规划取其强对偶把内层的min转化为max这样中层和外层就合并成了一个单层max问题。如果是混合整数第二阶段问题处理起来会更复杂一些需要线性化处理或引入大M法。我们这里的模型第二阶段全是连续变量储能充放电、机组出力所以可以直接用强对偶转化这也是这类模型最常见的假设条件。2.3 场景辨识的数学表达与注意点子问题转化后的对偶问题会引入对偶变量乘子同时会出现双线性项不确定量与对偶变量的乘积。这个双线性项一般用大M法引入辅助变量来处理把不确定量的连续取值离散化到若干个取值点。这里有一个实操中很容易踩的坑大M值的选取。M太小可能把真正的最恶劣场景排除在外M太大会造成数值病态求解器容易报“numerical issues”。我这边调试下来比较稳妥的方式是根据不确定量的实际物理边界比如光伏出力不可能超过装机容量×效率把M设定为该边界值的2到3倍然后逐步微调。子问题解出来的不确定量取值就是本轮“辨别”出的关键场景。把它传给主问题在主问题中加入对应的第二阶段变量和约束然后继续迭代。这整个过程在Matlab里的实现并不复杂关键是逻辑要清晰主问题越“壮”约束越多鲁棒性越强子问题找的“茬”越准收敛越快。3. 实操过程与核心环节实现Matlab代码框架3.1 环境配置Yalmip Cplex/Gurobi Matlab这套代码不是纯Matlab就能跑的必须要装优化求解器。我的建议配置是Matlab 2020a及以上版本R2019b也能跑但部分语法会有小坑Yalmip工具箱求解器用Cplex或者Gurobi。装Yalmip很简单去官网下载压缩包解压后添加到Matlab路径即可。Cplex或Gurobi相对麻烦一些需要注册学术账号下载安装后要把对应路径添加到Matlab环境变量里。我建议装完先跑一个小测试代码% 测试求解器是否可用 x sdpvar(1,1); optimize([x 0, x 1], x, sdpsettings(solver,cplex));如果这能正常求解说明环境配置没问题。如果报错说找不到求解器多半是路径没加对或者Matlab版本和求解器版本不兼容。另外我强烈建议把求解器详细的输出关掉不然Cplex那满屏的迭代日志看得人脑壳疼。在sdpsettings里设置verbose, 0只保留我们自己打印的迭代信息既清爽又能看清楚CCG的收敛过程。3.2 微网系统结构与参数设置这里以一个典型的交流微网为例来说明参数设置思路。系统中包含一台燃气轮机MT、一台柴油机DE、一组储能BESS、光伏电站PV、本地负荷以及与配网的联络线。关键的参数我建议按下表来设置实际数值可以根据自己的算例调整参数值说明MT容量1.5 MW燃气轮机最大出力DE容量1.0 MW柴油机最大出力BESS容量2.0 MWh储能容量BESS功率0.5 MW最大充/放电功率PV装机1.2 MW光伏峰值联络线功率1.0 MW与主网交换功率上限负荷峰值2.0 MW典型日负荷峰值调度周期24h时间分辨率1h不确定量参数设置如下光伏出力和负荷预测误差分别取±15%和±10%预算值Γ_pv和Γ_load都设为4意思是全天24个时段里最多有4个时段允许光伏和负荷同时取到预测区间边界值。成本参数方面MT的燃料成本系数设为0.6DE的设为0.75单位是元/kWh。储能充放电成本设为0.1元/kWh用于计及损耗向主网购电价格采用分时电价峰时1.2元/kWh、平时0.8元/kWh、谷时0.4元/kWh。注意购电价是离散分档的这是个比较实用的细节跟实际市场机制能对上。3.3 CCG主循环的Matlab实现下面贴一段主循环的核心伪代码具体细节做了一些精简但整体框架是能跑的% 主程序两阶段鲁棒CCG迭代 %% 初始化 MP_cost inf; % 主问题最优值 SP_cost 0; % 子问题最优值 UB inf; % 上界主问题提供 LB -inf; % 下界子问题提供 iter 0; max_iter 20; tol 1e-4; % 初始场景取预测值场景作为第一个关键场景 U_scenarios nominal_scenario; %% CCG主循环 while true iter iter 1; % ---------- 主问题求解 ---------- [MP_cost, x1] solve_MP(U_scenarios); % 主问题给的解必然可行是上界 UB min(UB, MP_cost); % ---------- 子问题求解 ---------- [SP_cost, u_new] solve_SP(x1, uncertainty_set); % 子问题是对抗性的在当前方案下找最恶劣场景 % 如果最恶劣场景下的调整成本大于当前UB说明方案还不够鲁棒 if SP_cost UB tol break; % 收敛了关键场景已经全部被考虑 else % 把这个新场景加入场景集合 U_scenarios [U_scenarios, u_new]; % 下界更新 LB max(LB, SP_cost); end if iter max_iter warning(达到最大迭代次数可能未完全收敛); break; end end % 输出结果 fprintf(CCG迭代次数: %d\n, iter); fprintf(最优成本: %.4f 元\n, MP_cost);注意第七行初始场景我取的是预测值场景也就是确定性场景。这个选择不是随便定的——从预测值场景出发CCG能够保证至少第一轮就能产生一个有意义的方案后续迭代在这个基础上逐步增强鲁棒性整体收敛会比较平稳。3.4 子问题求解强对偶与线性化上面代码里solve_SP这个函数是核心中的核心。它的任务是给定第一阶段的机组启停和储能计划求在不确定集合内让“运行调整成本”最大的场景。这一步难在如何处理第二阶段优化问题的内层min。第二阶段变量的形式如下第一阶段给出的燃气轮机启停状态二进制决定了第二阶段燃气轮机的出力可行域。在第二阶段我们要最小化的是“偏离计划调度后的惩罚成本”——包括储能调整成本、切负荷惩罚、弃光惩罚、与主网的实时交易不平衡惩罚等。子问题写成数学形式后把所有第二阶段约束取对偶得到一个以不确定量为自变量的最大化问题。由于第二阶段的目标函数和约束都是线性的对偶问题依旧线性。对偶问题的目标里会出现不确定量乘对偶变量的双线性项这个地方我专门说明一下处理方式。function [obj, u_out] solve_SP(x1, unc) % 输入第一阶段决策 x1不确定集合 unc % 输出子问题目标值 obj最恶劣场景 u_out % 定义不确定变量 u_pv sdpvar(24, 1); u_load sdpvar(24, 1); % 不确定集合约束带预算 Constraints []; Constraints [Constraints, -abs_dev_pv u_pv abs_dev_pv]; Constraints [Constraints, -abs_dev_load u_load abs_dev_load]; Constraints [Constraints, sum(abs(u_pv)) Gamma_pv]; Constraints [Constraints, sum(abs(u_load)) Gamma_load]; % 第二阶段变量在u给定时的调整量 y sdpvar(...); % 储能充放电、MT调整量、购电调整量等 % 对偶问题的目标函数会包含双线性项 % 使用大M法线性化 M 10; % 需要根据实际参数调整 % ... % 求解这个 max 问题 options sdpsettings(solver, cplex, verbose, 0); optimize(Constraints, -obj_expr, options); % Yalmip默认min取负号转为max obj value(obj_expr); u_out [value(u_pv), value(u_load)]; end这里要重点提醒一个细节Yalmip默认求解的是min问题要把max问题转成min需要在目标函数前面加负号。很多人第一次写这里会把符号搞反结果子问题一直在找“最不恶劣”的场景整个迭代完全跑偏。这类双线性项的线性化是代码里最烦人的部分。不过如果用的是Cplex或Gurobi的新版本支持直接在目标里处理部分二次项可以把双线性项直接写成u * lambda的形式交给求解器两个大厂求解器都能直接求解非凸二次问题Cplex的QP和Gurobi的bilinear。当然为了兼容性和稳定性我还是推荐自己动手做大M线性化尤其当不确定集合里含有绝对值约束时。3.5 主问题求解主问题的形式相对简单function [cost, x1] solve_MP(U_scenarios) % U_scenarios 是已经积累的关键场景集合矩阵 % 每一列是一个场景或者一个场景一个时序向量 % 第一阶段变量 u_MT binvar(24, 1); % MT启停 u_DE binvar(24, 1); % DE启停 p_MT sdpvar(24, 1); % MT出力 p_DE sdpvar(24, 1); % DE出力 p_BESS_ch sdpvar(24, 1); % 储能充电 p_BESS_dis sdpvar(24, 1); % 储能放电 p_grid sdpvar(24, 1); % 购售电正购负售 Constraints []; % 对所有已知场景添加约束 for k 1:size(U_scenarios, 2) u_pv U_scenarios(1:24, k); % 第k个场景的光伏出力 u_load U_scenarios(25:48, k); % 第k个场景的负荷 % 功率平衡约束每个场景都必须满足 Constraints [Constraints, ... p_MT p_DE p_BESS_dis - p_BESS_ch p_grid u_pv ... u_load]; % 储能动态约束每个场景都要满足 Constraints [Constraints, ... SOC(:, k1) SOC(:, k) eta_ch * p_BESS_ch - p_BESS_dis / eta_dis]; % ... 其他约束按需添加 end % 目标函数包含第一阶段成本 所有场景下的第二阶段成本期望/最坏情况成本 obj sum(price_fix .* u_MT c_MT .* p_MT ... price_fix_DE .* u_DE c_DE .* p_DE ... c_grid .* p_grid c_BESS .* (p_BESS_ch p_BESS_dis)); % 求解 options sdpsettings(solver, cplex, verbose, 0); optimize(Constraints, obj, options); cost value(obj); x1 value([u_MT, u_DE, p_MT, p_DE, p_BESS_ch, p_BESS_dis, p_grid]); end这段代码里隐藏了一个很关键的设计决策主问题中所有第一阶段变量对所有场景共享但每个场景有自己独立的第二阶段变量SOC等。这意味着储能系统在每个场景下会有不同的运行轨迹但在第一阶段决策启停状态、容量配置上保持一致。这在CCG中是标准的处理方式也是保证收敛正确性的前提。如果你在复现时发现主问题和子问题来回震荡不收敛建议优先检查主问题中是否为每个场景都正确添加了独立的储能SOC变量和功率平衡约束——很多复现代码的错误就是这里把不同场景的SOC混在了一个约束里。4. 常见问题与排查技巧实录4.1 不收敛或收敛极慢这是两阶段鲁棒复现中最常遇到的问题。我之前调试一个类似模型时CCG迭代了30多次还在振荡最后排查出来是因为子问题里的预算约束写错了绝对值求和忘了加abs()导致不确定集合实际上无界子问题每次都往无穷大跑。还有一个常见原因是主问题目标函数没有把子问题的“最坏情况成本”包含进去。正确的做法是主问题目标应当等于第一阶段成本加上所有已辨识场景的第二阶段成本的线性组合每个场景一个权重一般是等权或者取最坏值。如果只加了第一阶段成本主问题就会疯狂压低第一阶段成本然后子问题每次都找一个巨大的调整成本导致上下界差距一直很大。排查技巧把每次迭代的主问题成本UB和子问题最大成本SP分别打印出来。正常情况下UB应该是单减的主问题约束越加越多可行域越来越小成本不会上升而SP应该是大概率渐增或波动的。如果UB在增加说明主问题约束写错了可能是把不同场景的约束耦合错了。4.2 求解器报数值问题或“Infeasible”表现为Cplex报infeasible或者Gurobi报Numerical trouble。check这几处常见原因排查要点解决办法大M值不匹配M设置过小导致本质可行解被错误剪掉逐步增大M观察解是否稳定量纲不一致功率单位kW、MW混用成本单位元/万元混用统一量纲建议全用p.u.或全用kW、元储能SOC上下界与功率约束矛盾SOC最小上限和最大下限之间没有可行区间检查是否满足SOC_min ch_energy SOC_max 之类的耦合约束对偶转化中遗漏约束第二阶段的某些约束没有取对偶导致子问题不可行仔细核对强对偶条件补充所有约束的对偶数值问题的经验Cplex对数字精度极其敏感我建议在sdpsettings里加上cplex.barrier.tol, 1e-7这类参数或者直接用Gurobi的gurobi.NumericFocus, 1来增强数值鲁棒性。但治本的办法是统一量纲——不要一股脑用MW也别一股脑用kW混合使用最容易出事。4.3 结果比确定性模型贵太多保守度过高如果算出来的总成本比确定性模型高了30%以上大概率是保守度设置得不对或者不确定性集合定义得过宽。我的调试经验是先用Γ1跑一遍看看成本和方案跟确定性相比变化多少然后逐渐增大Γ观察成本和储能充放电策略的变化。如果某个Γ值下成本突变剧烈很可能是方案从“平时充电、峰时放电”变成了“为了防止最坏情况而全天保留裕度”——这种方案虽然鲁棒但经济性很差实际运营中往往不可取。保守度、预算值的经济学解读很有意思预算值Γ其实就是“运维人员的风险偏好”。你愿意为最坏情况多准备多少预算本质上是在为“风险规避”定价。这个认识在实际工程中非常有用——拿同一套代码跑几个不同Γ的方案交给运营人员选比直接给定一个方案要专业得多。4.4 Yalmip函数维度不匹配Matlab复现这类代码时sdpvar的维度声明是个高频错误来源。尤其是repmat、kron这两个函数在构建多时段约束时如果维度不匹配Yalmip有时不会立刻报错而是会悄悄生成错误的约束矩阵最终导致结果完全不符合物理逻辑。我的建议是每定义一个关键约束后立即用size(Constraints)检查维度。如果约束矩阵的行数不是期望的24假设24时段马上回溯是哪一步出了问题。不要等整个模型构建完再调试那时候排查的复杂度会成倍增加。5. 从复现到改进这个模型还能怎么延展5.1 分布鲁棒优化DRO扩展两阶段鲁棒最被人诟病的一点是结果偏保守。如果你算出来的方案在实际运营中“太浪费”可以考虑升级为分布鲁棒优化。核心思路是不需要精确概率分布只需要知道不确定量的一阶矩和二阶矩信息均值和协方差然后构建一个包含所有可能分布的模糊集在最坏分布下做优化。分布鲁棒的Matlab实现并不比两阶段鲁棒复杂太多关键在于模糊集的建模方式矩模糊集或Wasserstein球模糊集第二阶段的处理方式跟CCG类似也是强对偶线性化。这个方向作为后续扩展相当顺滑。5.2 多微网互联与需求响应联动现在的微网很少单打独斗多微网互联、微网与配网的互动调度是热门方向。两阶段鲁棒框架天生适合处理“多个微网各自有不确定性但共享联络线容量”的场景——第一阶段决定各自机组启停第二阶段在互联约束下做联合调整关键场景辨别会从“单微网最恶劣场景”变成“多微网组合最恶劣场景”计算复杂度会显著上升但算法框架不变。5.3 考虑储能寿命衰减的调度策略储能电池的充放电循环会带来寿命衰减这在两阶段鲁棒里是一个很容易被忽略、但实际中非常重要的问题。我见过不少复现代码把储能当成“永动机”来调度——每一轮都满充满放算出来的方案看着很美好但实际储能两年就报废了。比较实用的改进方案是在目标函数中加入储能SOC偏移惩罚项让储能尽量避免长期保持在满充或全放的状态或者在第二阶段约束中加入循环次数上限。这两种方式都能在一定程度上反映寿命因素且不破坏两阶段鲁棒的模型结构。5.4 并行计算加速多场景评估如果后续把模型扩展到多微网或者更长时间尺度比如168小时的周调度CCG中每轮子问题求解的时间会显著增加。一个可行的加速方案是把子问题按场景拆分对每个时段或每个子系统的子问题用parfor并行求解。Matlab的并行计算工具箱在这里非常好用我试过在8核机器上跑能把子问题求解时间压到串行的四分之一左右。实操中我自己的一点体会最后说一个个人经验层面的东西。两阶段鲁棒微网调度这套模型论文里看着高深真正落地时难的不是数学推导而是对物理模型细节的把握和对求解器特性的理解。我建议第一次复现的人先跑一个确定性模型作为基准把储能SOC变化、机组出力等曲线画出来确认每个物理过程都合理再逐步引入不确定性集合和CCG迭代。这样即使后面迭代出现问题也能快速定位到是哪一层约束导致的不一致。另外一个小技巧调试CCG时可以把每一轮辨识出的关键场景画出来看看这些场景的时序特征。通常第一轮是“光伏最低负荷最高”的极端场景第二轮可能是“光伏最高负荷最低”的逆调峰场景后续几轮则是介于两者之间的过渡场景。看到这些场景的形状基本就能判断迭代是否正确——如果关键场景看起来全都是同一种形态多半是约束写错了导致搜索空间被错误压缩了。我自己的习惯是在代码里加上一段把每轮迭代的关键场景写入mat文件的逻辑跑完之后用load把这些场景导出来画图。这个对排查问题的帮助比打印一堆数字大得多。