两阶段鲁棒优化这套东西我前前后后折腾了小半年从最开始连 min-max-min 结构都绕不清楚到后面把 CCG 算法和大M法在 Matlab 里跑顺中间踩的坑实在太多。最近刚好有朋友问我要这个“计及风、光、负荷不确定性两阶段鲁棒优化”的代码思路索性把整个项目从建模到落地完整总结一遍给想入坑鲁棒优化调度、或者正在被主子问题迭代折磨的同学一份可以直接抄作业的参考。这个项目解决的核心问题很明确风电、光伏、负荷都有预测误差传统确定性调度模型在极端场景下容易翻车。两阶段鲁棒优化就是在第一阶段先做决策第二阶段针对最恶劣的风光出力场景调整运行方案保证不管实际情况怎么变系统都能安全经济运行。实际投产前我是建议先拿小系统把 CCG列与约束生成跑熟再往大系统上迁移否则调试成本会很高。1. 问题建模min-max-min 是怎么拆开的1.1 为什么不能用单阶段期望值模型很多人一开始会问直接用预测值做确定性优化或者用场景法做随机优化不好吗好但不够稳。确定性模型只认一个点如果实际风光出力偏离预测值约束就可能被打破出现弃风、切负荷甚至潮流越限。随机优化要预先给大量场景场景怎么抽、给多少概率本身又是一堆坑而且计算量暴涨。鲁棒优化的思路是“我承认我预测不准但你随便怎么偏我都要有办法兜住”。两阶段鲁棒模型写成标准形式就是min_x c^T x max_{u∈U} min_y d^T y s.t. A x ≥ b G y ≥ h − E x − M u第一阶段决策变量 x 是开机组合、储能充放电状态这类“拍板了就不能反悔”的变量第二阶段变量 y 是各机组出力、储能功率这类可以根据实际场景“事后调整”的变量u 是不确定量代表风电、光伏实际出力和负荷波动。目标函数里第一段是启停成本第二段是运行成本中间那个 max 就是在找“最恶劣的场景”。1.2 第一阶段到底要定什么我在建模时把第一阶段变量定为机组启停状态0-1 变量、储能充放电状态0-1 变量也可以加系统备用容量等长期决策。这里有个关键认知第一阶段决策必须在不确定量实现之前敲定所以它不能依赖 u 的具体取值。第二阶段变量包括常规机组出力 P_g、储能充放电功率 P_ch/P_dis、弃风弃光量、切负荷量等。第二阶段变量是 u 的函数也就是场景变了这些变量可以跟着变。实际代码里每次 CCG 迭代生成一个新的场景 u_k*就要对应引入一组新的 y 变量和约束这一步非常容易写错后面讲算法的时候我会专门说。1.3 不确定性集合的构建是“保守度”的调节旋钮不确定性集合 U 是整个鲁棒优化的灵魂。我用的标准做法是“盒式集合 预算约束”风电、光伏、负荷各自的预测误差都限制在一个区间里同时用预算 Γ 约束总偏差量。以风电为例设预测出力为 P_wt_pred实际出力为 P_wt P_wt_pred Δw其中 Δw 满足−Δw_max ≤ Δw ≤ Δw_max光伏、负荷同理。预算约束写为(|Δw| / Δw_max) (|Δpv| / Δpv_max) (|Δload| / Δload_max) ≤ ΓΓ 越大集合越大代表允许越多的预测误差同时“对着干”结果越保守、成本越高。Γ 0 就退化成确定性模型Γ 取到总不确定维数时就是最极端盒式。我第一次跑的时候把 Γ 设得很大得到的结果保守到机组全开、成本上天后来逐步调整才发现这其实是个很好的灵敏度分析工具。2. 大M法把双线性项“翻译”成 MILP2.1 双线性项藏在几个地方两阶段鲁棒优化本身是一个三层优化问题没法直接用商业求解器求解。CCG 算法的思路是把主问题变成 MILP把子问题变换成单层优化再通过迭代逼近最优解。变换过程中最难啃的就是双线性项。子问题在对偶化之后目标函数里会出现 u^T α 这种形式其中 u 是连续不确定量α 是对偶变量两项都是变量乘积是非凸的。另外如果用 KKT 条件把子问题转为单层优化那互补松弛条件 λ ⊥ (G y − h E x M u) 也是典型的非线性条件。大M法就是把这些非线性条件用“逻辑判断 大常数 M”来线性化。它的思想类似于“一个开关”引入一个布尔变量 z当 z 1 时执行某个约束z 0 时放松该约束用足够大的 M 来“关掉”不需要生效的约束。2.2 三种最常用的线性化模板第一种逻辑或or约束。比如要表达“要么支路潮流不超过上限要么切负荷”可以写f(x) ≤ Limit M·z f(x) ≥ −Limit − M·zz 为 0 时约束生效z 为 1 时约束被 M 放松失去限制力。第二种连续变量与 0-1 变量相乘。比如储能充放电功率 P 与充放电状态 z 相乘0 ≤ P_ch ≤ z_ch · P_ch_max这其实是“变量有界、由状态决定”的写法不需要真的把 z·P 拆开但要写成两个约束P_ch ≤ M·z_ch P_ch ≤ P_ch_max P_ch ≥ 0这里是“大M”最常见的应用场景也是初学者最容易漏写第三个约束的地方。第三种互补松弛条件。KKT 条件里要求 λ ≥ 0g(x) ≥ 0且 λ·g(x) 0可以用下面这组约束替代λ ≤ M·z g(x) ≤ M·(1 − z) λ ≥ 0, g(x) ≥ 0, z ∈ {0,1}这套模板我项目里用了无数遍必须记牢。2.3 M 值怎么选选错有多致命大M最忌讳的就是“随便给个 1e6 完事”。M 太小会错误地截断可行域导致漏掉最优解甚至得到不可行问题M 太大求解器数值稳定性崩溃出现各种匪夷所思的结果。我的经验是先估算模型里的量级。比如功率平衡约束里火电最大出力可能是 300 MW机组台数 5 台那 M 取 1000 就够如果涉及成本项成本系数可能在几十到几百考虑总成本后 M 取 1e6 也合理。你可以算一下各约束里变量的系数范围再按“比最大可能取值大 10~100 倍”的标准来取。实操中更稳妥的做法是先用一个较大的 M 求解观察结果里有没有变量被 M 过度放松的痕迹比如某变量取值明显超出物理边界再逐步缩小 M 直到结果稳定。这个方法虽然笨但比拍脑袋强得多。3. CCG 算法主子问题迭代的核心流程3.1 主问题 MP先“猜”几个最坏场景CCG 算法的核心思路是把无穷多个不确定场景 U 离散成有限集合每轮迭代往集合里加入一个“当前最恶劣”的场景反复逼近真实最坏情况。主问题形式为min_{x, y_j, θ} c^T x θ s.t. A x ≥ b θ ≥ d^T y_j, 对 j 1, 2, ..., k G y_j ≥ h − E x − M u_j*, 对 j 1, 2, ..., k其中 u_1*, u_2*, ..., u_k* 是已经生成的场景。每加入一个场景就要引入一组新的第二阶段变量 y_j 和对应约束。θ 代替了原本的 max 部分相当于把“最恶劣场景的运行成本”用一个标量 θ 来逼近。主问题的解给出第一阶段决策 x* 和一个下界 LB c^T x* θ*因为我们只用了一个场景子集真实最优成本一定不低于这个值。3.2 子问题 SP在固定决策下找最恶劣场景子问题是固定第一阶段决策 x*求解不确定场景下的最恶劣运行成本Q(x*) max_{u∈U} min_{y} d^T y s.t. G y ≥ h − E x* − M u这是一个 max-min 两层问题。标准解法是取内层 min 的对偶将问题转为单层 max 问题Q(x*) max_{u, α} (h − E x* − M u)^T α s.t. G^T α ≤ d, α ≥ 0, u ∈ U问题来了目标函数里出现 u^T (M^T α) 这个双线性项。处理方式有两种一是用大M法把这个双线性项线性化得到一个 MILP二是枚举不确定集合 U 的极点因为线性目标函数的最优解一定在极点取到。实际代码里当 U 是盒式集合时不确定变量的极值组合数量是 2^NN 为不确定变量维数N 大时枚举不现实所以主流还是大M线性化。子问题求解结果给出当前 x* 下的最恶劣场景 u_new* 和运行成本 Q(x*)上界更新为 UB min{UB, c^T x* Q(x*)}。3.3 迭代终止条件与初始化技巧算法在上下界间隙小于阈值时终止(UB − LB) / UB ≤ εε 我习惯取 1e-3 或 1e-4太小会导致迭代次数增加太大则结果不可信。初始化时通常会取预测场景即 u 取预测值或某个极端场景作为第一个场景。我实际测试下来初始场景取预测场景可以让前几轮迭代更平稳但整体迭代次数差别不大。每次迭代的流程求解 MP得到 x_k* 和 LB固定 x_k*求解 SP得到 Q(x_k*) 和 u_k1*更新 UB min{UB, c^T x_k* Q(x_k*)}如果 (UB − LB) / UB ≤ ε停止否则将 u_k1* 加入 MP 的场景集合k k1回到第 1 步。4. Matlab 代码实现建模、循环与性能调优4.1 建模层选择别自己写单纯形法Matlab 里求解这种问题我强烈建议用 YALMIP 或 CVX 做建模层再加 Gurobi 或 Cplex 做底层求解器。YALMIP 对 0-1 变量、大M法、双线性项的建模支持很好binmodel 命令还能自动做部分线性化。安装组合我推荐Matlab YALMIP Gurobi学术版免费。YALMIP 配置 Gurobi 只需要几分钟关键是把 Gurobi 的安装目录添加到 Matlab 路径里。建议直接下载 Gurobi 的 Matlab 接口包里面有现成的 setup 脚本。4.2 代码框架主函数、MP、SP 三块结构我的代码分成三个部分main_run.m数据初始化、参数设置、CCG 主循环build_mp.m构建主问题 MILP 并求解build_sp.m构建子问题 MILP 并求解。主循环的伪代码如下% 初始化 LB -1e6; UB 1e6; k 1; u_list{1} u_pred; % 初始场景取预测值 eps 1e-3; while (UB - LB) / abs(UB) eps % 1. 求解主问题 [x_k, theta_k, LB_new] build_mp(u_list); LB max(LB, LB_new); % 2. 求解子问题 [Q_k, u_new] build_sp(x_k); UB min(UB, c*x_k Q_k); % 3. 判断是否加入新场景 if (UB - LB) / abs(UB) eps k k 1; u_list{k} u_new; end end4.3 主问题 MP 的 YALMIP 写法主问题中每加入一个场景要动态增加一组变量和约束。YALMIP 支持在循环里累加约束用 cell 数组保存变量比较高效x_start binvar(n_gen, 1, full); % 机组启停 theta sdpvar(1, 1); % 最恶劣场景成本 Constraints []; for j 1:length(u_list) y{j} sdpvar(n_y, 1); % 每个场景一组连续变量 Constraints [Constraints, theta d*y{j}]; Constraints [Constraints, G*y{j} h - E*x_start - M*u_list{j}]; % 其他场景约束... end Constraints [Constraints, A*x_start b];运行时我会用 optimizer 命令预编译问题避免每轮迭代重复建模能省不少时间。4.4 子问题 SP 的 KKT 与大M线性化子问题我采用“先对偶、再大M”的方案。内层 min 是线性规划写出对偶后目标函数为max_{u, α} (h − E x* − M u)^T α约束 G^T α ≤ d, α ≥ 0, u ∈ U。其中 u^T α 是双线性项。大M法处理思路为引入辅助布尔变量 z 表示 u 是否取到上边界把连续变量 u 写成u u_min (u_max − u_min)·z但这里 u 是连续向量而 z 是 0-1乘积项仍然存在。更实用的做法是利用“线性规划最优解在极点”的性质直接枚举盒式集合的极点2^N 个或者用 YALMIP 的 binmodel 自动处理双线性项。YALMIP 会自动引入大M约束但需要保证 u 和 α 都有界。α 的有界性一般可以通过约束 G^T α ≤ d 和 α ≥ 0 在工程参数下自然保证实际调试时如果求解器警告无界我会额外加一个 α ≤ 1e4 的约束。另一个常用思路是用 KKT 条件将 SP 一步转化为单层 MILP。内层问题 min_y d^T y 的 KKT 条件是G y ≥ h − E x* − M u λ ≥ 0 G^T λ ≤ d λ ⊙ (G y − h E x* M u) 0其中 λ 是内层问题的对偶变量按 ≥ 约束的标准形式。最后一条互补松弛条件用大M线性化% 互补松弛条件线性化 for i 1:length(lambda) z{i} binvar(1, 1); Constraints [Constraints, lambda(i) M_big * z{i}]; Constraints [Constraints, slack(i) M_big * (1 - z{i})]; end注意这里要额外定义松弛变量 slack G y − h E x* M u ≥ 0。5. 常见问题与调试经验5.1 迭代上下界不收敛或震荡我最开始跑 CCG 时遇到最常见的问题就是上下界不闭合曲线像锯齿一样跳。排查思路第一检查 SP 是否真的求解到了全局最优。如果 SP 是非凸的比如双线性项没处理干净得到的“最恶劣场景”可能不是真正的最高成本场景导致 UB 偏低。第二检查 MP 中每个新增场景的变量是否独立。如果多个场景共用一组 y 变量约束会互相干扰下界会失真。第三检查初始场景。把 u_list{1} 设为预测值通常比较平稳但如果你特别关注极端场景也可以从某个角点场景启动。5.2 对偶问题方向搞反对偶这步容易出错特别是约束里有等式和不等式混合的时候。我的经验是把内层问题先化成标准形式min d^T y s.t. G y ≥ h − E x* − M u y ≥ 0然后对偶得到max λ^T (h − E x* − M u) s.t. G^T λ ≤ d, λ ≥ 0如果内层约束有等式对偶变量就是自由变量无正负约束这点特别容易写错。建议每写完一个典型对偶就往 YALMIP 里塞一个随机小算例做验证数值对比能快速暴露问题。5.3 求解器提示无界或数值病态Gurobi 报 unbounded 绝大多数是 M 值太大或边界条件缺失。我项目里曾经给 M 取到 1e8结果求解器直接警告数值精度丢失约束被畸形的 M 扭曲得到一堆 NaN。后来我把所有约束的系数做了归一化M 值按各类变量最大可能取值估算问题立刻消失。补充一点如果模型里同时有特别大和特别小的系数例如成本和功率相差 6 个数量级建议把目标函数分成“单位一致”的项或者对成本项做缩放比如全部除以 1000能显著降低求解器病态概率。5.4 不确定集合参数怎么调才合理Γ 预算值需要根据实际系统的预测精度来定。我做过一组实验Γ 0 时成本最低但风光全部取预测值Γ 从 0 增加到 1/3 总不确定维度时成本上升 5% 左右但极端场景下的失负荷量显著改善Γ 再往上就只剩下成本增加、收益甚微。建议你把这个“成本-保守度”曲线跑出来由调度人员根据风险偏好选点。5.5 一个容易被忽略的细节子问题可能不可行在某些第一阶段决策下第二阶段调度可能无解。这时需要给 SP 加松弛变量如虚拟切负荷并在目标函数里加一个足够大的惩罚系数。否则 CCG 会在某轮迭代突然报“infeasible”整个流程直接崩掉。我在负荷波动较大的算例里加过松弛量效果很好代价是要多调试一个惩罚系数的量级。问题现象可能原因解决方案UB 反复跳、不下降SP 双线性项处理错误检查对偶推导用 binmodel 或枚举验证最值LB 一直不增长MP 新增场景约束没生效检查每轮新增 y 变量和 u 索引是否对应求解器报 infeasible子问题缺松弛变量加虚拟切负荷/弃风惩罚项结果出现明显不合理值M 值太大或太小估算量级M 取最大量级 10~100 倍迭代次数很多才收敛Γ 过大或容差过严调整不定集合预算适当放宽 ε最后分享一个我个人的实战经验刚上手别急着完整复现 118 节点系统。先用 6 节点小系统风电光伏各一台负荷两个节点手写每一轮迭代的 MP、SP 输出把上下界变化曲线画出来。看到曲线一步一步收敛你对 CCG 的理解才算真的建立起来。然后再上大系统时你才知道哪些问题是模型本身的问题哪些是求解器的数值问题。大M法在项目里其实只占很小一段代码却承担了整个“非线性转化为线性”的关键任务它的调试成本远高于它占的代码量。建议你专门封装一个linearize_with_bigm.m函数把三种线性化模板都收进去后续换系统、换约束时直接复用省心很多。