☰
两阶段鲁棒微网优化调度:关键场景辨别算法与MATLAB实现
2026/9/29 16:22:37 网站建设 项目流程

这两年做微电网优化调度的朋友,应该都对“不确定性”又爱又恨。风电光伏出力飘忽不定,负荷预测总有偏差,你要是把备用留少了,最坏场景一来就可能越限甚至失负荷;留多了,运行成本又压不下来。我最近把一套基于关键场景辨别算法的两阶段鲁棒微网优化调度MATLAB代码完整跑通了,这里把思路、建模、代码结构和踩坑记录都整理出来,给正在做两阶段鲁棒、不确定优化或者微网调度的同学一个可复现的参考。这套代码的核心特点是:不依赖风电/光伏/负荷的概率分布,只要给预测区间即可;调度策略偏保守但能保证在最坏场景下系统不越限、不失负荷;在求解层面用关键场景辨别算法减少了传统C&CG(列与约束生成)框架下子问题的求解次数,实测收敛速度和稳定性都明显好于逐顶点枚举的做法。

1. 两阶段鲁棒微网调度到底在优化什么

1.1 先从“不确定性”这个痛点说起

微网调度相比传统电网调度,最大的麻烦是源荷两头都不确定。分布式风电出力取决于风速,光伏出力取决于辐照度,负荷也有很强的随机性。如果你只用一个确定性的预测值去做日前计划,第二天实际运行的时候大概率会出现偏差。偏差小时可以用储能、联络线功率修正,偏差大时就得切负荷或者弃风弃光,经济性和可靠性都不好看。

解决这个问题的思路无非两条:随机优化和鲁棒优化。随机优化需要知道不确定量的概率分布,然后用蒙特卡洛抽样或者场景树把随机问题离散化。优点是解出的是期望意义上的最优,缺点是概率分布往往拿不准,而且样本数量一大,求解规模会成倍膨胀。鲁棒优化换了一个角度,它不关心不确定量具体取多少,只关心“最坏情况下系统能不能扛得住”。只要给每个不确定量一个上下界,它就能在这个范围内把你需要防御的最恶劣场景找出来。微网调度里,风电光伏的区间预测比精确分布好做得多,所以两阶段鲁棒在微网场景下特别合适。

1.2 两阶段,指的是哪两个阶段

这里的两阶段是决策时序上的两阶段,不是物理上的两个调度时段。第一阶段(here-and-now)做的是要在不确定量实现之前就必须定下来的决策,比如机组启停、联络线日前购售电计划、储能是否参与某时段的充放电等。这些决策一般带整数变量,定了之后不能轻易改,否则会产生额外的启停成本或违约风险。

第二阶段(wait-and-see)是等风电、光伏、负荷的实际出力都“揭晓”之后,在第一阶段决策给定的前提下,做最省钱的运行调整。比如调整火电机组出力、储能实际充放电功率、与主网的实时交互功率,目的就是让总运行成本最低,同时满足功率平衡、爬坡、储能SOC等约束。把两阶段写在一起,就成了一个 min-max-min 结构:第一阶段最小化启停成本和第二阶段最坏情况下的运行成本;中间的 max 是在不确定集合里挑那个让运行成本最高的场景;最内层的 min 是在这个场景下做最优经济调度。代码里所有复杂的东西,本质上都是在处理这个 min-max-min。

1.3 谁适合参考这套代码

如果你正在做微电网、主动配电网、综合能源系统的优化调度;如果你在写两阶段鲁棒相关的论文但一直卡在子问题求解;或者你只懂确定性优化,想快速上手不确定优化,这套Matlab代码都值得对照着改一版。代码框架是通用的,换算例、换目标函数、加新能源类型都不会太费劲。

2. 关键场景辨别算法:把无穷场景压成有限几个

2.1 不确定集怎么建模

鲁棒优化的第一步是定义不确定集。代码里用的是一般的盒式不确定集加预算约束:风电出力偏差、光伏出力偏差、负荷偏差各自落在预测值加减一个δ的区间内,同时限制“同时偏离预测值的时间段/变量个数”不超过预算参数Γ。这个Γ很关键,Γ=0时结果等于确定性优化,Γ越大保护程度越高,解也越保守。Γ取一个适中的值,可以在保证鲁棒性的同时不牺牲太多经济性。

