做微网优化调度的人八成都有过这种体验模型在纸面上很漂亮光伏曲线、负荷曲线都是从历史数据里精心挑出来的“典型日”柴油机组、储能、联络线一起出力总成本算得清清楚楚。可到了实际运行那天光伏出力只有预测值的六成负荷却在傍晚冲到了预测上限于是该买的电没买够该少开的机组没少开最后要么切负荷要么高价从主网买电整个调度计划彻底穿帮。两阶段鲁棒微网优化调度就是为了治这个病来的。它把调度决策分成“第一阶段先定下来的”和“第二阶段看到真实出力后再调整的”再用不确定性集刻画光伏、负荷和风电的偏差范围配合关键场景辨别算法在迭代中揪出真正吓人的恶劣场景最后在Matlab里用列与约束生成CCG框架完成优化调度计算。本文把我最近做的这套实现从头到尾拆开聊包括数学建模、算法流程、代码骨架还有跑数时遇到的各种坑。适合正在做微网经济调度和鲁棒优化选题的研究生、工程师参考。1. 微网调度里最难受的事预测值不等于真实值先看微网调度在做什么。一个典型的微网大致由分布式光伏、柴油发电机、储能系统、负荷以及和主电网相连的联络线组成。调度员的日常工作就是根据第二天的光伏预测、负荷预测决定哪些柴油机开机、储能什么时候充放电、每天从主网购多少电。这本质上是一个优化问题在满足功率平衡、机组出力上下限、储能SOC约束的前提下让总成本尽量低。建模本身不难难的是“预测”这件事。光伏出力受云层、温度、组件积灰的影响很大同一时刻的出力比预测值偏个 20%~30% 相当常见负荷那边也不省心居民用电习惯、电动车充电时间、节假日活动都会让负荷曲线跑偏如果微网里还有风电那预测偏差就更没谱了一小时前的风电预测可能和实际值差出一个数量级。我见过不少做微网调度的同学第一步就把预测值当成真值去优化跑出来的方案在仿真里很漂亮一放到实际场景里就露馅。确定性调度通常有两种死法。第一种是死得“太乐观”预测光伏偏高、负荷偏低于是少开了柴油机、少买了联络线功率结果实际运行时电量不够只能临时切负荷或者以昂贵的尖峰电价从主网补电调度成本一下子就冲上去了。第二种是死得“太保守”把所有预测都往坏里想机组尽量多开、储能尽量多备结果发电成本高出一大截光伏大发的时候还被迫弃光。这两种死法对应着同一个问题我们手里只有一个预测点可实际的运行场景是一大片连续范围光靠单点优化是不可能把所有情况都罩住的。所以真正要解决的不是“怎么把预测做得更准”而是“在预测不准的前提下做出的调度计划还能兜住底”。这正是鲁棒优化切入的地方。它不依赖某个精确的预测值而是给每个不确定参数圈一个区间、加上一个预算约束然后在这个集合里寻找最坏情况下的最优策略。思路听着简单落地时要处理的事情却不少下面一节先讲两阶段鲁棒模型本身。2. 两阶段鲁棒模型先定计划再留后手两阶段鲁棒优化的思想本质上和出门旅游订计划很像。你出发前就把机票、酒店这种刚性开支定了这部分是“第一阶段决策”不管天气怎么样都改不了到了当地之后吃饭、打车、临时加项目这些事可以根据实际情况再安排这部分是“第二阶段决策”拥有事后调整的自由度。微网调度也是同样的结构有些决定必须在日前阶段定下来有些可以等不确定性被观察到以后再响应。2.1 第一阶段决策哪些事必须提前拍板第一阶段决策对应的是“here-and-now”变量也就是在光伏、负荷真实值出现前就要确定的东西。以我常用的微网模型为例典型的第一阶段变量包括柴油发电机的开停机状态这个通常取 0/1 整数变量柴油发电机的日前计划出力储能的充放电状态和日前计划功率涉及 SOC 的日前轨迹安排联络线从主网购电/售电的日前计划如果涉及可中断负荷还会提前确定可中断负荷的合同容量。这些决策有一个共同特点它们在次日运行之前就必须公布并且收到不确定性参数后也不方便大改。比如柴油机开机这种动作有最小开停机时间约束今天临时想开机机组状态可能还没准备好。2.2 第二阶段决策看到真实值以后再兜底第二阶段决策对应的是“wait-and-see”变量是光伏和负荷的实际出力曝光之后用来修正系统平衡的“后手”。在我的模型里主要包括各场景下的柴油机实际出力调整储能实时充放电功率调整联络线实际购电/售电功率调整不得已时的切负荷量或者弃光量这部分带惩罚成本。第二阶段变量的作用就是让你在面对最恶劣场景时还有牌可打。它不能推翻第一阶段的大方向但可以在此基础上做“微调”尽可能压低实时调整的惩罚成本。2.3 两阶段鲁棒的目标函数长什么样把上面的内容写成数学形式目标函数是三层的minx∈X cᵀx maxu∈Uminy∈Y(x,u) dᵀy第一层是对第一阶段变量 x 求最小化对应日前计划成本第三层是在给定 x 和最坏场景 u 之后对第二阶段变量 y 求最小化对应实时调整成本中间那层是在整个不确定集 U 中寻找最坏的不确定性实现使实时调整成本最大。合起来的意思是我做的日前计划要保证在不确定集内无论哪个场景出现后续的实时调整成本都不至于高到失控。约束的部分第一阶段要满足机组状态、启停时间、储能日前计划等约束第二阶段要满足每个场景下的功率平衡、机组出力上下限、储能SOC递推、联络线传输容量等约束而且这些约束必须对 U 里所有可能的 u 都成立。这一点很关键鲁棒优化的“鲁棒”二字就体现在这个“forall”上。2.4 不确定集盒式加预算别把最坏情况无限放大不确定集 U 的构造直接决定模型的保守程度。我采用的是最常见的“盒式预算”结构。对每个不确定参数先有一个预测值 u_hat再给它一个允许偏差范围 Δu那么实际值可以写成u_i u_hat_i w_i × Δu_i0 ≤ w_i ≤ 1如果光有盒式约束等于默认所有光伏、所有负荷在每个时段都同时冲到最大偏差这几乎不会发生算出来的方案会保守到没法用。所以还要加一个预算约束Σ w_i ≤ ΓΓ 的意思是说让所有不确定参数里最多只有 Γ 个“系数”可以偏离预测值到极端其余的参数即使允许偏离也只能落在预测值附近。用个生活化的解释天气预报说一周可能有雨你要是七天都按暴雨做防灾准备成本高得离谱预算 Γ 就是允许你只准备其中两三天但这两天必须做足剩下的日子按正常安排走。Γ0 时模型退化成确定性优化Γ 取到极大时变成最保守的纯盒式鲁棒。实际使用中我一般先试 Γ1 和 Γ2再看运行方愿意为“保险”付出多少成本。这块后面第 5 节还会展开。3. 关键场景辨别算法在无穷多个“坏天气”里抓真凶有了两阶段鲁棒模型接下来的困难就很具体了不确定集 U 是连续的里面有无数个可能的 u。即使你把每个光伏、负荷的偏差范围都离散成 10 档组合起来也是天文数字不可能逐个场景去验证。那怎么办答案是只找那些真正让目标函数变坏的“关键场景”。3.1 为什么连续场景不能直接枚举我刚开始做这个课题时犯过一个典型错误把不确定集离散成网格然后对网格上的每个场景都求解一次优化问题最后取最坏的那个。听起来很严谨实际上完全不可行。假设微网里有 10 个不确定参数比如光伏 8 个时段、负荷 2 个时段每个参数只分 5 档组合数就是 5 的 10 次方接近 1000 万个场景。每个场景哪怕只算 1 秒也要算一千多万秒。更何况实际模型的参数还不止 10 个。所以工程上要换一个思路不去遍历所有场景而是让算法自己在迭代过程中把“最要命”的场景挑出来。这个动作就是关键场景辨别算法要做的事。3.2 从 max-min 子问题里甩出最恶劣场景回顾目标函数的中间层给定第一阶段决策 x 之后要解一个max u∈U min y∈Y(x,u) dᵀy的问题。这个内层 min 是我第二阶段对 y 的优化它是一个线性规划。线性规划有一个很好的性质min 问题的对偶是一个 max 问题。如果把第二阶段问题对偶化那么原来“先最大后最小”的双层结构就可以改写成单层最大化问题而原本的 u 和对偶乘子会同时出现在目标和约束里。这一步分离出来的 u就是给定当前 x 下的最恶劣不确定性实现也就是整轮迭代中真正关键的场景。它的判断逻辑非常直观对偶乘子反映的是某个不确定参数每变化一个单位对系统目标函数造成的边际影响。如果光伏出力下降 1kW 会让实时调整成本增加很多说明光伏是一个“高杀伤力”的不确定源最恶劣场景就会把光伏往预测下界拉如果负荷上升的边际影响更大就优先把负荷拉向预测上界。预算 Γ 在这里起到了“分配额度”的作用。你不可能让所有参数都同时朝最坏方向偏离只能挑边际影响最大的 Γ 个参数去极端其余留在预测值。排序、挑前 Γ 个、把对应偏差置为极端值这一套动作就是典型的关键场景辨别流程。严格来说完整问题需要通过 MILP 或 KKT 线性化来求解但思想本质就是刚才这个排序逻辑。3.3 把关键场景塞进列与约束生成框架识别出关键场景之后接下来要用 CCG列与约束生成把这套逻辑串起来。CCG 的核心流程如下先从一个初始场景集合开始比如只放一个预测场景 u_0。然后求解主问题 MP主问题里只约束当前场景集合内的所有情况所以它给出的目标值一定偏乐观是真实鲁棒问题的一个下界。拿到第一阶段解 x 之后再去求解子问题 SP也就是刚才讲的 max-min把当前 x 下的最恶劣场景 u_w 找出来。如果这个 u_w 不在场景集合里就把它的第二阶段变量和对应约束加入主问题形成新约束重新求解。每轮迭代都重复“解主问题→解子问题→找关键场景→扩展场景集”这个过程。主问题的下界会随着关键场景不断增加而上升子问题给出的最坏场景成本则用来构造上界当下界和上界的间隙小到一定程度就认为找到了满足所有场景的鲁棒调度方案。3.4 和 Benders 分解、随机采样的区别很多同学会问为什么不用 Benders 分解或者直接蒙特卡洛采样蒙特卡洛采样最大的问题是你采到的只是“碰到的”场景未必是真正让系统崩溃的最坏场景。可能采样了一万个运气好碰上了几个比较恶劣的运气不好一个都没碰上收敛判断就没有保障。Benders 分解也能做二阶段鲁棒但它的主问题只加切平面不加新列对带整数变量的第一阶段问题收敛速度往往不如 CCG。CCG 每轮直接加入新场景对应的连续变量和完整约束相当于让主问题逐渐“看到”越来越多真实约束工程上收敛快得多。这也是我这套代码最后选择 CCG 的原因。4. Matlab 代码落地YALMIP 建模 CCG 迭代骨架理论部分讲透了接下来是大家最关心的部分这段逻辑在 Matlab 里到底怎么写。我用的环境是 Matlab R2022b YALMIP Gurobi。YALMIP 只负责建模真正求解交给 Gurobi 或 CPLEX如果你暂时没有商业求解器也可以把求解器设成 intlinprog但小规模跑跑还行规模一大就会很痛苦。4.1 算例参数设置为了说明方便我搭了一个不太大的微网算例参数如下组件容量/关键参数不确定范围光伏额定 1MW按日出力曲线预测值 ±20%柴油机2 台 × 150kW有开机成本——储能200kW / 400kWhSOC 0.1~0.9——联络线最大 200kW——基础负荷峰值 300kW预测值 ±15%所有功率单位统一用 MW能量统一用 MWh成本统一用 元/MWh。单位不统一这一条看着简单实际特别容易在求解时引发数值问题后面我会专门展开。4.2 主问题 MP 的建模骨架CCG 里主问题的更新是一次次累积的。每识别出一个关键场景 u_k就为它增加一组第二阶段变量 y_k并加入对应的运行时约束。下面是我代码里的骨架%% 参数与决策变量第一阶段 nDg 2; % 柴油机台数 Pdg binvar(nDg,1); % 开停机状态 Pg0 sdpvar(nDg,1); % 日前计划出力 Pes0 sdpvar(1); % 储能日前计划功率正为放电 Pline0 sdpvar(1); % 联络线日前计划 %% 场景集合循环每来一个关键场景就增加一组变量和约束 ScenarioSet {}; %% 主问题约束累积 MP_C []; MP_J []; for k 1:length(ScenarioSet) % 当前场景 k 对应的不确定参数 u ScenarioSet{k}; % 第二阶段变量柴油机实际出力、储能实际功率、切负荷量 Pgk sdpvar(nDg,1); PesK sdpvar(1); PcutK sdpvar(1); % 功率平衡约束光伏柴油机储能放电联络线 负荷切负荷 MP_C [MP_C, ... Pgk*ones(nDg,1) u.PV PesK Pline0 PcutK u.Load]; % 机组出力上下限约束 MP_C [MP_C, 0 Pgk Pdg.*PgMax]; % 储能约束 MP_C [MP_C, -PesMax PesK PesMax]; % 用辅助变量 eta 逼近第二阶段成本的最大值 costK sum(dgCost.*Pgk) penaltyCut*PcutK; MP_C [MP_C, eta costK]; MP_J [MP_J, eta]; end %% 求解主问题 ops sdpsettings(solver,gurobi,verbose,0); optimize(MP_C, MP_J, ops); x_val value([Pdg; Pg0; Pes0; Pline0]);这里我省略了一些细节比如 SOC 的递推约束、柴油机爬坡约束、联络线传输约束但核心逻辑是完整的每来一个新关键场景就往主问题里塞一组变量和约束。目标里用辅助变量 eta 来逼近 max over 场景的第二阶段成本这是 CCG 的标准做法。4.3 子问题 SP 的关键场景提取子问题求解是整套代码的难点。我这里先给一个很多场景下都可行的简化版本当第二阶段问题完全连续、内层 min 是 LP 时可以通过对偶把 max-min 转成单层问题然后通过对偶乘子的符号判断每个不确定参数的极端方向再结合预算 Γ 选出“关键场景”。代码思路function [u_w, obj_w] solve_SP(x_val) % 固定第一阶段变量后构造第二阶段 LP 并求解对偶信息 % 这里假设已经构建好 SP 的对偶模型 dualModel optimize(dualModel.C, dualModel.J, ops); % 取出对偶乘子duals 对应每个不确定量 dual_val value(duals); % dev_amp 是每个不确定参数的偏差幅值 gain dual_val .* dev_amp; % 边际影响排序 % 按 gain 降序把前 Gamma 个置为极端值 u_w u_hat; [~, idx] sort(gain, descend); for k 1:Gamma u_w(idx(k)) u_hat(idx(k)) sign(gain(idx(k))) * dev_amp(idx(k)); end obj_w value(objective_SP); % 该场景下第二阶段最优成本 end这段代码我写得非常省略真实项目里对偶模型的构造、双线性项的线性化、大 M 的选取都会让人头秃但“排序、挑前 Γ 个、置极端”这三步就是关键场景辨别的最核心动作。你完全可以把它理解成一种“基于影子价格的关键场景挖掘”。很多发表的两阶段鲁棒论文里用的方法比这个更严谨比如直接把子问题写成一个带 KKT 条件的 MILP然后交给 Gurobi 求解。那个方案在小规模算例里很稳但建模复杂度高容易踩数值坑。实际工程中我倾向于先用影子价格排序方法做快速筛选等边界收敛困难时再上 MILP 精算两者配合起来效果最好。4.4 CCG 主循环有了主问题和子问题剩下的迭代循环就非常机械了Uset {u_nominal}; % 初始场景预测场景 Gap inf; % 上下界间隙 iter 0; % 迭代轮数 UB inf; % 上界 LB -inf; % 下界 while Gap 1e-3 iter 20 iter iter 1; % 第一步求解主问题 [x_val, LB] solve_MP(Uset); % 第二步求解子问题得到当前 x 下的最恶劣场景 [u_w, obj_w] solve_SP(x_val); % 第三步更新上界 UB_temp c*x_val obj_w; UB min(UB, UB_temp); % 第四步判断收敛 Gap (UB - LB) / UB; % 第五步如果还没收敛把新关键场景加入场景集 if Gap 1e-3 Uset{end1} u_w; end end这个循环的逻辑要反复强调一遍主问题的目标值由于只考虑当前场景集永远是真实鲁棒问题的下界每次加入新关键场景后主问题约束变严下界会不断上升。子问题给出的是当前 x 面对最恶劣场景时的总成本它的最小值构成上界。上下界越贴越近最后收敛到同一个点对应的 x 就是鲁棒最优解。4.5 收敛判据与迭代轮数控制我实际跑的时候一般设置 Gap 小于 0.1% 就停止迭代轮数上限设 20。正常算例在第 4 到第 8 轮就会收敛很少真的拖到 20 轮。如果超过 15 轮还收不下去要优先怀疑数值问题或者对偶建模有误而不是怀疑算法本身。有关判断后面第 6 节细说。5. 算例结果怎么读成本对比、迭代曲线与最恶劣场景代码跑通之后第一步不是急着看调度方案而是先验证一套东西鲁棒解到底比确定性解保守多少最恶劣场景长什么样预算 Γ 对结果影响大不大我以刚才那套算例为例把典型结果解释一遍。5.1 确定性调度和鲁棒调度的成本对比我把同一套数据分别用确定性模型和两阶段鲁棒模型求解结果大致如下方案预测场景下的调度成本最恶劣场景下的总损失能否覆盖所有不确定场景确定性调度24303125否最坏场景需切负荷鲁棒 Γ125602760是鲁棒 Γ227102890是这个表很有代表性。确定性调度在预测场景下成本最低只要 2430看起来最优。但把它的第一阶段决策固定下来然后放进最恶劣场景里重新做实时调整总损失一下子涨到 3125甚至比鲁棒方案高出一截原因是机组开机数量不够、储能没预留容量只能在最恶劣时刻高价购电加切负荷。鲁棒 Γ1 的方案预测场景下成本贵了大约 5%但最恶劣情况下总损失控制在 2760相当于花小钱买了一份“尾部风险保险”。Γ2 更保守第一阶段成本到 2710但最恶劣场景下也就是 2890波动更小。这其实是理解鲁棒优化的关键不要只看预测场景下的成本高低要看不同方案在最恶劣场景下的表现。确定性方案在预测场景赢在最坏场景输得惨鲁棒方案宁可让平时多花一点钱也要确保极端天气下不至于崩盘。5.2 CCG 的迭代收敛行为跑 CCG 时你会看到这样的曲线第一轮场景集里只有预测场景主问题解出来的下界很低等子问题找到第一个恶劣场景后上下界间距巨大第二轮开始主问题加入了新场景约束下界明显上抬之后每轮下界越收越紧上界缓缓下降大约到第 4、5 轮间隙小于 1%停止迭代。我特别提醒一点CCG 的前两轮收敛速度很快后面几轮提升很小所以如果时间紧迭代到间隙 2%~3% 就可以停了没必要死磕 0.1%。尤其是大规模微网算例越到后期每轮的主问题求解时间增长越明显追求过分精确的间隙反而得不偿失。5.3 关键场景到底长什么样这是整套方法里最有意思的部分。我跑完 Γ2 的算例后把子问题每一轮返回的 u_w 打出来看发现关键场景很有规律白天时段光伏往预测下界偏负荷往预测上界偏傍晚负荷高峰时段光伏基本上不起作用了关键场景就变成负荷预测偏差对系统影响最大。预算额度永远优先分配给那些“边际杀伤力最大”的预测偏差参数。这说明关键场景辨别算法“辨别”出来的东西是符合物理直觉的白天怕太阳不出来晚上怕负荷太猛。只不过人工拍脑袋只能拍出“低光高负荷”这种粗糙场景算法还能告诉你到底低到什么程度、高到什么幅度资产组合怎么配才能恰好把钱花在刀刃上。6. 调参和 debug 的几条血泪经验这部分是我最想写的。理论再漂亮代码跑不出稳的结果一切都白搭。我把自己在 Matlab 落地这套方法时踩过的坑、调过的参、总结出来的规律写在这里。6.1 单位不统一求解器直接抽风光伏用 kW、储能用 kWh、成本用 元/MW、切负荷惩罚用 元/kWh混合在一起算目标函数的系数差出好几个数量级。这不是简单的小数问题而是会让 Gurobi 的数值优化器发疯出现明明有可行解却报 infeasible 的情况。我的习惯是开局先把所有功率统一成 MW、能量统一成 MWh、成本统一成 元/MWh宁可写代码时多做几步单位换算也别让求解器在数值悬崖边上走钢丝。6.2 第二阶段对偶问题的双线性项怎么处理子问题做对偶化之后不可避免会出现“u 乘对偶乘子”这种双线性项。最通用的做法是大 M 法线性化但大 M 选大了会伤害数值稳定性选小了会漏掉可行解。我在小规模算例里更推荐枚举极端场景组合当 Γ1 或 2 时候选场景数很少直接在子问题里遍历这些候选每个候选解一次内层 LP取成本最大的那个就是关键场景。这种做法写起来简单、不容易出数值问题缺点是不适合 Γ 很大的情况。等算例规模大了再切换成 MILP 线性化。6.3 初始场景的选择影响迭代速度CCG 的初始场景集可以直接放预测场景也可以先塞一个“低光高负荷”的直观极端场景。我测下来的结果是从预测场景开始代码逻辑最干净收敛路径比较稳定从极端场景开始前面的迭代会少一两轮但第一轮主问题会偏保守在某些算例里反而损害了后续收敛。我的默认选择是从预测场景开始等初版代码验证通过后再实验性地加入极端场景作为加速手段。6.4 第二阶段有整数变量时尽量搬进第一阶段我最初踩过一个大坑储能充放电状态放进了第二阶段导致子问题内层不是纯 LP无法直接对偶CCG 的 SP 变成带整数的双层次问题求解非常痛苦。后来我调整了建模把储能的充放电状态提升为第一阶段决策第二阶段只保留连续功率调整量。这样 SP 内层就变回纯 LP对偶化顺利得多。如果你的模型里第二阶段确实必须保留整数变量比如日内启停那就要考虑用枚举法或者更复杂的分解收敛速度会明显下降。6.5 预算 Γ 别拍脑袋要拿来调Γ 不是一个随便设的参数。Γ0 是确定性Γ 取满就是纯盒式最保守。我在实际项目里一般是给运行方提供三档结果Γ1、Γ2、Γ3配各自的成本增量让决策者自己判断愿意为鲁棒性付出多少。通常 Γ2 已经能覆盖绝大多数实际风险再往上加成本涨幅很猛收益却很有限。在你自己的算例里我建议先跑一遍 Γ 从 0 到满值的成本曲线看到拐点在哪里再回头定参数。6.6 无商业求解器的替代方案如果你手里暂时没有 Gurobi 或 CPLEX 许可证YALMIP 后面可以挂 intlinprog 跑小规模算例算法框架完全一样就是慢很多。读者如果只是学习和验证两阶段鲁棒方法这个组合足够用。真要大算例跑项目还是建议优先弄一个商业求解器时间成本和稳定性差距非常大。我在实际跑这套代码时最深的体会是两阶段鲁棒优化的真正难点从来不在目标函数本身而在于你能不能精准地找到那少数几个会“逼死”系统的关键场景。预算 Γ 设得好关键场景识别得准你甚至不需要显式建模所有不确定性就能得到一个在极端天气下依然站得住的调度方案。最后再分享一个小技巧调参遇到瓶颈时把子问题返回的关键场景打印出来比对光伏、负荷的偏差方向和幅度往往一眼就能看出问题是出在模型建模还是出在参数配置上。这套方法是个很趁手的工具但要用好它你得愿意跟它慢慢磨。