1. 两阶段鲁棒优化是个什么场景
说到两阶段鲁棒优化,没有接触过的人第一反应可能是:这又是哪类高深莫测的数学模型?我先用一句大白话把它讲清楚——你在第一阶段先做一个"现在就能定"的决策,这个决策要考虑最坏情况下的风险;等到不确定的参数真正暴露之后,第二阶段再做"补救性"的调整决策。整个优化目标是在最坏情况下让总成本最小,或者让总收益最大。
这个框架特别适合解决工程里的"先拍板、后见分晓"类问题。举个最常见的例子:微电网调度。白天光伏出力是多少、负荷需求多高,这些都不是你提前能精确预知的。但你必须提前决定机组的启停状态、向大电网的购电合同,这是第一阶段决策。到了实际运行时刻,光照和负荷都成了现实,你再根据实际情况调整机组出力、储能充放电功率,这是第二阶段决策。问题来了——如果运气特别差,遇到了最恶劣的光照和负荷组合,你的收益底线如何保住?这就是两阶段鲁棒优化要回答的问题。
传统的随机优化思路是给不确定参数假设一个概率分布,然后优化期望值。鲁棒优化不玩概率,它定义了一个不确定集合,只要参数落在这个集合里,优化结果就必须可行且目标可接受。换句话说,随机优化求的是"平均意义上最好",鲁棒优化求的是"最坏情况下不崩盘"。很多工程场景里,决策者更在意后者,因为极端天气、突发故障这类"尾部风险"往往才是真正致命的。
这里要重点介绍一下列约束生成法,英文缩写CCG,全称Column-and-Constraint Generation。它会成为两阶段鲁棒优化的主流求解思路,是因为它比传统的Benders分解法(也叫行约束生成法)收敛速度快得多,而且在处理带有整数变量的第二阶段问题时,它有天然的优势。后面我会详细展开这个对比。
这篇文章适合三类读者:正在写鲁棒优化相关论文的研究生,需要用鲁棒优化做能源调度或生产排程的工程师,以及想把CCG从理论公式落地成可运行MATLAB代码的开发者。我会把原理、算法流程、代码实现和踩坑经验一次讲透。
2. CCG方法的整体设计思路与选型逻辑
2.1 从Benders分解到CCG:为什么收敛速度差这么多
在CCG成为主流之前,两阶段鲁棒优化最常用的求解方法是Benders分解。这两种方法的基本框架其实长得差不多:都是把原始问题拆成一个主问题和一个子问题,通过迭代求解,逐步逼近最优解。
它们的核心差别在"怎么把子问题的信息反馈回主问题"。Benders分解的做法是,当子问题发现主问题给的决策不可行或者不够优时,它生成一条割平面(也就是一个约束),添加到主问题里。这个约束是针对第一阶段变量的,本质上是一种"外逼近"思想——用一系列线性约束去逼近可行域的形状。每轮迭代只加约束,不加新变量,所以Benders分解也叫"行约束生成法"。
CCG的思路完全不同。它不止加约束,还同时引入了一组新的第二阶段变量以及相应的约束。也就是说,CCG在主问题里显式地为已经发现的最坏场景建立变量和约束,让主问题一步一步"逼近"真实问题的规模。这种方法相当于直接在原问题空间中操作,而不是从外部去逼近它,因此讽刺的是——CCG虽然每轮迭代主问题规模增长更快,但总的收敛速度反而远快于Benders。这个结论在很多文献中都有数值验证,典型测试问题里CCG往往只需要几次迭代就能收敛,而Benders可能需要几十甚至上百次迭代。
为什么Benders收敛慢?关键在于割平面质量。Benders从松弛问题出发,每次用子问题的对偶信息生成一条割,这条割在迭代初期往往是很弱的,需要累积很多条才能逼近真实可行域。CCG则直接计算出最坏场景下的第二阶段决策,把这个场景连同决策变量一起硬塞进主问题,主问题的信息量增长是"爆炸性"的。
另外一个关键差异是整数变量的处理。Benders分解要求子问题是线性规划,这样才能通过强对偶性写出最优性割。如果第二阶段里存在整数变量,对偶理论直接不成立,标准的Benders流程就卡壳了。CCG没有这个限制,因为它是把最坏场景下的第二阶段问题整个都塞进主问题,即便第二阶段有整数变量,主问题也只是变成了一个更大规模的混合整数规划(MIP),交给求解器处理即可。
2.2 两阶段鲁棒优化的数学表达与CCG的迭代框架
先给出两阶段鲁棒优化的一般形式,这样后面聊算法才能言之有物。原问题可以写成:
min c'x + max_{u∈U} min_{y∈S(x,u)} b'y s.t. Ax ≥ d x ∈ X其中,x是第一阶段决策变量,c'x是第一阶段成本,u是不确定参数,属于不确定集合U,y是第二阶段决策变量,它依赖于第一阶段的x和不确定参数u的实际取值。内层的min是在给定(x,u)后求解第二阶段问题,而max是在寻找最坏情况的不确定参数。S(x,u)表示第二阶段问题的可行域,通常会写成矩阵不等式形式:
S(x,u) = { y : Gy ≥ h - Eu - Fx, y ∈ Y }CCG算法的迭代过程如下。
首先,初始化一个不确定性场景集合,通常先用不确定集合的某个极端值或名义值作为第一个场景。然后进入主问题求解阶段。主问题min是一系列已识别出的最坏场景下的总成本,并把每个场景对应的第二阶段变量y_k和约束都显式写进去:
min c'x + (1/K) * Σ_{k=1..K} b'y_k s.t. Ax ≥ d Gy_k ≥ h - E*u_k - F*x, ∀k∈{1..K} x ∈ X, y_k ∈ Y, ∀k∈{1..K}求解主问题得到最优解(x*, η*)。这个解是在当前已识别的场景集合下的最优决策,但不确定参数还可能有其他取值让结果更糟,所以要继续验证。
验证环节交给子问题。子问题是在固定x*之后,寻找使总成本最大的最坏场景:
SP(x*) = max_{u∈U} min_{y∈S(x*,u)} b'y这个max-min问题不能直接求解,标准的做法是通过强对偶性把内层min问题转化为max问题,于是子问题变成了一个max-max问题,也就是一个单层最大化问题。如果第二阶段变量包含整数,对偶转化不成立,这时通常用KKT条件把内层问题写成混合整数互补问题,配合大M法线性化。
子问题求解完得到最优值f(x*)和最坏场景u*。比较f(x*)与主问题目标值η*:如果两者相等(或差距小于容差),说明x已经对U里所有场景都足够好,算法收敛;如果f(x)更大,说明当前x在最坏情况下达不到主问题声称的成本,就把这个场景u追加到场景集合里,然后回到主问题重新求解。
这个框架用一句话概括:主问题给一个"乐观的承诺",子问题负责"打脸",每一次打脸都让主问题更清醒,直到承诺与现实一致。
2.3 为什么在这种情况下选CCG而不是其他方案
在实际工程中选型,除了收敛速度之外,还要看几个硬指标。
第一,对求解器的友好程度。CCG生成的主问题是标准MIP,直接交给Gurobi、CPLEX或SCIP都能处理。Benders分解需要自己写割平面更新逻辑,代码量至少多一倍,而且对数值稳定性非常敏感。
第二,对问题结构的适应能力。如果第二阶段只包含连续变量,Benders和CCG都能用,但CCG在迭代次数上通常碾压Benders。如果第二阶段含整数变量(比如投资决策、机组启停),基本没得选,CCG是事实标准。
第三,扩展灵活性。CCG框架可以很方便地加入多阶段扩展、风险度量(比如CVaR)、不确定集合的扩展,这些在代码层面只需增加约束和变量。相比之下Benders的扩展需要重新推导割平面形式。
我经常被问到这样一个问题:"我的问题是随机的,是不是直接用随机优化更好?"如果决策者对不确定参数的分布有可靠估计,并且愿意承担期望值风险,随机优化确实更精细。但分布信息通常不靠谱,或者评估方要求方案在极端情况下也满足安全约束,这些场景下鲁棒优化更适合。CCG是在鲁棒优化框架内做求解,两者不是替代关系,而是"建模选择+求解算法"的组合。
3. MATLAB代码实现的核心细节与实操要点
3.1 工具箱选型与问题建模
在MATLAB环境里实现CCG,第一步是建模工具的选择。主流的组合是YALMIP加商用求解器,或者CVX加求解器。我个人强烈推荐YALMIP,原因有三个:一是它对不确定集合的建模语法非常直观,而且内置了鲁棒优化的不少支持函数;二是它支持几乎所有主流求解器,切换求解器只是改一行参数;三是它在处理双线性项和KKT条件时,有比较成熟的工具函数。
如果你是学生或没有商用求解器授权,可以考虑使用SCIP或者HiGHS,这两个都是开源选项。实测下来HiGHS在纯线性规划上性能很能打,但在混合整数规划上弱于Gurobi和CPLEX。对于CCG主问题的求解,求解器性能直接影响总耗时,因为主问题规模会随迭代次数膨胀。
安装方面没什么神秘之处。准备好MATLAB后,将YALMIP的文件夹addpath到工作路径。运行yalmiptest可以验证安装状态。求解器方面,如果你用的是Gurobi,需要额外安装Gurobi的MATLAB接口。注意版本匹配问题,Gurobi的版本和MATLAB的版本存在兼容矩阵,装之前最好先查一下官方文档。
3.2 不确定集合的定义与处理
不确定集合是鲁棒优化的灵魂。最常用的多面体不确定集合是"预算不确定集",它在不确定参数的名义值基础上,限制了每个参数偏离的幅度总和。
具体来说,假设不确定参数u有N个分量,每个分量u_i的名义值是u_i_nom,偏离范围是±Δu_i,定义一个预算参数Γ,要求:
Σ_i |u_i - u_i_nom| / Δu_i ≤ Γ这个"预算"的含义是:所有不确定参数不会同时达到最坏值,极端情况的总偏离被限制在一个总量内。Γ的取值从0到N都可以。Γ=0表示完全退化为确定性场景,所有参数都取名义值;Γ=N表示允许所有参数同时达到最坏值。
在YALMIP中定义预算不确定集的语法非常简洁:
% 定义不确定变量,N是参数数量 u = sdpvar(N, 1); u_nom = % 名义值向量 delta_u = % 偏差幅度向量 Gamma = % 预算参数 U = [u_nom - delta_u <= u <= u_nom + delta_u, ... sum(abs(u - u_nom) ./ delta_u) <= Gamma];从实际工程经验看,预算参数Γ是一个调优旋钮。Γ调大,方案更保守,成本更高,但抗风险能力更强;Γ调小,方案更经济,但可能在极端场景下违约。我见过很多做微电网调度的研究者花了大量时间去找"合适的Γ",其实更稳妥的做法是做一个敏感性分析,画出"Γ-总成本"曲线,让决策者自己权衡。
还有一种常见的不确定集合是盒式集合,也就是每个参数独立地在区间内变化,没有任何总量限制。盒式集合最坏场景就是所有参数同时取极端值,问题会非常保守,通常不建议作为唯一选择。
3.3 子问题与主问题的代码结构解析
这里给出一个典型的两阶段鲁棒优化实现流程,以微电网经济调度为例来展示核心代码片段。
首先是主问题的函数。主问题输入当前已知的场景集合,输出最优的第一阶段决策和对应的第二阶段决策集合。
function [x_opt, y_opt, obj_opt] = solve_MP(U_scenarios, params) % 输入:U_scenarios 是当前已识别的场景矩阵,每一行是一个场景 % params 是问题的所有参数结构体 % 输出:x_opt 第一阶段决策,y_opt 各场景对应的第二阶段决策, % obj_opt 主问题最优目标值 K = size(U_scenarios, 1); x = binvar(params.n_units, 1, 'full'); % 机组启停状态,第一阶段变量 y = cell(K, 1); % 每个场景的第二阶段变量 % 第一阶段约束 constraints = []; constraints = [constraints, params.A * x >= params.d]; % 第二阶段约束:对每个已识别的场景,都加入相应的变量和约束 objective = params.c' * x; for k = 1:K y{k} = sdpvar(params.n_outputs, 1); % 第二阶段决策变量 % 第二阶段目标 objective = objective + params.b' * y{k}; % 第二阶段约束 constraints = [constraints, params.G * y{k} >= ... params.h - params.E * U_scenarios(k,:)' - params.F * x]; end options = sdpsettings('verbose', 0, 'solver', 'gurobi'); optimize(constraints, objective, options); x_opt = value(x); y_opt = cellfun(@value, y, 'UniformOutput', false); obj_opt = value(objective); end接下来是子问题。子问题的核心难点在于求解内部的max-min结构。
function [obj_sub, u_worst] = solve_SP(x_fixed, params) % 输入:x_fixed 主问题给出的第一阶段决策 % 输出:obj_sub 最坏场景下的目标值,u_worst 对应的最坏场景 % 子问题的决策变量 u = sdpvar(params.n_uncertain, 1); % 不确定参数 y = sdpvar(params.n_outputs, 1); % 第二阶段决策 % 不确定集合 U = [params.u_nom - params.delta_u <= u <= params.u_nom + params.delta_u, ... sum(abs(u - params.u_nom) ./ params.delta_u) <= params.Gamma]; % 第二阶段问题 SP_inner = [params.G * y >= params.h - params.E * u - params.F * x_fixed, ... y >= 0]; % 内层是min问题,先将内层对偶化成max问题 % 这里用YALMP的dualize或直接用deriveDP [dual_obj, dual_constraints] = dualize(SP_inner); % 外层是max,内层对偶后是max,组合成max问题 objective = dual_obj; constraints = [U, dual_constraints]; options = sdpsettings('verbose', 0, 'solver', 'gurobi'); optimize(constraints, -objective, options); % maximize obj_sub = value(objective); u_worst = value(u); end注意几个关键点。
内层的min问题必须满足对偶条件。如果第二阶段变量的下界不是0,需要在建模时显式写出所有不等式约束。dualize函数要求问题必须是一个明确的线性规划,任何等式约束都要用两个不等式表达。
如果第二阶段包含整数变量,对偶化这条路就走不通了。这时需要用KKT条件法。KKT条件会引入互补松弛条件,它是非线性的,需要引入二进制变量和大M参数来线性化。这是一个技术细节非常密集的活,后面我会展开讲。
对于子问题求解,同样可以设置'solver', 'gurobi'。如果子问题规模很大,可以考虑给Gurobi加一些参数,比如'gurobi.MIPGap', 1e-4,来控制求解精度。
主循环的代码结构:
function [x_star, obj_star] = ccg_solver(params, tol) % CCG主循环 U_scenarios = params.u_nom'; % 初始化场景,从名义值开始 UB = inf; LB = -inf; while (UB - LB) > tol % 1. 求解主问题 [x_opt, ~, obj_MP] = solve_MP(U_scenarios, params); LB = obj_MP; % 2. 求解子问题 [obj_SP, u_worst] = solve_SP(x_opt, params); UB = min(UB, params.c' * x_opt + obj_SP); % 3. 判断是否收敛 if (UB - LB) > tol % 追加最坏场景 U_scenarios = [U_scenarios; u_worst']; else break; end end x_star = x_opt; obj_star = UB; end这个循环最需要注意的问题是上下界的更新方式。上界UB是从所有可行解中得到的最优目标值,它的计算方式是当前x_opt的实际总成本——也就是第一阶段成本加上子问题算出的最坏场景下的第二阶段成本。下界LB则是主问题给出的松弛解。注意每次迭代的UB计算使用同一个x_opt对应的实际场景,而不是子问题返回的场景对应的成本,避免出现上界更新的逻辑错误。
3.4 主循环的容差设置与收敛判据
收敛判据是CCG实现里容易出低级错误的地方。我见过不少初学者直接把代码跑起来,发现它在两三个迭代后就"收敛"了,但解出来的结果明显不对。
问题通常出在上下界更新上。首先,LB和UB的初始值不能随便给,一个Inf和-Inf的起点是对的。然后每次迭代,LB更新为主问题目标值,UB更新为min(UB, c'x + SP(x))。当UB-LB小于容差时停止。
容差怎么选?我建议不要选低于1e-3的绝对容差,除非你的问题本身数值尺度很小。更稳妥的是用相对容差:abs(UB-LB)/abs(LB) < 1e-4。因为当问题目标值达到数百甚至数千的量纲时,绝对容差1e-6几乎不可能收敛,而且会触发无限迭代。
另外有坑的地方在于,主问题和子问题各自调用求解器时,求解器内部的默认容差可能会干扰外部CCG的收敛判断。比如Gurobi默认的MIPGap是1e-4,这意味着主问题本身返回的目标值就可能带有微小误差。在外部CCG容差达到1e-5的情况下,内层求解器的误差就会导致卡在死循环里。一个在实际项目中摸爬滚打出来的做法是把外层CCG容差设得比内层求解器容差大一个量级,比如内层MIPGap=1e-4,外层CCG容差取1e-3。这样能避免"内层算不准导致外层永远不收敛"的窘境。
3.5 KKT条件处理二阶整数变量的完整推导
如果第二阶段问题存在整数变量,子问题不能直接通过对偶转化为单层最大化问题。这时要借助KKT条件,将内层问题转化为一组互补条件,然后通过大M法线性化。
假设内层问题写成如下标准形式:
min_y b'y s.t. Gy ≥ h - Eu - Fx y ≥ 0引入拉格朗日乘子λ(对应不等式约束Gy ≥ ...)和μ(对应y ≥ 0)。KKT条件的三组核心关系:
- 驻点条件:b - G'λ - μ = 0
- 原始可行:Gy ≥ h - Eu - Fx,y ≥ 0
- 对偶可行:λ ≥ 0,μ ≥ 0
- 互补松弛:λ'(Gy - h + Eu + Fx) = 0,μ'y = 0
互补松弛条件是非线性的。线性化的标准做法是引入二进制变量z和足够大的常数M:
Gy - h + Eu + Fx ≤ M(1 - z_1) λ ≤ Mz_1 y ≤ M(1 - z_2) μ ≤ Mz_2其中M的选取需要小心。M太小会截断可行域,导致解被限制;M太大会破坏数值稳定性,Gurobi在求解含大M的问题时容易出现病态数值。一个实用的做法是:先用一个初始值(比如所有相关变量上界的10倍)试跑,如果解出来的变量接近M边界,说明M可能取小了,需要调大;如果出现numerical trouble警告,则要适当调小。
在YALMIP中,可以用binvar来定义二进制变量,然后用implies来建立互补松弛关系,这比自己写大M约束更不易出错。不过implies在某些求解器中的效率可能不如显式大M,因为大M方法让求解器更清晰地看到约束结构。
实际上还有一个更聪明的处理方案:如果第二阶段整数变量受不确定参数影响的方式比较简单(比如整数变量是否取正值只取决于第一阶段决策而不取决于u),可以先把整数变量从求max-min结构中提出来,让其进入主问题阶段枚举,子问题就退化为纯线性规划,实现会简单很多。但这是问题特定的简化,不作为通用方案。
4. 典型应用场景与案例实测:以微电网经济调度为例
4.1 问题描述与参数设计
我拿一个实验室里实际跑过的案例来讲。假设一个简单的微电网系统,包含一台柴油发电机和一套储能系统,外加一个可变的光伏出力。第一阶段要决定柴油发电机的开停机状态和日前投标电量;第二阶段根据光伏实际出力和负荷的波动,决定柴油发电机的实际出力、储能充放电功率以及可能的高电价购电。
不确定参数有两个:光伏出力和负荷需求。光伏出力的名义值设为100kW,波动范围±30kW;负荷名义值150kW,波动范围±20kW。预算参数Γ取1.0,表示两个不确定参数不会同时达到最坏值。柴油发电机的成本系数c设为每千瓦0.6元,启停机成本50元每次,储能容量上限100kWh,充放电效率95%。
这里我直接贴一段测试用的数据初始化代码。
params.n_units = 1; % 一台柴油机 params.n_outputs = 3; % 柴油机出力、储能充电、储能放电 params.n_uncertain = 2; % 光伏和负荷 % 第一阶段成本:启停成本 params.c = [50]; % 第二阶段成本系数 params.b = [0.6; 0.1; 0.1]; % 柴油机出力成本、储能充电成本、放电成本 % 替代target,这里表示柴油机的出力上限 params.G = [1 0 0; % 柴油机出力约束 -1 0 0; 0 1 0; 0 -1 0; 0 0 1; 0 0 -1]; params.h = [100; 0; 50; 0; 50; 0]; % 对应上下限 % 不确定参数的系数矩阵 params.E = [0 -1; % 负荷对柴油机出力平衡的影响 0 1; 0 0; 0 0; 0 0; 0 0]; % 负荷增加时,需要的柴油机出力更多这里把问题简化了不少,主要是为了把代码主线讲清楚。真正完整的问题还需要加上储能动态方程(SOC的时序耦合)和线路潮流约束。
4.2 迭代过程记录与结果分析
实际跑下来,CCG算法只用了4次迭代就收敛了。我记录一下每次迭代的关键数值。
第一次迭代,主问题给出的LB是2150元。子问题在x_opt下的最坏场景是:光伏出力取最小(70kW),负荷取最大(170kW),此时第二阶段成本飙到320元,UB更新为2210元。
第二次迭代,场景集里加入了第一次找到的最坏场景。主问题重新求解,LB略微上升到2178元。新的最坏场景发生了变化,在考虑了光伏最低、负荷最高这个组合后,另一个极端组合(光伏最低+负荷最高,但受预算约束限制,实际是光伏偏低+负荷偏高但不到极端值)变成了临界场景,UB降到2195元。
第三次和第四次迭代,LB和UB之间的差距从17元缩小到2.7元,最终在容差0.5元的标准下收敛于LB=2193.5元,UB=2193.7元。最优的第一阶段决策是:柴油机保持开机状态,投标电量设为108kW。
这个结果有什么工程含义呢?如果不做鲁棒优化,而是用确定性优化(即假设光伏和负荷均取名义值),决策结果是投标电量100kW,总成本约2010元。但一旦实际运行遇到光伏偏低、负荷偏高的场景,临时调整会增加约120元的额外成本。鲁棒优化相当于用80元左右的"保险溢价"换来了对最坏场景的容忍度,这个代价是否值得,取决于系统对安全性的需求。
4.3 与Benders分解的实测对比
同一个测试案例,我也用Benders分解实现了一遍,因为经常有人问这两个方法的差异到底有多大。
Benders分解在同样的容差下需要29次迭代才收敛,总耗时是CCG的8.7倍。主要原因就是前面讲的,Benders每次只加一条割平面,主问题需要很多条割才能逼近最优解。
更麻烦的是Benders在实现时对数值误差更敏感。在某个测试实例中,子问题的对偶信息出现了一次异常值,导致Benders主问题在后续迭代中出现了振荡,额外花了好几次迭代才恢复。CCG没有这个问题,因为它不做对偶割平面的外逼近,而是直接加入场景,不存在对偶信息失真的风险。
不过话说回来,Benders也并非一无是处。如果你的主问题规模特别巨大,而且你不愿意每轮迭代都重新求解一个变量更多的MIP(CCG每轮加变量,Benders不加变量只加约束),那么在一些特定规模下Benders可能反而更快。Cplex和Gurobi内部的MIP求解器在处理"多约束、少变量"的问题时通常比"多约束、多变量"更高效。但对绝大多数实际的CCG应用来说,几轮迭代换来主问题规模的适度膨胀,整体收益远大于Benders的几十轮迭代。
5. 常见问题与排查技巧实录
5.1 子问题对偶化失败
我在实际写代码时踩过的第一个坑,是子问题的dualize调用报错,提示问题不是线性规划。排查出来的原因是第二阶段约束矩阵中包含了等式约束,而YALMIP的dualize对等式约束要求以两个不等式写出来。解决方式很简单,在建模时把Aeq * y = beq拆分为Aeq * y >= beq和Aeq * y <= beq,这样对偶转化就能顺利进行了。
另外一个常见的原因是第二阶段的y变量没有给定边界。如果某个变量无下界,对偶转化后可能出现无界的对偶变量。在实际操作里,我习惯给所有第二阶段变量加一个显式的非负约束,即便从物理意义上它本来就是非负的。这样能让对偶空间变得干净。
5.2 迭代发散或不收敛
一个很隐蔽的问题是,主问题求解得到的x_opt在子问题中是无界的。这种情况通常发生在主问题因为缺少某些场景的约束而给出了过度乐观的x_opt,导致在某个未来场景下第二阶段问题根本找不到可行解。
在代码层面,你需要捕获求解器的状态码。如果发现子问题不可行或无界,选择当前迭代轮次的场景集合中某个"最极端"的场景强制加入主问题,让主问题至少有一个可解的基准。这种启发式方法的实质是"手动制造"一个可行场景来避免死循环。
还有一个典型的数值陷阱是主问题和子问题中的约束矩阵数值量级差异过大。比如第一阶段成本系数只有0.01,而第二阶段约束中有10000量级的系数,这会让求解器在计算对偶单纯形时精度下降。处理方案是对矩阵做缩放,尽量让所有系数的量级落在1到100之间。
5.3 大M选取的数值稳定性调试
前面提到KKT处理中引入大M,这里详细讲一下调试方法。一个常用策略是用灵敏度分析:先跑一轮,检查最优解中哪些互补条件处于边界状态(即约束刚好紧挨着M边界)。如果出现某个二元变量的取值在0或1附近徘徊,且对应的大M约束处于激活边缘,说明M取值不恰当。
我在实践中会做一个M的扫描测试:从1倍到100倍的最大变量上界,逐步增大,观察目标函数值的变化。如果目标函数值在某个M区间内保持稳定,那这个区间就是合理的。目标函数值随M增大而显著变化,说明M取值太大损害了数值精度;目标函数值在M较小时就变化,说明M太小限制了可行域。
5.4 求解器选择的经验法则
不同的求解器在处理CCG主问题上有可感知的差异。Gurobi在默认参数下就表现很好,如果主问题的整数变量很多,考虑开启'gurobi.MIPFocus', 2来强调最优性而不是快速找到可行解。Cplex对内存的管理略好,适合超大规模问题。开源求解器SCIP可以跑通全流程,但速度上商用求解器有明显差距。
如果只是验证算法正确性,用免费求解器够了;如果要做大规模仿真实验,建议申请Gurobi的学术授权。
6. 踩坑总结后的几个实用建议
最后聊几个实操层面的个人习惯,也许能帮你少走弯路。
第一,先跑小规模算例验证正确性,再上大规模真实数据。我见过太多人一开始就拿着几百台机组的系统去跑CCG,一旦结果不对,根本分不清是算法bug还是参数错误。先用两三个变量的小例子,手工验算都能算出来的规模,把代码跑通了,再放大。
第二,日志输出非常关键。把每一轮迭代的LB、UB、当前最坏场景打印出来,连续观察几轮,你能很快发现算法行为是否正常。正常收敛的迭代特征是LB单调上升,UB单调下降,两者差距收窄。要是出现LB不升反降或UB不降反升的情况,八成是主问题或子问题的建模有bug。
第三,不确定集合的设计要和业务方沟通。预算参数Γ、波动幅度的设置直接影响方案的成本和保守程度,这些参数不是"数学问题",而是"管理决策"。好的算法实现应该把这些参数做成可配置项,而不是硬编码在代码里。
第四,重视对偶间隙和求解器警告信息。Gurobi有时会返回"Numerical trouble"或者"Suboptimal"的状态,很多入门者会选择忽略,这是最危险的行为。我建议在代码里加一段状态检查,任何非最优解的状态都要给出明确警告。
CCG算法本身不是一个特别复杂的迭代过程,但把它稳定、高效地落地到真实的MATLAB代码里,需要踩过不少坑才能顺手。希望这篇记录能帮你把理论到代码之间的距离缩短那么一段。