以风电为例,某个时段风机出力的不确定集合可以写成:

[ u \in U=\left{ u \ \middle| \ u^{\text{pre}} - \delta \le u \le u^{\text{pre}} + \delta, \ \sum |u - u^{\text{pre}}| \le \Gamma \right} ]

实际编程时,我会把各时段的风电、光伏、负荷偏差向量拼成一个大的不确定向量,统一建集合。这样做的好处是主问题和子问题的矩阵结构一致,后面加割平面的时候不用反复改索引。

2.2 为什么不能直接枚举所有顶点

两阶段鲁棒的子问题是一个 max-min 问题。理论上有有限顶点定理:如果不确定集是有界多面体,且内层问题对u是线性的,那么外层max的最优解一定可以取在不确定集的某个顶点上。很多人第一反应是“那就把所有顶点组合枚举一遍呗”。但组合数是2的N次方。比如一台风机24个时段,顶点组合就是2的24次方种,还不算光伏、负荷。逐顶点枚举在PPT上讲得通,工程上根本跑不动。

即便你只取每小时最恶劣的单一场景,也会漏掉多个时段联合恶劣的情况。风电在凌晨大发、光伏在中午大发、负荷在晚高峰飙升,这些在不同的时段各自取到极值,组合起来才是真正的最坏场景。如果只按单个变量独立最坏来处理,得到的结果太乐观,起不到鲁棒保护的作用。

2.3 场景辨别算法的核心逻辑

关键场景辨别算法做的事情,说白了就是:不枚举所有顶点,而是利用对偶乘子的符号和大小,直接定位“哪个不确定变量取哪个方向的极值,会让子问题目标函数增长最快”。具体到我实现的版本,思路分三步:

第一步,把内层min问题用对偶方法转成max问题,得到一个新的max-max结构,目标函数里出现对偶变量与不确定量的乘积项(λ^T B u),这是唯一的双线性项。

第二步,由于在最优解处u可以取区间端点,判定u取上界还是下界,取决于乘积项的系数符号。于是给每个不确定变量算一个“系数绝对值排名”,按排名分配有限的偏离预算。说白了就是:预算只有那么多,谁对目标函数影响最大,谁优先用掉这个偏离额度。

第三步,把识别出的关键场景作为新的割平面加入主问题,进入下一轮C&CG迭代。这样一来,每次迭代只需要识别一个(或少数几个)关键场景,而不是对所有顶点做一遍子问题,计算量会成数量级下降。这也是标题里“关键场景辨别算法”的实际作用。

我补充说明一下,这套思路是基于常见实现的逻辑补全。如果你拿到的代码用的是KKT条件或者对偶变换加大M线性化,本质上等价,只是在某些细节上(比如整数变量的处理)会有差异。

3. Matlab代码结构与实现细节

3.1 整体框架:别再把所有代码塞进一个脚本

我先说一个经验:这种两阶段鲁棒代码,千万不要把所有逻辑写在一个大脚本里。C&CG本身就有主问题、子问题、割平面更新三层循环,中间还有数据读入和结果绘图,全堆在一起,调试的时候你会疯掉的。我推荐按职责拆分成下面这几个文件:

main_ccg.m 主循环,调度主程序 init_parameters.m 全局参数、基础数据读入 build_uncertainty.m 不确定集与关键场景预计算 solve_master_problem.m 求解第一阶段+割平面主问题 solve_subproblem.m 求解子问题并提取关键场景 plot_results.m 画图(收敛曲线、调度结果)

main_ccg.m 的流程其实很固定:初始化上下界LB=-inf、UB=inf;先用预测场景解一次确定性模型,得到初始主问题解;进入C&CG循环;求解主问题更新LB;固定主问题第一阶段解,调用solve_subproblem.m求解最坏场景,得到UB;如果UB和LB的间隙小于收敛阈值(比如1e-3),退出;否则把子问题识别出的关键场景对应的第二阶段约束和变量加入主问题,继续循环。

整个过程可以用下面这段伪代码表达:

