开年第一篇不写花架子直接聊干货。最近在做一个计及风电、光伏、负荷不确定性的两阶段鲁棒经济调度项目模型本身不算特别新鲜但把鲁棒优化、大M法、CCG算法这几样东西串在一起用Matlab实现中间踩的坑、绕的弯确实值得拿出来复盘一遍。这篇文章从问题建模开始讲到CCG求解框架、大M法处理双线性项的具体操作最后给出可直接参考的Matlab实现结构和调试心得。适合正在做电力系统调度、微电网优化、储能配置这类方向并且对鲁棒优化和分解算法有一定理论基础、但卡在代码实现上的朋友。1. 问题建模为什么两阶段鲁棒在哪里1.1 第一阶段确定性决策与第二阶段调整我们面对的问题是标准的日前调度加实时调整场景。第一阶段决策指的是机组启停、日前出力和备用容量安排等必须在风光和负荷实际值观测之前敲定约束里只包含确定性参数和常规机组的技术约束。第二阶段是实时调整阶段当风光实际出力、负荷实测值逐渐明确之后调度系统可以在一定范围内调整机组出力、切负荷或弃风弃光以弥补日前计划与实际情况之间的偏差。这一段的决策是在不确定性揭示之后做出的因此也叫“等待与看到”决策。两阶段鲁棒优化的核心结构数学上写成[ \min_{x} (c^T x \max_{u \in U} \min_{y \in F(x,u)} d^T y) ]其中 (x) 是一阶段决策变量(u) 是风光出力和负荷的不确定参数(y) 是第二阶段调整变量。最外层最小化表示总成本最小中间层的 max 表示在所有可能的情况下取最恶劣的场景内层 min 则表示在给定最恶劣场景下尽可能用最经济的调整手段去补救。三层结构叠在一起就是鲁棒优化的标准姿态我不赌某个特定场景而是要求作出一个“即使碰上最坏情况也不至于崩溃”的方案。1.2 风光、负荷不确定集合的建模方式不确定集合的选择直接决定了模型保守程度和求解复杂度。最常用的三种集合类型数学表达优点缺点箱式集合( |u - \bar{u}|_\infty \le \alpha )建模简单最保守盒式/区间集合( u \in [\underline{u}, \bar{u}] )直观仅限独立变量预算约束集合( |u - \bar{u}|1 \le \Gamma, |u - \bar{u}|\infty \le \alpha )可调保守度需调参数我实际用的是预算约束budget集合风电、光伏和负荷分别加了自己的不确定性预算参数。它的物理意义很直观不可能所有风电场的预测误差同时达到最坏值预算参数限制了“同时偏移”的场站个数总和。(\Gamma) 取0就是确定性模型(\Gamma) 取满就是最保守的集合通过调节它可以画出一条鲁棒性和经济性的帕累托曲线这部分在论文里非常受欢迎。这里有个细节我最开始没处理好——风电和光伏的不确定资源量级不一样如果共用一个预算参数会出现光伏的偏移量基本不影响总预算的情况。实践中的处理方式是对各类资源分别设 (\Gamma)再在总预算上做一个归一化加权这样物理意义更清晰也方便调整保守程度。1.3 为什么不用随机规划或场景法聊到这里肯定有做随机优化的朋友要问直接用蒙特卡洛抽场景不好吗我的观点两阶段鲁棒优化和随机规划解决的是不同层面的问题。随机规划需要知道不确定参数的精确概率分布场景法也依赖大量采样来逼近分布计算量和场景数的平方级增长很容易让模型变成一头推不动的巨兽。更关键的是概率分布本身也存在估计误差万一真实分布跟模型分布偏差较大随机规划的解很容易出现“理论上最优、实际里翻车”的情况。鲁棒优化不依赖精确概率分布只需要确定不确定参数的上下界和预算水平虽然解会偏保守但它的可行性保障是非常硬的——只要实际值落在我定义的不确定集合内方案一定可行。这种“牺牲一点经济性换来极高的可靠性”的特点在电力系统这种对安全性极其敏感的行业里非常受用。2. CCG算法精讲迭代生成极端场景2.1 从Benders分解到CCG的演化逻辑处理两阶段鲁棒优化这类大规模问题直接用商业求解器单层求解几乎不可能尤其当网络规模扩展到时节点配电网决策变量动辄几十万个。分解算法是必然选择。Benders分解处理两阶段问题的思路是把第二阶段问题对偶得到一组割平面反馈给主问题。这思路成熟但在鲁棒双层结构中会碰到一个麻烦——第二阶段本身里面嵌套了一层max-min对偶之后影子变量的几何维度很高割平面的收敛速度往往不理想遇到整数变量更是直接卡死。CCG算法列与约束生成是2013年Zeng Bo和Zhao Long提出的改进方案核心思想非常直接不在主问题中添加割平面而是把第二阶段决策变量和约束“显式”地加入主问题每次迭代加入一个对应最恶劣场景的新“列”。这样做的好处有三点主问题的规模是可控增长的每一步只增加一组变量和约束上下界的收敛速度快通常十几轮就可以稳定支持第二阶段含整数变量适用范围比Benders广得多。2.2 主问题与子问题的分工CCG把原问题拆成两部分主问题MP包含一阶段决策变量和若干个已生成的极端场景对应的二阶段变量是一个相对紧凑的混合整数线性规划。[ \min_{x, y_k} ; c^T x \alpha ][ s.t. ; A x \le b, \quad D x E y_k \le f G u_k^*, \quad k 1, 2, ..., K ]其中 (u_k^*) 是迭代过程中子问题返回的极端场景(K) 是当前迭代轮数。主问题求出的目标值是原问题的一个下界对应原问题的 min 部分。子问题SP固定主问题得到的一阶段解 (x^*) 后求解第二阶段的 max-min 问题找到当前方案下最恶劣的不确定性场景[ Q(x^*) \max_{u \in U} \min_{y} ; d^T y ][ s.t. ; D x^* E y \le f G u ]子问题求得的目标值加上一阶段成本构成原问题的一个上界。然后就是两边轮着来主问题给子问题传一阶段解子问题把最恶劣场景传回主问题直到上下界间隙小于预设容忍度。迭代过程中我观察到前几轮上下界收敛极快往往第三轮就能落下80%的间隙后面几轮属于精调阶段。2.3 收敛判据与迭代终止条件CCG的迭代终止条件用相对间隙而不是绝对间隙否则量级大的模型很难收敛。我常用[ \frac{UB - LB}{UB} \le \epsilon ](\epsilon) 取1%或0.5%。有一点要特别注意如果问题本身不可行或求解器数值不稳定可能会出现上下界震荡不收敛的情况这时候不要盲目调大 (\epsilon)应该回去检查子问题的对偶间隙和主问题的约束松弛。2.4 关于Benders与CCG的决策参考对比维度Benders分解CCG算法收敛速度慢割平面信息弱快加入完整场景约束二阶整数变量不支持支持实现复杂度较低中等主问题规模增长只加一行割加一组变量约束适用场景纯线性大规划两阶段鲁棒优化/随机优化一个常见误区是下意识把所有两阶段问题都套CCG。如果二阶段问题本身没有不确定性或不确定性只是对等式右侧常数进行扰动Benders可能是更顺手的工具。CCG最大的优势在于它天然把“场景”这个语义单元引入主问题和鲁棒优化中“极端场景”的概念一拍即合。3. 大M法处理子问题对偶的绊脚石3.1 双线性项从哪来子问题内部是一个 max-min 结构。处理它的标准套路是内层取对偶把 min 问题变成 max 问题于是得到 max-max 即一个 max 问题。但在做对偶时会因为不确定参数与对偶变量相乘而产生双线性项 (u^T \pi)。以风电不确定集合为例令风电出力 (p_w \bar{p}_w \Delta p_w \cdot z_w)其中 (z_w \in [-1, 1])。对偶之后目标里会出现 (\pi \cdot \Delta p_w \cdot z_w)(\pi) 是对偶变量(z_w) 是不确定缩放因子两项相乘就是双线性。这问题数学上是非凸的不能直接丢给Gurobi或CPLEX。主流处理方案有两种大M法把 (z_w) 替换为0-1离散变量组合通过引入M值线性化双线性项KKT条件法用内层问题的KKT条件替换但引入互补松弛条件需要 SOS1 约束或大M法配合本质上复杂度类似。MATLAB环境下我用的是大M法将不确定集合离散化配合YALMIP建模效果稳定且容易调试。3.2 大M的引入与线性化过程以电流运行的线性化为例假设子问题内层对偶问题中包含乘积项 (z \pi)其中 (z \in [0,1]) 是连续变量。引入二进制变量 (v) 和大M常数令 (w z\pi)等价约束为[ 0 \le z - v \le M(1 - \beta) ][ 0 \le \pi - w \le M(1 - \beta) ][ 0 \le w \le M \beta ][ \beta \in {0, 1} ]这样就把非凸的双线性项转化成了混合整数线性约束。更常见的简化版是离散化 (z) 取值到若干个档位每个档位对应一个二进制变量再让对偶变量与二进制变量的乘积通过大M松弛为线性约束。3.3 大M取值过大会出问题大M的取值经验是“够大就行”这句话害死人M取得过大直接导致数值病态求解器会出现线性约束失效、误判可行性的问题。我踩过的坑是M取 (10^6)Gurobi报numerical trouble松弛解大量违约。实际做法是通过一次预求解给每个双线性项单独计算M值。以大M法的典型约束 (w \le M\beta) 为例M的取值上限可以通过对偶变量 (\pi) 的理论边界来确定。对偶变量的经济含义是资源的影子价格在机组出力上下限约束中对偶变量的量级基本与边际成本同级别因此M取该系统最大边际成本的5~10倍即可。直接用一个全局统一的10万级别的M本质上是在惩罚求解器的数值稳定性。M取值经验速查先解一次不含双线性项的子问题记录对偶变量最大绝对值取该值的5~10倍作为初始M如果后续求解出现不可行逐步放大1.5倍而不要一步跳一个数量级每次求解后检查关键线性化约束的松弛余量若余量过大说明M不合适。3.4 另一种思路对偶线性化的替代方案如果大M法调M调到怀疑人生可以试试外逼近法把非线性项用一阶泰勒展开迭代逼近。但这个方法在主问题每次都会引入新的辅助变量实现复杂度不低。另一种思路是用强对偶直接消内层。如果内层是线性规划且满足Slater条件那么 min 和它的对偶 max 在最优值处相等可以把内层对偶直接“嵌入”外层得到一个单层max问题目标函数里的双线性项仍然需要处理。本质上绕不开但大M法在当前商业求解器框架下已经是工程上最成熟的方案。4. Matlab代码实现从建模到迭代求解4.1 系统参数输入与数据结构设计Matlab实现的第一步是看数据怎么组织。我的经验是把系统参数、不确定集合参数、求解器参数分开三个struct存放后续调试会省心很多。% 系统参数定义示例 sys.NG 6; % 常规机组数 sys.NW 2; % 风电场数 sys.NP 1; % 光伏场数 sys.T 24; % 调度时段 % 常规机组参数出力和爬坡约束生成矩阵 % 注意这里只是示意真实数据建议用Excel读入再初始化 sys.pmin [100; 80; 60; 40; 30; 20]; sys.pmax [300; 250; 200; 150; 120; 100]; sys.ramp [100; 80; 60; 40; 30; 20]; % 不确定性预算参数 unc.alpha_w 0.2; % 风电归一化波动范围 unc.alpha_p 0.15; % 光伏归一化波动范围 unc.alpha_l 0.05; % 负荷归一化波动范围 unc.Gamma_w 4; % 风电预算 unc.Gamma_p 2; % 光伏预算一个容易出错的点是起步阶段决策变量的下标映射。常规机组出力、风电出力、光伏出力、负荷和各节点相角全部要写成统一的索引映射函数否则在构建主问题的场景变量时会乱套。我习惯直接定义idx (t, g) (t-1)*NG g;这样的匿名函数做索引转换简单粗暴但有效。4.2 主问题建模与YALMIP接口Matlab环境下我推荐用YALMIP做建模层求解器底层接Gurobi或CPLEX。YALMIP的symbolic建模方式在调试时非常直观变量约束一眼就能看明白。举个例子主问题的核心约束构建% 主问题建模框架 % x表示一阶段变量y_k表示第k个场景对应的二阶段变量 % 需要定义K个场景对应的变量单元 ops sdpsettings(solver, gurobi, verbose, 0); MP.x sdpvar(sys.NG, sys.T); % 一阶段机组出力 MP.alpha sdpvar(1, 1); % 辅助变量表示二阶段成本 % 场景循环加入变量和约束 for k 1:K MP.y{k} sdpvar(sys.NG, sys.T); % 场景k下的调整出力 MP.theta{k} sdpvar(sys.NB, sys.T); % 场景k下的相角 MP.ls{k} sdpvar(sys.NL, sys.T); % 场景k下的切负荷量 % 功率平衡约束 Constraints [Constraints, ... sum(MP.x MP.y{k}, 1) ... sum(unc.w_real{k}, 1) sum(unc.p_real{k}, 1) ... - MP.ls{k} load_real{k}]; end这里 YALMIP 做变量管理的好处是你不需要手动拼接大规模稀疏矩阵约束就按业务逻辑一行行写由YALMIP去翻译成求解器的标准形式。4.3 子问题构建与对偶化子问题在给定一阶段解 (x^*) 后构建对偶化处理是我封装成一个独立函数solve_subproblem(x_fixed)来做的。function [obj_sp, u_best] solve_subproblem(x_fixed) % 定义对偶变量 pi_balance sdpvar(sys.NB, sys.T); % 功率平衡对偶 pi_gen_lo sdpvar(sys.NG, sys.T); % 出力下界对偶 pi_gen_hi sdpvar(sys.NG, sys.T); % 出力上界对偶 % 双线性项线性化引入z的离散化和大M z_w binvar(sys.NW, sys.T); % 极端场景选择 w_aux sdpvar(sys.NW, sys.T); % 辅助变量w z*pi % 大M约束 M_val 1000; % 具体值通过预求解校准 for i 1:sys.NW for t 1:sys.T Constraints [Constraints, ... 0 w_aux(i,t) M_val * z_w(i,t), ... pi_balance(unc.map(i),t) - M_val*(1-z_w(i,t)) w_aux(i,t), ... w_aux(i,t) pi_balance(unc.map(i),t)]; end end % 目标函数最大化对偶目标 Obj_SP - sum(sum(pi_gen_min .* pmin)) ...; optimize(Constraints, -Obj_SP, ops); end有用的一个细节子问题的优化方向。因为是max-min结构转对偶后变成纯max问题在Matlab里调用optimize时要对目标取负号否则会默认求解最小化。这是我第一次犯的低级错误结果子问题一直在求下界整个迭代逻辑完全乱套。4.4 CCG主循环迭代求解主循环是整个程序的核心骨架逻辑上并不复杂但每一行都关键LB -inf; UB inf; K 0; constraints_MP []; scenarios {}; x_init init_solution(); % 怎么获取初始可行解后面专门讲 while (UB - LB)/abs(UB) tol iter max_iter % 1. 解主问题得到一阶段解和LB [x_opt, alpha_val] solve_MP(constraints_MP, scenarios); LB c*x_opt alpha_val; % 2. 固定x_opt解子问题得到极端场景和上界增量 [obj_sp, u_new] solve_subproblem(x_opt); UB min(UB, c*x_opt obj_sp); % 3. 如果间隙不收敛将新场景加入主问题 if (UB - LB)/abs(UB) tol K K 1; scenarios{K} u_new; % 在主问题中新增一组二阶段变量和约束 constraints_MP [constraints_MP, new_scenario_constraints(u_new)]; end fprintf(iter%d, LB%.2f, UB%.2f, gap%.4f%%\n, ... iter, LB, UB, (UB-LB)/abs(UB)*100); iter iter 1; end这段逻辑看起来简单但每一步之间都有坑。主问题返回的alpha是辅助变量的值对应的是目前所有已生成场景下二阶段成本的最大值这个是LB计算的核心。子问题返回的obj_sp是给定当前一阶段方案后的最恶劣二阶段成本它和一阶段成本加在一起应该和UB比较。新手最容易出错的点就是用错了上下界的构成方式导致gap怎么迭代都不收敛。4.5 初始可行解怎么给CCG的每一步都需要一个当前一阶段解来启动子问题。第一轮迭代时这个解怎么来是有讲究的。最简单的办法是解一个确定性场景下的模型用预测值替代所有不确定参数得到初始解。但有一个风险这个解在不确定极端情况下可能不可行子问题会报infeasible。CCG算法理论上是允许这种情况的但工程实现上会让代码的鲁棒性变得很难维护。更稳妥的做法是给主问题增加松弛变量让初始迭代不会因为不可行而崩掉% 在主问题中增加松弛变量 s_slack sdpvar(sys.NB, sys.T); Constraints [Constraints, sum(MP.x,1) ... s_slack load_pred]; Constraints [Constraints, s_slack 0]; % 目标函数中加惩罚项 Objective c*MP.x MP.alpha 10000 * sum(sum(s_slack));等迭代稳定后逐步把惩罚系数调大最终实现无松弛的严格可行解。这不是纯CCG该干的但工程上很多产品化的代码都这样干是实用主义的选择。5. 调试实录与常见问题排查5.1 子问题对偶不成立CCG中这一步最容易出的问题是内层min问题不满足强对偶条件。模型中有离散变量、或者约束矩阵不满足线性规划约束规格如存在冗余约束导致Slater条件不满足对偶就会出问题。我当时排查的过程是固定一阶段解后先单独求解原min问题记录目标值再求解对偶max问题比较两者是否一致。如果不一致逐条检查内层约束最终发现是一处潮流约束写重复了导致对偶间隙非零。处理办法是删掉冗余约束同时用预求解器做一步模型化简。5.2 上下界不收敛或间隙震荡这个问题表现为LB和UB在几轮迭代后始终在一个区间内抖动降不下去。可能原因按概率排序现象常见原因解决办法gap卡在5%左右不动大M取值偏小线性化松弛过紧倍增M重新求解gap在第3轮开始反弹子问题返回的场景不是全局最恶劣检查子问题目标函数方向gap一直很大且不下降主问题缺少部分场景约束检查场景索引是否错位数值警告numerical troubleM值过大逐项校准M避免全局统一大M我实际卡得最久的是第二种情况子问题目标函数符号写反了返回的“最恶劣”场景变成了“最温和”场景主问题不断加入温和场景上界自然收不下去。这个错误的隐蔽性在于求解器不会报错日志看起来也正常只能靠人工对比每个场景的目标值来判断。5.3 预算参数Gamma影响分析预算参数 (\Gamma) 对结果的影响非常直观推荐大家都跑一遍灵敏度分析。在风电、光伏和负荷三种不确定性共存的情况下我测试的结果是(\Gamma) 较小时如2以内鲁棒解和确定性解的成本差距在1%~3%符合理论推断(\Gamma) 增大到中段时成本缓慢上升但约束满足的极端场景覆盖比例显著增高(\Gamma) 接近上限时成本跳跃式上升这是预算约束在“最坏场景”附近变得异常严苛导致的。我的建议是实际工程应用中(\Gamma) 的大小应该基于历史数据的最大同时偏差次数统计来确定而不是为了论文好看拍脑袋取个大值。否则鲁棒优化耗散的经济性可能会让项目的投资回报率很难看。5.4 Matlab具体函数的几个坑sdpvar定义三维变量时YALMIP打包后约束构建速度会明显变慢建议用循环逐个时段构建约束虽然代码看着啰嗦但效率更高。Gurobi在Linux和Windows上对MIP的默认参数有些许差异如果换平台后结果对不上先检查相关参数是否一致。Matlab R2020b之前对YALMIP的兼容性有一些已知问题建议至少R2021a以上版本。大模型求解时会内存溢出记得加gurobi的NodeFileStart参数让节点缓存落盘。optimize返回的solveyalmipproblem状态信息不要忽略YALMIP实际上提供了很多调试信息只是太多人直接忽略info字段。5.5 调试节奏与日志设计最后分享一下调试节奏的经验。不要一上来就跑完整24时段、31节点的大规模模型几乎必翻车。我的习惯是先用一个3节点、1台机组、1个风电场的迷你模型验证算法逻辑把主问题变量打印出来手算目标函数和自己核对一遍再升级到6节点、3机组的小系统验证不确定集合预算参数的效果最后才是完整系统。主循环的每次迭代日志我也会打出详细的信息包括主问题求解时间、子问题求解时间、当前场景的目标值、新增约束数量。这些信息在论文里写“算法收敛速度比较”的时候都是最直接的数据支撑。6. 从代码到论文结果分析与扩展方向6.1 结果可视化的几个推荐角度分析结果时下面这些图是做鲁棒优化项目必出的不同 (\Gamma) 值下总成本变化曲线展示鲁棒性和经济性的权衡极端场景下的机组出力时序图展示第二阶段调整的合理性不确定性集合中选取的极端场景分布图验证场景选取逻辑常规机组爬坡压力对比确定性方案和鲁棒方案在极端场景下的表现对比。Matlab的绘图体系在这里很好用我在代码末尾封装了一个plot_results.m把所有关键图自动输出成高分辨率PNG直接可用在论文或汇报材料中。6.2 模型的扩展方向两阶段鲁棒优化的框架非常通用扩展方向也很多加入储能系统后第二阶段变量会包含储能的充放电决策不确定性集合需要同时考虑荷电状态的时间耦合约束考虑需求响应时负荷侧的弹性资源也可以描述成第二阶段调整手段交通网和电网耦合的场合电动汽车集群不确定性来源又多了一类此时CCG的求解框架依然适用只是子问题规模会进一步膨胀从两阶段推广到多阶段鲁棒优化rolling horizon这时要引入模型预测控制思想CCG的循环结构依然可以复用只是约束添加和数据管理变得更复杂。6.3 关于求解规模的现实判断最后说点很多人关心但文献不太提的问题Matlab YALMIP Gurobi这一套到底能撑起多大规模的问题我自己的实测数据供参考30节点以下配电网系统、几十台机组规模的调度模型这个组合可以比较轻松地处理CCG迭代时间通常在几十秒内上百节点、大规模机组组合问题时主问题规模增长会很快这时建议考虑把主问题和子问题的求解器分开设置主问题用Gurobi跑MIP子问题用Cplex或Gurobi跑LP同时加上热启动机制如果模型规模再大一个数量级纯Matlab的建模层会成为瓶颈这时需要考虑直接用Python Pyomo/JuMPJulia等更底层的建模框架甚至考虑将CCG算法写成并行版本因为不同场景的可行性检验天然是并行化的。从工程落地角度说我个人的感受是CCG算法在前几轮就能给出质量很高的解真正限制效率瓶颈的往往不是算法本身而是建模层和求解器之间的数据交互效率。所以代码实现的精细程度有时候比算法选型更影响最终效果。每次跑完这个程序看到它在十几轮迭代内把上下界逐步收敛到千分之一的间隙还是会觉得这套数学机制确实精巧。但也提醒自己模型不收敛、数值病态、场景选取出错每一个坑都真实发生过。希望这篇拆解能帮后来者少走一些弯路。最后再分享一个小建议不妨把CCG实现封装成一个通用的求解器类输入不需要限定在电力调度任何可以写成“先决策—后调整”形式的两阶段问题都可以直接复用这个方法框架。