做电力系统调度优化的朋友应该都体会过那种“模型算出来挺漂亮一到实际就跑偏”的无力感。风电预测给的是100MW实际出力可能只有60MW光伏中午预报大晴天下午一片云飘过来出力直接腰斩负荷就更不用说了居民区和工业区的用电习惯根本不是一条平滑曲线。传统的确定性调度模型把这些预测值当成板上钉钉的常数算出来的机组组合、备用容量、潮流分布一旦遇到风光出力波动轻则经济性大幅恶化重则出现切负荷甚至潮流越限。这也是为什么近些年鲁棒优化在电力系统调度里越来越受重视——它不赌预测值的准确性而是把不确定性“圈”在一个集合里保证集合内所有可能场景下方案都可行。这篇文章想分享的就是一套我实际跑通的计及风、光、负荷不确定性的两阶段鲁棒优化完整实现建模用YALMIP求解器用Gurobi核心算法是CCG列与约束生成中间的双线性项用大M法做线性化。内容会从问题建模一路讲到Matlab代码细节和踩坑经验适合正在做微电网调度、园区综合能源优化、电力系统日前计划这类课题的研究生和工程师参考。1. 为什么确定性调度到了含新能源的系统里就不够用了1.1 确定性模型的“赌徒心态”传统经济调度模型里风电、光伏、负荷都被写成了确定数值。模型长这样目标是最小化总运行成本约束包括功率平衡、机组出力上下限、爬坡速率、网络安全等。求解器给出的结果是一个“在预测值完全准确时才最优”的方案。问题在于新能源的预测误差不是小概率事件而是常态。风电短期预测误差的均方根误差普遍在10%~20%之间极端天气下甚至能到30%以上。如果日前计划按预测值100MW安排火电出力实际风电场发到60MW那40MW的功率缺口就得靠实时市场的备用、储能或者紧急爬坡来补。要是备用容量不够就得限负荷。这种“赌预测准”的思路在传统火电为主的系统里问题不大因为负荷预测相对精准机组调节能力也强。但风光占比提高之后系统的不确定性来源从一个变成了三个预测误差的叠加效应让确定性方案的实际运行成本远高于模型计算值。1.2 鲁棒优化的核心思想不赌但做好准备鲁棒优化的思路完全不同。它不依赖预测值的概率分布而是给每个不确定参数划定一个变化范围构造一个不确定集合Uncertainty Set然后寻找一个对所有集合内场景都可行的方案。说句大白话确定性优化是“我赌明天不下雨所以只带一把伞”鲁棒优化是“我不知道明天下不下雨所以我包里常备雨衣真下大了也不慌”。代价是雨衣占地方、有点重——对应到调度里就是鲁棒方案的经济性会比确定性方案差一些多出来的成本就是为了换取“不管风光怎么波动系统都安全运行”的保障。这个保守度是可以调节的。通过预算约束Budget Constraint我们可以限制“最坏情况下同时偏移到极端的不确定源数量”。比如三个不确定源预算Γ取1意味着优化只保证最多一个源同时取到极端值Γ取3则要求所有源同时取极端值时方案依然可行。Γ越大方案越保守成本越高。1.3 两阶段结构先定下来再调整两阶段鲁棒优化对应的是电力系统调度的实际时序逻辑第一阶段日前决策机组启停、备用容量预留、储能充放电计划等必须提前确定下来。这些决策变量在不确定参数实现之前就要敲定特点是“定了就改不了”。第二阶段实时调整等风光出力、负荷的实际值揭晓后系统在给定第一阶段决策的基础上进行再调度通过调整火电出力、储能充放电、切负荷等手段以最小代价保证功率平衡和网络安全。所以模型的目标函数由两部分构成第一阶段的启动成本和基础运行成本加上第二阶段在最坏场景下的调整成本。外层是第一阶段决策的min问题内层是一个max-min嵌套——先枚举所有可能的不确定场景max再求每个场景下的最小再调度成本min。2. 两阶段鲁棒优化的数学模型拆解2.1 不确定集合的构造与离散化风、光、负荷的不确定集合最常用的是盒式集合加预算约束。先把每个不确定源写成名义值加偏差的形式P_w ∈ [P_w^0 − ΔP_w, P_w^0 ΔP_w] P_v ∈ [P_v^0 − ΔP_v, P_v^0 ΔP_v] P_l ∈ [P_l^0 − ΔP_l, P_l^0 ΔP_l]其中P_w^0是风电预测值ΔP_w是允许的最大偏差三个不确定源的范围共同构成一个超立方体。实际建模时我不建议直接在这个连续区间上处理。原因有两个一是连续区间与第二阶段对偶变量的乘积会构成复杂的双线性项线性化非常麻烦二是从工程角度看真正让调度员头疼的是“某个源出力取到上界还是下界”中间值的影响往往是单调的。所以实操中更常见的做法是引入0-1变量把每个不确定源离散成两个极端状态u_i u_i^0 Δu_i · z_i, z_i ∈ {0, 1}z_i 0表示该不确定源取下界z_i 1表示取上界。再叠加预算约束Σ z_i ≤ ΓΓ是预算参数控制同时取极端的不确定源数量上限。离散化之后不确定集合变成了一个包含有限个极端场景的组合集合这正好为后续CCG算法中“把最坏场景带回主问题”提供了天然的离散场景支持。三个不确定源的具体参数化在算例里常这样设置以下数值仅为示范实际应基于预测数据统计得到不确定源名义值偏差范围离散状态风电出力80 MW±16 MW (20%)64 MW / 96 MW光伏出力50 MW±10 MW (20%)40 MW / 60 MW负荷需求200 MW±20 MW (10%)180 MW / 220 MW2.2 第一阶段决策模型第一阶段决策变量主要包括机组启停状态u_g0-1变量机组启动动作变量用于计算启动成本预留的旋转备用容量储能的充放电状态与SOC初值如果有储能第一阶段模型可以写成min Σ (启动成本 空载成本) ηs.t. 启停逻辑约束、最小启停时间约束、机组出力上下限约束、备用容量约束这里的η是一个辅助变量在CCG迭代中用来逼近第二阶段的期望最坏场景成本这也是两阶段鲁棒优化中“第二阶段成本被第一阶段决策耦合”的标准处理手法。2.3 第二阶段再调度模型给定第一阶段决策x机组启停状态、备用等之后第二阶段模型针对某一具体的不确定场景u求解最小再调度成本min Σ (机组出力调整成本 切负荷惩罚 弃风弃光惩罚)s.t. 功率平衡约束 线路潮流约束直流潮流或交流潮流线性化 机组爬坡速率约束 切负荷量不超过该节点负荷需求 弃风弃光量不超过该时刻风光出力写成紧凑的矩阵形式会更方便后续的对偶推导Q(x, u) min_y b^T y s.t. D y ≥ h − A x − E u y ≥ 0其中y代表第二阶段的连续决策变量机组出力调整量、切负荷量、弃风弃光量等h是常数向量A、D、E是系数矩阵u是不确定向量。2.4 主问题与子问题的耦合结构CCG算法把完整的两阶段鲁棒问题拆成主问题MP和子问题SP交替求解主问题在第k次迭代时min c^T x η s.t. 第一阶段约束 η ≥ b^T y_j, j 1, 2, ..., k D y_j ≥ h − A x − E u_j^*, j 1, 2, ..., k y_j ≥ 0这里u_j^*是前k次迭代中子问题找到的“最坏场景”。每迭代一次主问题就会多一组变量y_j和对应的约束。这个机制通俗点理解就是主问题一开始“乐观”地只考虑已知的场景子问题每次迭代都会找出一个让主问题方案吃瘪的坏场景把它“塞”回主问题逼着主问题改进决策。子问题给定第一阶段决策x^*后Q(x^) max_{u ∈ U} min_y { b^T y | D y ≥ h − A x^− E u, y ≥ 0 }这就是经典的max-min双层结构。外层的max寻找最坏场景内层的min求解该场景下的最小再调度成本。这个双层结构不能直接求解需要借助强对偶定理把它压平成单层优化问题。3. 大M法把双层嵌套压平的关键线性化3.1 强对偶转换子问题中内层的min是一个线性规划满足强对偶条件原问题可行且有界时对偶问题最优值相等。我们对内层min取对偶得到Q(x^) max_{u, λ} λ^T (h − A x^− E u) s.t. D^T λ ≤ b, λ ≥ 0 u ∈ U这是一个纯max问题变量是对偶变量λ和不确定变量u。到这里max-min的双层嵌套变成了单层max可以直接扔给求解器。3.2 双线性项的出现与大M线性化但注意目标函数λ^T (h − A x^* − E u) 里藏着麻烦展开后会有λ^T E u这样的项即对偶变量与不确定变量的乘积。λ是连续变量u经过离散化后是0-1变量λ_i × z_j这类项是双线性项非线性非凸Gurobi和CPLEX都直接处理不了。处理思路是利用大M法把“连续变量 × 0-1变量”的乘积线性化。引入辅助变量v_ij λ_i · z_j用以下一组线性不等式刻画v_ij的行为0 ≤ v_ij ≤ M · z_j λ_i − M(1 − z_j) ≤ v_ij ≤ λ_i M(1 − z_j)这组约束的含义是当z_j 0时第一个约束把v_ij强制压成0当z_j 1时第二个约束把v_ij限制在λ_i的邻域内当M足够大时v_ij就等于λ_i。可以把这个机制理解成“开关”——z_j是开关打开时λ_i通过关闭时乘积归零。项目实践中我曾经绕过离散化、直接用连续u变量做线性化需要引入更复杂的多面体近似实现难度高而且不稳定。相比之下0-1离散化大M法是目前最成熟、最稳的方案特别是与预算约束天然搭配所以我一直推荐这条路。3.3 大M取值的实操经验大M法看起来简单真正翻车最多的地方就是M的取值。M取得太小会错误地切掉可行解导致子问题求出来的目标值偏小CCG迭代出来的结果不对M取得太大会造成数值病态Gurobi求解时会出现“Model may be unreliable”的警告矩阵条件数恶化求解速度大幅下降甚至出现错误的整数解。根据我的实操经验M的取值可以参考以下逻辑先不加入大M线性化约束求解一个松弛的子问题观察对偶变量λ的范围M取λ_i可能最大值的2到5倍一般情况下电力系统调度模型里M取1e3到1e4是比较安全的区间尽量不要超过1e6超过1e6之后Gurobi的数值稳定性会明显变差如果模型做完标幺化处理所有功率除以基准值λ的量级通常在0到几十之间M取500到1000就够用了。这个细节直接决定求解稳定性值得多花几分钟做一次敏感性测试。4. CCG算法主问题和子问题怎么互相“喂料”4.1 为什么选CCG而不是Benders分解两阶段鲁棒优化早期常用Benders分解每次迭代只给主问题加一条割平面约束。CCG的区别在于它不仅加约束还向主问题中直接添加一组完整的第二阶段决策变量和对应的场景约束。每轮迭代主问题的规模会增加但换来的是割平面质量显著提升迭代次数大幅减少。打个比方Benders分解像是一个学生每次考试只被告诉“总成绩不及格”需要自己猜哪道题错了CCG则直接告诉他“这次考试哪道题错了、错在哪、正确答案应该是什么”。虽然每次批改的工作量更大但学生进步快得多。对于两阶段鲁棒优化这种场景CCG已经是事实上的标准算法收敛速度比Benders快一个数量级。4.2 算法伪代码完整流程如下初始化UB ∞LB −∞k 1选取初始场景u_1^*通常取名义值即所有不确定源取预测值求解主问题得到最优解(x_k^, η_k^)更新下界LB max(LB, c^T x_k^* η_k^*)固定x_k^求解子问题得到最坏场景u_{k1}^和最优值Q(x_k^)更新上界UB min(UB, c^T x_k^ Q(x_k^*))如果(UB − LB) / UB ≤ εε通常取0.01或更小则停止迭代输出当前解否则将新场景u_{k1}^*对应的变量y_{k1}和约束η ≥ b^T y_{k1}、D y_{k1} ≥ h − A x − E u_{k1}^*加入主问题k k 1回到第3步这个流程里有个微小但很重要的细节更新UB时每一轮都要取历史所有子问题值中的最小值而不能只取当前值。原因是子问题在给定x_k^*下求出的最坏场景成本本来就是一个可行方案的真实成本上界历史上任何一轮的x_k^*都是一个可行的第一阶段决策对应一个上界。取最小值能得到更紧的上界有助于提前收敛。4.3 收敛判据的工程处理理论上UB和LB会随着迭代逐渐逼近实际代码里我还会加两个辅助收敛条件如果连续三轮迭代中新场景u_k^*与已有场景的差异小于某个阈值比如所有场景变量的汉明距离为0说明子问题已经找不到新信息了提前终止避免无谓的迭代开销设置最大迭代次数比如30次防止数值问题导致的死循环子问题求解时有一个常见的坑如果对偶问题不可行或有界性异常多半是原问题第二阶段约束A、D矩阵写错了或者是M取值过小导致线性化约束过度收紧。排查时建议先把预算Γ设为0退化为确定性场景如果此时子问题正常再逐步增加Γ缩小问题定位范围。5. Matlab YALMIP代码实现与踩坑记录5.1 环境与求解器配置我的环境是Matlab R2022b YALMIP Gurobi 10.0。Gurobi求解MILP的能力业界领先CPLEX也可以但Gurobi在鲁棒优化这类含大量0-1变量和辅助变量的问题上表现更稳。YALMIP的安装不用多说重点说Gurobi参数配置。在CCG迭代中主问题和子问题都是MILP有些参数值得重点关注ops sdpsettings(solver, gurobi, ... verbose, 0, ... gurobi.mipgap, 1e-4, ... gurobi.numericsfocus, 1, ... gurobi.timelimit, 300);numericsfocus设为1或2可以在出现数值警告时显著提升求解稳定性代价是求解时间略微增加在鲁棒优化这种数值敏感性较高的问题里非常值得。timelimit一定要设否则遇到病态模型Gurobi可能跑几小时不出结果。5.2 主问题代码骨架主问题用YALMIP建模的核心逻辑如下% 第一阶段变量 u_start binvar(n_g, 1); % 机组启停 start_action binvar(n_g, 1); % 启动动作 eta sdpvar(1, 1); % 第二阶段成本逼近变量 % 第一阶段约束启停逻辑、最小启停时间、备用容量等 Constraints []; % 历史割集合每轮迭代追加 % u_history 是 n_unc × k 的矩阵存储历史最坏场景 for j 1:k y_j sdpvar(n_y, 1); % 第j轮添加的第二阶段变量 Constraints [Constraints, ... eta b * y_j, ... D * y_j h - A * u_start - E * u_history(:, j), ... y_j 0]; end objective c * u_start eta; optimize(Constraints, objective, ops);这段代码的关键在于u_history(:, j)是已知常数是子问题在历史迭代中求出来的最坏场景。主问题的决策变量是u_start、eta和所有历史添加的y_j。每轮迭代y_j会多一组主问题规模线性增长。5.3 子问题代码骨架子问题经过强对偶大M线性化之后YALMIP代码结构如下% 输入当前第一阶段解 u_start_star lambda sdpvar(n_con, 1); % 对偶变量 z binvar(n_unc, 1); % 不确定变量的0-1表示 % 对偶可行域 Constraints [D * lambda b, lambda 0]; % 不确定集合与预算约束 Constraints [Constraints, sum(z) Gamma]; % 大M线性化v lambda * z v sdpvar(n_con, n_unc); M 1000; % 根据λ的量级调整 for i 1:n_con for j 1:n_unc Constraints [Constraints, ... 0 v(i, j) M * z(j), ... lambda(i) - M * (1 - z(j)) v(i, j) lambda(i) M * (1 - z(j))]; end end % 目标函数λ^T(h - A*x*) - sum(sum(E .* v)) % YALMIP默认求解min因此max问题取负号 Objective lambda * (h - A * u_start_star) - sum(sum(E .* v)); optimize(Constraints, -Objective, ops); % 提取最坏场景 u_star u_nominal delta_u .* value(z); Q_star value(Objective);这代码里有几个关键点值得多说一嘴。大M法的核心是v(i, j)这个辅助变量n_con代表对偶约束数量n_unc代表不确定源数量通常等于3乘以时段数如果按每个时段的风、光、负荷都独立不确定来建模的话。注意这里我没有用双循环以外的向量化写法来简化展示实际项目里建议用repmat或kron做向量化当n_con和n_unc变大时双循环会非常慢。另一个容易踩的坑是目标函数里E .* v的符号。E矩阵是第二阶段约束中u的系数矩阵它可能是负的也可能有正有负取决于原问题的结构。目标函数整理成λ^T(h − A x^* − E u)之后展开项−λ^T E u对应到辅助变量就是−sum(sum(E .* v))。符号反了子问题求出来的“最坏场景”就变成了“最好场景”CCG会完全失效。我第一次实现时就在这里卡了一整天最后是通过对比确定性场景Γ0下的子问题和直接优化结果才发现的。5.4 迭代主循环主循环的写法相对直接Gamma 2; % 预算参数 max_iter 30; epsilon 1e-3; UB inf; LB -inf; k 0; % 初始场景名义值 u_his u_nominal; while k max_iter k k 1; % 求解主问题 optimize(mp_constraints, mp_objective, ops); x_k value(u_start); eta_k value(eta); LB max(LB, value(mp_objective)); % 求解子问题 [u_new, Q_k] solve_subproblem(x_k, Gamma, ops); UB min(UB, c * x_k Q_k); % 收敛判断 if (UB - LB) / abs(UB) epsilon break; end % 向主问题添加新场景 u_his [u_his, u_new]; add_new_cut_to_mp(u_new); % 追加 y_k 变量和相关约束 end注意一个容易忽略的细节LB初始化用−inf但如果主问题在第一轮就不可行说明第一阶段约束写得有问题需要回头检查是不是某个约束把所有机组都禁用了或者备用容量要求过严。实际项目中我遇到过主问题第一轮就不可行的情况排查下来是备用容量约束里留的余量太大超出了所有机组的调节能力总和。5.5 高频问题排查表每次在社区里看到有人问“为什么我的CCG不收敛”或者“为什么结果明显不对”基本都是下面这几个原因之一症状可能原因排查方法子问题目标出现负值强对偶转换后符号弄反退化为确定性场景对比直接优化的结果主问题不可行第一阶段约束过严或冲突去掉备用约束测试逐一放回排查UB和LB长期不收敛M值过小或过大输出λ的范围按3.3节调MGurobi警告数值问题量纲不一致或M过大做标幺化减小M设置numericsfocus子问题耗时爆炸双循环大M约束效率低向量化v约束或减少不确定源个数结果与确定性模型相同预算Γ写成了0或场景初始化有问题检查Gamma赋值和u_his初始化5.6 初始化技巧一个容易被人忽略但非常有效的技巧是用确定性场景的结果初始化CCG。具体做法是先令所有不确定源取名义值求解一个确定性两阶段问题把得到的机组启停方案作为CCG主问题第一轮迭代的初始可行解。这样第一轮子问题就能给出一个质量相当不错的上界能明显减少迭代次数。如果模型中包含储能系统建议把储能的初始SOC也纳入第一阶段决策变量处理。因为SOC的初值直接影响后续所有时段的充放电可行域如果放到第二阶段再决定会破坏“日前决策实时调整”的时序逻辑。6. 从能跑到能用的进阶建议6.1 预算Γ的调参逻辑预算参数Γ是鲁棒优化模型提供的一个“保守度旋钮”。Γ0时模型退化为确定性优化成本最低但对不确定性没有任何抵御能力Γ等于不确定源总数时模型最保守要求所有源同时取极端值时方案依然可行成本最高但安全性最强。实际工程中我会这样给决策者分析跑一组Γ从0到最大值的结果画出成本-保守度曲线。曲线上每个点的斜率可以理解为“再买一单位安全性的边际成本”。大多数情况下Γ从小往大增加时成本会先陡增后趋缓。拐点位置的Γ值往往就是推荐的折中选择——既能覆盖绝大部分风险又不会让成本高到离谱。6.2 与大M法相关的几个数值稳定性经验大M法是最容易引起数值隐患的环节。除了前面说的M取值还有几个经验值得分享所有变量尽量统一量纲功率用MW成本用元不要混着用kW和万元如果模型中同时存在数量级相差1e6以上的系数比如M取1e5而其他约束系数在个位数Gurobi的内置预求解器有时会自行简化模型反而导致线性化约束被错误消去。这种情况可以在Gurobi参数里关掉Presolve的部分功能或者把M缩小大M法的M值即使是同一个模型不同约束上也可以不一样。比如某个对偶变量λ_i理论上限是100那么对应的M取300就足够另一个λ_j上限是1000M可以取3000。分别给值比统一一个大M要稳得多6.3 算例验证的推荐做法验证代码时我建议按以下顺序推进先在修改的6节点系统上测试三个不确定源分别对应风机节点、光伏节点和负荷节点。预测值随便设一组合理的数偏差范围按10%~20%取。先跑Γ0的确定性情形确认主问题和子问题在名义场景下求解正常再逐步增加Γ。然后对比确定性方案和鲁棒方案在极端场景下的表现把确定性方案代入子问题得到它在最坏场景下需要付出的高额再调度成本可能伴随切负荷再把鲁棒方案代入成本会低得多。这个对比是论文里最常用的结果展示也是验证鲁棒优化价值最直观的方式。等6节点算例跑通后再移植到IEEE 14节点或118节点系统。移植时唯一要注意的是网络约束矩阵的维度和稀疏性YALMIP里用稀疏矩阵构造约束比用全矩阵快很多内存占用也小很多。套用这套框架实际做过的项目里6节点系统一般5轮以内收敛118节点系统大概需要10到20轮单轮求解时间视Gurobi性能从几秒到几分钟不等。如果迭代超过30轮还没收敛几乎可以肯定是代码里有bug别怀疑是算法本身的问题。6.4 扩展方向这套两阶段鲁棒优化的框架本身是非常通用的。换一个不确定集合的定义方式比如用数据驱动的多面体集合代替盒式集合换一个目标函数比如加入碳排放约束或者把第二阶段的线性规划换成混合整数规划比如考虑储能充放电状态的实时调整大框架都不需要动只需要改对应的约束和变量定义。更进一步如果想考虑不确定性之间的时空相关性可以考虑用椭球不确定集合或者数据驱动的不确定集合这时大M法的角色会弱化但CCG主-子问题迭代的框架依然成立。最后说一点个人体会。这套方法表面上看是数学建模和代码实现的功夫真正吃透之后最有价值的其实是它对“不确定性成本”的量化能力。当决策者问“鲁棒方案为什么比确定性方案贵这么多”用一条Γ-成本曲线就能讲清楚你愿意花多少钱换来对多少不确定性的保障。这种把工程决策翻译成优化语言的能力是做电力系统优化工作的人最值钱的本事。 两阶段鲁棒优化在电力系统调度中的实战：计及风光负荷不确定性
做电力系统调度优化的朋友应该都体会过那种“模型算出来挺漂亮一到实际就跑偏”的无力感。风电预测给的是100MW实际出力可能只有60MW光伏中午预报大晴天下午一片云飘过来出力直接腰斩负荷就更不用说了居民区和工业区的用电习惯根本不是一条平滑曲线。传统的确定性调度模型把这些预测值当成板上钉钉的常数算出来的机组组合、备用容量、潮流分布一旦遇到风光出力波动轻则经济性大幅恶化重则出现切负荷甚至潮流越限。这也是为什么近些年鲁棒优化在电力系统调度里越来越受重视——它不赌预测值的准确性而是把不确定性“圈”在一个集合里保证集合内所有可能场景下方案都可行。这篇文章想分享的就是一套我实际跑通的计及风、光、负荷不确定性的两阶段鲁棒优化完整实现建模用YALMIP求解器用Gurobi核心算法是CCG列与约束生成中间的双线性项用大M法做线性化。内容会从问题建模一路讲到Matlab代码细节和踩坑经验适合正在做微电网调度、园区综合能源优化、电力系统日前计划这类课题的研究生和工程师参考。1. 为什么确定性调度到了含新能源的系统里就不够用了1.1 确定性模型的“赌徒心态”传统经济调度模型里风电、光伏、负荷都被写成了确定数值。模型长这样目标是最小化总运行成本约束包括功率平衡、机组出力上下限、爬坡速率、网络安全等。求解器给出的结果是一个“在预测值完全准确时才最优”的方案。问题在于新能源的预测误差不是小概率事件而是常态。风电短期预测误差的均方根误差普遍在10%~20%之间极端天气下甚至能到30%以上。如果日前计划按预测值100MW安排火电出力实际风电场发到60MW那40MW的功率缺口就得靠实时市场的备用、储能或者紧急爬坡来补。要是备用容量不够就得限负荷。这种“赌预测准”的思路在传统火电为主的系统里问题不大因为负荷预测相对精准机组调节能力也强。但风光占比提高之后系统的不确定性来源从一个变成了三个预测误差的叠加效应让确定性方案的实际运行成本远高于模型计算值。1.2 鲁棒优化的核心思想不赌但做好准备鲁棒优化的思路完全不同。它不依赖预测值的概率分布而是给每个不确定参数划定一个变化范围构造一个不确定集合Uncertainty Set然后寻找一个对所有集合内场景都可行的方案。说句大白话确定性优化是“我赌明天不下雨所以只带一把伞”鲁棒优化是“我不知道明天下不下雨所以我包里常备雨衣真下大了也不慌”。代价是雨衣占地方、有点重——对应到调度里就是鲁棒方案的经济性会比确定性方案差一些多出来的成本就是为了换取“不管风光怎么波动系统都安全运行”的保障。这个保守度是可以调节的。通过预算约束Budget Constraint我们可以限制“最坏情况下同时偏移到极端的不确定源数量”。比如三个不确定源预算Γ取1意味着优化只保证最多一个源同时取到极端值Γ取3则要求所有源同时取极端值时方案依然可行。Γ越大方案越保守成本越高。1.3 两阶段结构先定下来再调整两阶段鲁棒优化对应的是电力系统调度的实际时序逻辑第一阶段日前决策机组启停、备用容量预留、储能充放电计划等必须提前确定下来。这些决策变量在不确定参数实现之前就要敲定特点是“定了就改不了”。第二阶段实时调整等风光出力、负荷的实际值揭晓后系统在给定第一阶段决策的基础上进行再调度通过调整火电出力、储能充放电、切负荷等手段以最小代价保证功率平衡和网络安全。所以模型的目标函数由两部分构成第一阶段的启动成本和基础运行成本加上第二阶段在最坏场景下的调整成本。外层是第一阶段决策的min问题内层是一个max-min嵌套——先枚举所有可能的不确定场景max再求每个场景下的最小再调度成本min。2. 两阶段鲁棒优化的数学模型拆解2.1 不确定集合的构造与离散化风、光、负荷的不确定集合最常用的是盒式集合加预算约束。先把每个不确定源写成名义值加偏差的形式P_w ∈ [P_w^0 − ΔP_w, P_w^0 ΔP_w] P_v ∈ [P_v^0 − ΔP_v, P_v^0 ΔP_v] P_l ∈ [P_l^0 − ΔP_l, P_l^0 ΔP_l]其中P_w^0是风电预测值ΔP_w是允许的最大偏差三个不确定源的范围共同构成一个超立方体。实际建模时我不建议直接在这个连续区间上处理。原因有两个一是连续区间与第二阶段对偶变量的乘积会构成复杂的双线性项线性化非常麻烦二是从工程角度看真正让调度员头疼的是“某个源出力取到上界还是下界”中间值的影响往往是单调的。所以实操中更常见的做法是引入0-1变量把每个不确定源离散成两个极端状态u_i u_i^0 Δu_i · z_i, z_i ∈ {0, 1}z_i 0表示该不确定源取下界z_i 1表示取上界。再叠加预算约束Σ z_i ≤ ΓΓ是预算参数控制同时取极端的不确定源数量上限。离散化之后不确定集合变成了一个包含有限个极端场景的组合集合这正好为后续CCG算法中“把最坏场景带回主问题”提供了天然的离散场景支持。三个不确定源的具体参数化在算例里常这样设置以下数值仅为示范实际应基于预测数据统计得到不确定源名义值偏差范围离散状态风电出力80 MW±16 MW (20%)64 MW / 96 MW光伏出力50 MW±10 MW (20%)40 MW / 60 MW负荷需求200 MW±20 MW (10%)180 MW / 220 MW2.2 第一阶段决策模型第一阶段决策变量主要包括机组启停状态u_g0-1变量机组启动动作变量用于计算启动成本预留的旋转备用容量储能的充放电状态与SOC初值如果有储能第一阶段模型可以写成min Σ (启动成本 空载成本) ηs.t. 启停逻辑约束、最小启停时间约束、机组出力上下限约束、备用容量约束这里的η是一个辅助变量在CCG迭代中用来逼近第二阶段的期望最坏场景成本这也是两阶段鲁棒优化中“第二阶段成本被第一阶段决策耦合”的标准处理手法。2.3 第二阶段再调度模型给定第一阶段决策x机组启停状态、备用等之后第二阶段模型针对某一具体的不确定场景u求解最小再调度成本min Σ (机组出力调整成本 切负荷惩罚 弃风弃光惩罚)s.t. 功率平衡约束 线路潮流约束直流潮流或交流潮流线性化 机组爬坡速率约束 切负荷量不超过该节点负荷需求 弃风弃光量不超过该时刻风光出力写成紧凑的矩阵形式会更方便后续的对偶推导Q(x, u) min_y b^T y s.t. D y ≥ h − A x − E u y ≥ 0其中y代表第二阶段的连续决策变量机组出力调整量、切负荷量、弃风弃光量等h是常数向量A、D、E是系数矩阵u是不确定向量。2.4 主问题与子问题的耦合结构CCG算法把完整的两阶段鲁棒问题拆成主问题MP和子问题SP交替求解主问题在第k次迭代时min c^T x η s.t. 第一阶段约束 η ≥ b^T y_j, j 1, 2, ..., k D y_j ≥ h − A x − E u_j^*, j 1, 2, ..., k y_j ≥ 0这里u_j^*是前k次迭代中子问题找到的“最坏场景”。每迭代一次主问题就会多一组变量y_j和对应的约束。这个机制通俗点理解就是主问题一开始“乐观”地只考虑已知的场景子问题每次迭代都会找出一个让主问题方案吃瘪的坏场景把它“塞”回主问题逼着主问题改进决策。子问题给定第一阶段决策x^*后Q(x^) max_{u ∈ U} min_y { b^T y | D y ≥ h − A x^− E u, y ≥ 0 }这就是经典的max-min双层结构。外层的max寻找最坏场景内层的min求解该场景下的最小再调度成本。这个双层结构不能直接求解需要借助强对偶定理把它压平成单层优化问题。3. 大M法把双层嵌套压平的关键线性化3.1 强对偶转换子问题中内层的min是一个线性规划满足强对偶条件原问题可行且有界时对偶问题最优值相等。我们对内层min取对偶得到Q(x^) max_{u, λ} λ^T (h − A x^− E u) s.t. D^T λ ≤ b, λ ≥ 0 u ∈ U这是一个纯max问题变量是对偶变量λ和不确定变量u。到这里max-min的双层嵌套变成了单层max可以直接扔给求解器。3.2 双线性项的出现与大M线性化但注意目标函数λ^T (h − A x^* − E u) 里藏着麻烦展开后会有λ^T E u这样的项即对偶变量与不确定变量的乘积。λ是连续变量u经过离散化后是0-1变量λ_i × z_j这类项是双线性项非线性非凸Gurobi和CPLEX都直接处理不了。处理思路是利用大M法把“连续变量 × 0-1变量”的乘积线性化。引入辅助变量v_ij λ_i · z_j用以下一组线性不等式刻画v_ij的行为0 ≤ v_ij ≤ M · z_j λ_i − M(1 − z_j) ≤ v_ij ≤ λ_i M(1 − z_j)这组约束的含义是当z_j 0时第一个约束把v_ij强制压成0当z_j 1时第二个约束把v_ij限制在λ_i的邻域内当M足够大时v_ij就等于λ_i。可以把这个机制理解成“开关”——z_j是开关打开时λ_i通过关闭时乘积归零。项目实践中我曾经绕过离散化、直接用连续u变量做线性化需要引入更复杂的多面体近似实现难度高而且不稳定。相比之下0-1离散化大M法是目前最成熟、最稳的方案特别是与预算约束天然搭配所以我一直推荐这条路。3.3 大M取值的实操经验大M法看起来简单真正翻车最多的地方就是M的取值。M取得太小会错误地切掉可行解导致子问题求出来的目标值偏小CCG迭代出来的结果不对M取得太大会造成数值病态Gurobi求解时会出现“Model may be unreliable”的警告矩阵条件数恶化求解速度大幅下降甚至出现错误的整数解。根据我的实操经验M的取值可以参考以下逻辑先不加入大M线性化约束求解一个松弛的子问题观察对偶变量λ的范围M取λ_i可能最大值的2到5倍一般情况下电力系统调度模型里M取1e3到1e4是比较安全的区间尽量不要超过1e6超过1e6之后Gurobi的数值稳定性会明显变差如果模型做完标幺化处理所有功率除以基准值λ的量级通常在0到几十之间M取500到1000就够用了。这个细节直接决定求解稳定性值得多花几分钟做一次敏感性测试。4. CCG算法主问题和子问题怎么互相“喂料”4.1 为什么选CCG而不是Benders分解两阶段鲁棒优化早期常用Benders分解每次迭代只给主问题加一条割平面约束。CCG的区别在于它不仅加约束还向主问题中直接添加一组完整的第二阶段决策变量和对应的场景约束。每轮迭代主问题的规模会增加但换来的是割平面质量显著提升迭代次数大幅减少。打个比方Benders分解像是一个学生每次考试只被告诉“总成绩不及格”需要自己猜哪道题错了CCG则直接告诉他“这次考试哪道题错了、错在哪、正确答案应该是什么”。虽然每次批改的工作量更大但学生进步快得多。对于两阶段鲁棒优化这种场景CCG已经是事实上的标准算法收敛速度比Benders快一个数量级。4.2 算法伪代码完整流程如下初始化UB ∞LB −∞k 1选取初始场景u_1^*通常取名义值即所有不确定源取预测值求解主问题得到最优解(x_k^, η_k^)更新下界LB max(LB, c^T x_k^* η_k^*)固定x_k^求解子问题得到最坏场景u_{k1}^和最优值Q(x_k^)更新上界UB min(UB, c^T x_k^ Q(x_k^*))如果(UB − LB) / UB ≤ εε通常取0.01或更小则停止迭代输出当前解否则将新场景u_{k1}^*对应的变量y_{k1}和约束η ≥ b^T y_{k1}、D y_{k1} ≥ h − A x − E u_{k1}^*加入主问题k k 1回到第3步这个流程里有个微小但很重要的细节更新UB时每一轮都要取历史所有子问题值中的最小值而不能只取当前值。原因是子问题在给定x_k^*下求出的最坏场景成本本来就是一个可行方案的真实成本上界历史上任何一轮的x_k^*都是一个可行的第一阶段决策对应一个上界。取最小值能得到更紧的上界有助于提前收敛。4.3 收敛判据的工程处理理论上UB和LB会随着迭代逐渐逼近实际代码里我还会加两个辅助收敛条件如果连续三轮迭代中新场景u_k^*与已有场景的差异小于某个阈值比如所有场景变量的汉明距离为0说明子问题已经找不到新信息了提前终止避免无谓的迭代开销设置最大迭代次数比如30次防止数值问题导致的死循环子问题求解时有一个常见的坑如果对偶问题不可行或有界性异常多半是原问题第二阶段约束A、D矩阵写错了或者是M取值过小导致线性化约束过度收紧。排查时建议先把预算Γ设为0退化为确定性场景如果此时子问题正常再逐步增加Γ缩小问题定位范围。5. Matlab YALMIP代码实现与踩坑记录5.1 环境与求解器配置我的环境是Matlab R2022b YALMIP Gurobi 10.0。Gurobi求解MILP的能力业界领先CPLEX也可以但Gurobi在鲁棒优化这类含大量0-1变量和辅助变量的问题上表现更稳。YALMIP的安装不用多说重点说Gurobi参数配置。在CCG迭代中主问题和子问题都是MILP有些参数值得重点关注ops sdpsettings(solver, gurobi, ... verbose, 0, ... gurobi.mipgap, 1e-4, ... gurobi.numericsfocus, 1, ... gurobi.timelimit, 300);numericsfocus设为1或2可以在出现数值警告时显著提升求解稳定性代价是求解时间略微增加在鲁棒优化这种数值敏感性较高的问题里非常值得。timelimit一定要设否则遇到病态模型Gurobi可能跑几小时不出结果。5.2 主问题代码骨架主问题用YALMIP建模的核心逻辑如下% 第一阶段变量 u_start binvar(n_g, 1); % 机组启停 start_action binvar(n_g, 1); % 启动动作 eta sdpvar(1, 1); % 第二阶段成本逼近变量 % 第一阶段约束启停逻辑、最小启停时间、备用容量等 Constraints []; % 历史割集合每轮迭代追加 % u_history 是 n_unc × k 的矩阵存储历史最坏场景 for j 1:k y_j sdpvar(n_y, 1); % 第j轮添加的第二阶段变量 Constraints [Constraints, ... eta b * y_j, ... D * y_j h - A * u_start - E * u_history(:, j), ... y_j 0]; end objective c * u_start eta; optimize(Constraints, objective, ops);这段代码的关键在于u_history(:, j)是已知常数是子问题在历史迭代中求出来的最坏场景。主问题的决策变量是u_start、eta和所有历史添加的y_j。每轮迭代y_j会多一组主问题规模线性增长。5.3 子问题代码骨架子问题经过强对偶大M线性化之后YALMIP代码结构如下% 输入当前第一阶段解 u_start_star lambda sdpvar(n_con, 1); % 对偶变量 z binvar(n_unc, 1); % 不确定变量的0-1表示 % 对偶可行域 Constraints [D * lambda b, lambda 0]; % 不确定集合与预算约束 Constraints [Constraints, sum(z) Gamma]; % 大M线性化v lambda * z v sdpvar(n_con, n_unc); M 1000; % 根据λ的量级调整 for i 1:n_con for j 1:n_unc Constraints [Constraints, ... 0 v(i, j) M * z(j), ... lambda(i) - M * (1 - z(j)) v(i, j) lambda(i) M * (1 - z(j))]; end end % 目标函数λ^T(h - A*x*) - sum(sum(E .* v)) % YALMIP默认求解min因此max问题取负号 Objective lambda * (h - A * u_start_star) - sum(sum(E .* v)); optimize(Constraints, -Objective, ops); % 提取最坏场景 u_star u_nominal delta_u .* value(z); Q_star value(Objective);这代码里有几个关键点值得多说一嘴。大M法的核心是v(i, j)这个辅助变量n_con代表对偶约束数量n_unc代表不确定源数量通常等于3乘以时段数如果按每个时段的风、光、负荷都独立不确定来建模的话。注意这里我没有用双循环以外的向量化写法来简化展示实际项目里建议用repmat或kron做向量化当n_con和n_unc变大时双循环会非常慢。另一个容易踩的坑是目标函数里E .* v的符号。E矩阵是第二阶段约束中u的系数矩阵它可能是负的也可能有正有负取决于原问题的结构。目标函数整理成λ^T(h − A x^* − E u)之后展开项−λ^T E u对应到辅助变量就是−sum(sum(E .* v))。符号反了子问题求出来的“最坏场景”就变成了“最好场景”CCG会完全失效。我第一次实现时就在这里卡了一整天最后是通过对比确定性场景Γ0下的子问题和直接优化结果才发现的。5.4 迭代主循环主循环的写法相对直接Gamma 2; % 预算参数 max_iter 30; epsilon 1e-3; UB inf; LB -inf; k 0; % 初始场景名义值 u_his u_nominal; while k max_iter k k 1; % 求解主问题 optimize(mp_constraints, mp_objective, ops); x_k value(u_start); eta_k value(eta); LB max(LB, value(mp_objective)); % 求解子问题 [u_new, Q_k] solve_subproblem(x_k, Gamma, ops); UB min(UB, c * x_k Q_k); % 收敛判断 if (UB - LB) / abs(UB) epsilon break; end % 向主问题添加新场景 u_his [u_his, u_new]; add_new_cut_to_mp(u_new); % 追加 y_k 变量和相关约束 end注意一个容易忽略的细节LB初始化用−inf但如果主问题在第一轮就不可行说明第一阶段约束写得有问题需要回头检查是不是某个约束把所有机组都禁用了或者备用容量要求过严。实际项目中我遇到过主问题第一轮就不可行的情况排查下来是备用容量约束里留的余量太大超出了所有机组的调节能力总和。5.5 高频问题排查表每次在社区里看到有人问“为什么我的CCG不收敛”或者“为什么结果明显不对”基本都是下面这几个原因之一症状可能原因排查方法子问题目标出现负值强对偶转换后符号弄反退化为确定性场景对比直接优化的结果主问题不可行第一阶段约束过严或冲突去掉备用约束测试逐一放回排查UB和LB长期不收敛M值过小或过大输出λ的范围按3.3节调MGurobi警告数值问题量纲不一致或M过大做标幺化减小M设置numericsfocus子问题耗时爆炸双循环大M约束效率低向量化v约束或减少不确定源个数结果与确定性模型相同预算Γ写成了0或场景初始化有问题检查Gamma赋值和u_his初始化5.6 初始化技巧一个容易被人忽略但非常有效的技巧是用确定性场景的结果初始化CCG。具体做法是先令所有不确定源取名义值求解一个确定性两阶段问题把得到的机组启停方案作为CCG主问题第一轮迭代的初始可行解。这样第一轮子问题就能给出一个质量相当不错的上界能明显减少迭代次数。如果模型中包含储能系统建议把储能的初始SOC也纳入第一阶段决策变量处理。因为SOC的初值直接影响后续所有时段的充放电可行域如果放到第二阶段再决定会破坏“日前决策实时调整”的时序逻辑。6. 从能跑到能用的进阶建议6.1 预算Γ的调参逻辑预算参数Γ是鲁棒优化模型提供的一个“保守度旋钮”。Γ0时模型退化为确定性优化成本最低但对不确定性没有任何抵御能力Γ等于不确定源总数时模型最保守要求所有源同时取极端值时方案依然可行成本最高但安全性最强。实际工程中我会这样给决策者分析跑一组Γ从0到最大值的结果画出成本-保守度曲线。曲线上每个点的斜率可以理解为“再买一单位安全性的边际成本”。大多数情况下Γ从小往大增加时成本会先陡增后趋缓。拐点位置的Γ值往往就是推荐的折中选择——既能覆盖绝大部分风险又不会让成本高到离谱。6.2 与大M法相关的几个数值稳定性经验大M法是最容易引起数值隐患的环节。除了前面说的M取值还有几个经验值得分享所有变量尽量统一量纲功率用MW成本用元不要混着用kW和万元如果模型中同时存在数量级相差1e6以上的系数比如M取1e5而其他约束系数在个位数Gurobi的内置预求解器有时会自行简化模型反而导致线性化约束被错误消去。这种情况可以在Gurobi参数里关掉Presolve的部分功能或者把M缩小大M法的M值即使是同一个模型不同约束上也可以不一样。比如某个对偶变量λ_i理论上限是100那么对应的M取300就足够另一个λ_j上限是1000M可以取3000。分别给值比统一一个大M要稳得多6.3 算例验证的推荐做法验证代码时我建议按以下顺序推进先在修改的6节点系统上测试三个不确定源分别对应风机节点、光伏节点和负荷节点。预测值随便设一组合理的数偏差范围按10%~20%取。先跑Γ0的确定性情形确认主问题和子问题在名义场景下求解正常再逐步增加Γ。然后对比确定性方案和鲁棒方案在极端场景下的表现把确定性方案代入子问题得到它在最坏场景下需要付出的高额再调度成本可能伴随切负荷再把鲁棒方案代入成本会低得多。这个对比是论文里最常用的结果展示也是验证鲁棒优化价值最直观的方式。等6节点算例跑通后再移植到IEEE 14节点或118节点系统。移植时唯一要注意的是网络约束矩阵的维度和稀疏性YALMIP里用稀疏矩阵构造约束比用全矩阵快很多内存占用也小很多。套用这套框架实际做过的项目里6节点系统一般5轮以内收敛118节点系统大概需要10到20轮单轮求解时间视Gurobi性能从几秒到几分钟不等。如果迭代超过30轮还没收敛几乎可以肯定是代码里有bug别怀疑是算法本身的问题。6.4 扩展方向这套两阶段鲁棒优化的框架本身是非常通用的。换一个不确定集合的定义方式比如用数据驱动的多面体集合代替盒式集合换一个目标函数比如加入碳排放约束或者把第二阶段的线性规划换成混合整数规划比如考虑储能充放电状态的实时调整大框架都不需要动只需要改对应的约束和变量定义。更进一步如果想考虑不确定性之间的时空相关性可以考虑用椭球不确定集合或者数据驱动的不确定集合这时大M法的角色会弱化但CCG主-子问题迭代的框架依然成立。最后说一点个人体会。这套方法表面上看是数学建模和代码实现的功夫真正吃透之后最有价值的其实是它对“不确定性成本”的量化能力。当决策者问“鲁棒方案为什么比确定性方案贵这么多”用一条Γ-成本曲线就能讲清楚你愿意花多少钱换来对多少不确定性的保障。这种把工程决策翻译成优化语言的能力是做电力系统优化工作的人最值钱的本事。