LB = -inf; UB = inf; k = 0; while abs(UB - LB) / abs(UB) > 1e-3 k = k + 1; % 求解主问题 [x_first, lb_k] = solve_master_problem(k); LB = max(LB, lb_k); % 固定第一阶段解,求解子问题 [u_k, ub_k] = solve_subproblem(x_first); UB = min(UB, ub_k); if abs(UB - LB) / abs(UB) < 1e-3 break; end % 把关键场景u_k加入主问题割平面 add_cutting_plane(k, u_k); end

3.2 主问题怎么建模:每轮割平面新增一组变量

主问题每一轮迭代都会增加一组“与关键场景相关的第二阶段变量和约束”。这样做的好处是主问题的规模虽然会增大,但始终只包含那些被识别出来的坏场景,而不是把所有场景一次性放进去。用YALMIP建模时,第一阶段变量里包含机组启停状态、联络线功率等,成本项用指示函数或二进制变量处理。

我习惯的做法是先把参数定义清楚,然后一次性写约束,避免在循环里反复拼接YALMIP表达式导致内存爆炸。主问题核心结构大致是这样:

% 主问题示例(示意) y_opt = binvar(nUnits, T, 'full'); % 机组启停 x_opt = sdpvar(nBuses, T, 'full'); % 第一阶段连续决策 objective = sum(sum(c_start .* y_opt)) ... % 启停成本 + sum(sum(rho_a .* p_import)) ... % 购电成本 + sum(sum(rho_b .* p_export)); % 售电收益 Constraints = [Constraints, A * y_opt + B * x_opt <= b_ub]; % 每轮关键场景加入后,新增一组第二阶段变量和约束 for k = 1:K y2{k} = sdpvar(nUnits, T, 'full'); x2{k} = sdpvar(nBuses, T, 'full'); % 关键场景下的功率平衡、爬坡、SOC约束 Constraints = [Constraints, C * y2{k} + D * x2{k} <= b_k(k)]; objective = objective + sum(sum(c_run .* y2{k})); end optimize(Constraints, objective, sdpsettings('solver','gurobi','verbose',0));

注意这只是伪代码示意,因为不同算例的系数矩阵差异很大,直接套用会出错。重点是理解“每一轮割平面只加一组新的第二阶段变量和约束”这个思想。每轮循环里新变量要放在元胞数组里,而不是覆盖同名变量,否则YALMIP会把上一轮的变量认错,导致约束串线。

3.3 子问题怎么处理对偶和双线性项

子问题是整个代码里最容易写错的地方。给定第一阶段解之后,原问题是一个 max-min 问题。我的实现方案是把内层min写成标准线性规划形式:

[ \min_{y} \ c^T y \quad \text{s.t.} \quad Ay \le b + Bu ]

如果第二阶段问题是连续线性规划,强对偶成立,就可以把它转成:

[ \max_{\lambda,u} \ \lambda^T(b + Bu) \quad \text{s.t.} \quad A^T\lambda \le c, \ \lambda \ge 0, \ u \in U ]

这样就把 max-min 变成了单层 max 问题,目标函数里只留下 λ^T B u 这个双线性项。处理这个双线性项有两条路:一是交给求解器自动线性化,但Gurobi处理双线性加整数变量的混合整数规划往往很慢;二是用关键场景辨别方法手动处理u——对每个u_i,看对应乘积项的系数符号决定取上界还是下界,有预算约束时,按系数绝对值排序来分配偏离额度。我强烈推荐第二种,因为快,而且逻辑完全可控。

代码里我用了一个小函数来生成关键场景:

function us = critical_scene(lambda_coef, delta, Gamma) % lambda_coef: 每个不确定变量对应的乘积项系数 % delta: 允许的最大偏差(正) % Gamma: 预算参数,即最多允许多少个变量同时取极值 n = length(delta); [~, idx] = sort(abs(lambda_coef), 'descend'); us = zeros(size(delta)); for i = 1:min(n, Gamma) pos = idx(i); us(pos) = delta(pos) * sign(lambda_coef(pos)); end end

这段代码看着简单,但它是整个加速的关键。同样一轮C&CG,如果你用枚举顶点的方式,可能要解2的24次方个子问题;用这个辨别函数,只需要解1个子问题,然后根据对偶乘子直接拼出关键场景。我实测下来,单轮子问题求解时间从“分钟级”降到了“秒级”,而且随着不确定变量维数增加,这个差距会越来越大。

3.4 收敛判据和停机条件

