两阶段鲁棒微网调度这套东西,最近在圈子里面是真的火。标题里“基于关键场景辨别算法”这几个字,我猜不少人是一眼扫过去就开始找代码了,但真正自己动手把两阶段鲁棒和场景辨识揉到一起,跑通一个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.6,DE的设为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(:, k+1) == 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把这些场景导出来画图。这个对排查问题的帮助,比打印一堆数字大得多。