两阶段鲁棒的上下界更新逻辑要特别注意。主问题解出来的是下界LB,因为主问题只包含了已经识别出的关键场景,真实的最坏场景可能还没出现;子问题解出来的是上界UB,因为它是在给定第一阶段解之后求出的最坏运行成本,这个成本在完整的可行域里是可实现的。

我一般设收敛间隙为1e-3,如果设为1e-6会白白多迭代很多轮。另一个心得是:第一轮子问题不要从随机场景开始,直接用预测场景先解一轮,把初始上下界拉近,能明显减少迭代次数。这个方法在文献里叫“warm start”,实际效果比想象中明显。

4. 参数设置、算例分析与调参心得

4.1 一个可复现的小算例

为了验证代码正确性,建议先跑一个小算例。我这里用的是一个6节点微网:2台柴油机组、1台风电机组、1座光伏电站、1套储能装置,通过联络线与上级电网连接。分时购售电价,储能容量按微网峰值负荷的20%配置。

参数名称数值
风机额定容量200 kW
光伏额定容量150 kW
柴油机组单台容量100 kW
储能容量200 kWh
储能最大充放电功率50 kW
储能效率0.95
联络线功率上限300 kW
购电价峰值/谷值1.2 / 0.4 元/kWh
售电价峰值/谷值0.9 / 0.3 元/kWh
风电预测偏差比例15%
光伏预测偏差比例10%
负荷预测偏差比例5%

预测曲线我直接用了一组典型日数据:风电夜间大、白天小,光伏中午尖峰,负荷早晚两个高峰。这样的算例能同时考验机组启停、储能充放电和联络线功率三个环节的协调能力。

4.2 确定性调度和两阶段鲁棒调度的结果对比

跑出来的结果很有意思。确定性模型在预测场景下的日运行成本大约是1200元,看起来很美。但如果把确定性解放到最坏场景下校核,系统会出现失负荷:晚高峰时段风电和光伏同时低于预测,柴油机组已经满发也补不上缺口,联络线功率顶到上限仍不够,只能切负荷。

同样一个系统,用两阶段鲁棒求解,预测场景下的运行成本大约是1340元,比确定性模型高了约12%。但把它放到同样的最坏场景下,功率平衡、储能SOC、联络线功率全部保持不越限。这12%就是鲁棒解的保护成本,也就是你为了避免小概率极端场景而花掉的保险费。这个成本高低是否可接受,取决于项目本身的可靠性要求。

4.3 关键场景辨别算法到底快了多少

这里给一组我实测的数据:在我的6节点算例里,如果用枚举顶点方式,每个时段风电、光伏、负荷各取上下界,单轮子问题要遍历8个顶点组合,C&CG整体迭代5轮左右才能收敛,总耗时大约4分钟。换成关键场景辨别算法后,每轮子问题只需要解1次内层LP,再用排序函数拼出关键场景,C&CG迭代3到4轮就收敛,总耗时压到了15秒以内。

提速的核心不在于C&CG迭代轮数减少,而在于每轮子问题的计算强度大幅下降。枚举顶点需要对每个顶点都解一次完整的max-min问题,等于把子问题重复了几十遍;辨别算法只解一次LP,再用一个O(n log n)的排序就拿到了最恶劣场景。随着不确定变量数量增加,这个差距是指数级的,算例规模一大,暴力枚举根本跑不动。

4.4 预算参数Γ怎么调

Γ本质上是“你愿意用多少预测精度换取鲁棒性”的旋钮。工程上建议先跑Γ=0的确定性模型,看风电光伏渗透率和联络线容量,然后逐步加大Γ,观察成本上升曲线。如果Γ从0加到6,成本只涨了5%,那说明系统有足够的调节冗余,你可以放心用较大的Γ;如果成本涨了20%以上,就说明系统灵活性不足,要么增加储能容量,要么适当放宽保护程度。

还有一种更实用的做法:用历史数据统计一下,过去一年里实际出力偏离预测超过三个时段的概率有多大。如果这种事件一年也就发生几次,Γ取3到4就够了,没必要让系统天天按最极端的情况来备着。

5. 常见问题与排查技巧实录

5.1 子问题一直不收敛或UB不下降

如果你发现UB一直不下降,或C&CG迭代了十几轮还在震荡,多半是主问题割平面加错了。常见的错误是:把关键场景对应的第二阶段约束加进主问题时,新增变量没有放到独立的元胞数组里,而是和已有变量共用了名字。YALMIP对重复变量很敏感,轻则约束错乱,重则变量维度对不上直接报错。排查方法很简单:在循环里打印每轮割平面的变量名和约束数量,一眼就能看出来是不是在覆盖。

还有一种情况是子问题求出来的关键场景和上一轮相同,导致割平面不产生新约束。这通常是关键场景辨别函数里排序取极值时,没有考虑多个变量系数相等的情况,只取了第一个。建议在排序时加上随机打破平局的处理,或者同时把排名前几位的场景都加入割平面,防止迭代卡死。

5.2 求解器报Infeasible或Unbounded怎么办

子问题对偶变换后如果报无界,通常是强对偶条件不满足。检查第二阶段模型里有没有整数变量,比如储能充放电状态的0-1变量。如果有,就不能直接对偶,需要把整数变量留在外层,或者用KKT条件加互补松弛来处理。很多初学两阶段鲁棒的同学在最开始都会踩这个坑,因为微网里储能SOC和充放电状态确实连在一起。

另一个常见原因是第二阶段问题本身不可行,即给定第一阶段决策后,在某些场景下没有可行解。这可能是联络线功率上限、爬坡约束、SOC边界设置不合理。排查思路:先固定第一阶段解,把最坏场景塞进去,单独解一个确定性LP,看看哪个约束被违反,再用约束的松弛变量定位问题。

5.3 YALMIP和Gurobi的版本兼容问题

新的MATLAB版本有些人装了YALMIP一直报错,90%是路径没配好,或者YALMIP版本对Gurobi接口版本不匹配。先运行yalmip('clear'),再用yalmiptest确认求解器能被识别。Gurobi的新版本对老版YALMIP不兼容时,需要升级YALMIP或者固定Gurobi版本。我用的组合是MATLAB R2023b + YALMIP R20230616 + Gurobi 10.0.3,实测下来很稳。如果你用的MATLAB版本很新,安装路径里带了空格,记得把路径加到MATLAB的set path里,不要直接用cd切过去。

5.4 Big-M的取值技巧

处理对偶问题的边界条件时,Big-M取太小会限制解空间导致次优,取太大会让数值病态。很多代码模板里直接写M=1e6,这在高维问题上容易让Gurobi的数值稳定性变差。建议先跑一次确定性模型,看看每条约束的对偶乘子量级,再按一个数量级放大去取M。比如对偶乘子最大是0.5,M取10左右就够,不要随手写1e6。

5.5 常见问题速查表

现象可能原因排查与解决
C&CG迭代不收敛割平面没有新增有效约束检查关键场景是否重复;加入多个候选场景
子问题报Unbounded对偶条件不满足或M过小检查第二阶段是否有整数变量;调整Big-M
主问题报Infeasible初始场景选择不合理先用预测场景初始化,再进入C&CG循环
UB不下降第二阶段约束写错单独固定第一阶段解做可行性校验
结果过于保守Γ设置过大根据历史数据统计偏离频率,取适中Γ
YALMIP识别不到Gurobi路径或版本不匹配重装匹配版本,运行yalmiptest确认

一个人把两阶段鲁棒微网调度从建模到代码完整跑下来,最大的感受是:关键场景辨别算法并不神秘,它只是把“枚举所有顶点”这种暴力思路换成了“基于对偶乘子排序定位关键顶点”的聪明思路。难的不是算法本身,而是它和C&CG框架的咬合。变量数组怎么组织、割平面怎么加、上下界怎么更新,任何一环出问题都会导致结果乱七八糟。如果你们运行代码时发现迭代次数异常,不要急着怀疑求解器,先从Γ和关键场景辨别函数查起,大概率是场景排序时符号处理出了问题。这套框架后续还可以扩展的方向不少,比如把静态两步推广到多阶段滚动优化,把关键场景辨别改成在线学习更新,或者和深度强化学习结合做分布鲁棒调度,都是可以继续打磨的路子。

需要专业的网站建设服务?

联系我们获取免费的网站建设咨询和方案报价,让我们帮助您实现业务目标

立即咨询