鲁棒优化入门(4)-两阶段鲁棒优化及行列生成算法(C&CG)超详细讲解(附matlab代码)

本文的主要参考文献：
Zeng B , Zhao L . Solving Two-stage Robust Optimization Problems by A Constraint-and-Column Generation Method[J]. Operations Research Letters, 2013, 41(5):457-461.

1.两阶段鲁棒优化问题的引入

鲁棒优化是应对数据不确定性的一种优化方法，但单阶段鲁棒优化过于保守。为了解决这一问题，引入了两阶段鲁棒优化(Two-stage Robust Optimization)以及更一般的多阶段鲁棒优化，其核心思想是将决策问题分为两个阶段。第一阶段是进行初步决策，第二阶段是根据第一阶段的决策结果制定更好的决策策略，应对数据不确定性的影响。这种方法可以降低保守性，提高鲁棒性。

假设一阶段和二阶段决策问题都是线性规划，并且不确定性集合U是一个有限的离散集合或者多面体集。使用y表示第一阶段决策变量，x表示第二阶段决策变量，表示不确定矢量。在此假设下的两阶段鲁棒优化的一般形式为：

其中：

向量c,b,d,h和矩阵A , G , E , M都是确定性的数值，不确定性体现在向量u上。注意到第二阶段优化的约束条件F(y,u)是关于不确定性u的线性函数。

原文献中提供了以运输问题作为算例，具体如下：

其中，yi为0-1变量，表示是否在i地建设仓库，zi表示仓库i储存的商品数量，xij表示从i仓库到j客户运送的商品数量，fi表示建设仓库i的固定成本，ai表示仓库i存储商品的单位成本，cij表示从i仓库到j客户运送单位商品的成本，ki表示仓库i的最大容量，dj表示客户j的需求。

不确定变量为客户的需求，表达方式如下：

具体算例参数：

根据上面的公式，我们可以写出各个参数矩阵以及变量的表达式：

用matlab代码表示：

%% 参数矩阵f = [400;414;326];a = [18;25;20];k = 800;C = [22, 33, 24;33, 23, 30;20, 25, 27];d_ = [206; 274; 220];d_wave = 40;gamma = [1.8,1.2];P = [1 1;1 1;1 0];%% 决策变量y = binvar(3,1);z = sdpvar(3,1);x = sdpvar(3,3,’full’);d = sdpvar(3,1);g = sdpvar(3,1);

可以尝试求解一下这个确定性优化问题，和后面的两阶段鲁棒优化进行对比：

%% 目标函数objective = f'*y + a'*z + sum(sum(C.*x));%% 约束条件Constraints = [];Constraints = [Constraints , z >= 0 , x >= 0 , g >= 0 , g <= 1];Constraints = [Constraints , z <= k*y];Constraints=[Constraints , sum(x) = d];Constraints=[Constraints ,d == d_ + g*d_wave];Constraints=[Constraints ,g'*P <= gamma];%% 设置求解器ops=sdpsettings('verbose', 3, 'solver', 'gurobi');sol=optimize(Constraints,objective,ops);

优化结果为：

进一步把算例写成两阶段鲁棒优化的形式：

针对这个两阶段鲁棒优化问题，可以分别采用Benders对偶割平面法和C&CG算法进行求解。

2.Benders对偶割平面法

2.1基本原理

Benders对偶割平面法可以用于解决两阶段鲁棒优化问题，首先将两阶段鲁棒优化问题分解为两部分：主问题(Master Problem,MP)和子问题(Subproblem,SP)。主问题包含第一阶段的决策变量y以及仅与y有关的约束和子问题返回的割，还包括辅助变量η，用于评估第二阶段目标函数的取值。子问题包含第二阶段的决策变量x和不确定变量u，旨在给出第二阶段目标函数值的一个界限值。针对式(1)中描述的两阶段鲁棒优化问题，其主问题MP可以写成：

主问题是一个线性规划问题。子问题SP则为：

而子问题是一个双层线性规划问题(如果不知道双层规划的概念，可以去看看我之前的几篇博客双层优化入门-CSDN博客)，其中上层优化的决策变量是u，下层优化的决策变量是x，而且在下层优化中，变量y和u的值都是确定的，可以视为参数。

对于双层优化形式的子问题的求解，主要有以下几种方式：

1.通过对偶变换将双层优化问题转为单层优化问题，再进行求解，可以使用智能优化算法、等价线性化、二次规划求解器(例如gurobi)等方式进行求解；

2.采用智能优化算法进行求解(可参考博客双层优化入门(3)—基于智能优化算法的求解方法)；

3.采用KKT条件进行求解(可参考博客双层优化基本原理与求解方法、基于yalmip的双层优化求解)。

在这里我们使用KKT条件来求解子问题，可以将双层优化的子问题转换为下列单层优化的形式：

（补充说明）

上面给出的式子是经过化简的，并不是最初的KKT条件，因此有读者和我反映，看起来有点懵圈。我在此把详细的转换步骤写一下：

首先由于子问题的内层优化问题是min形式，需要把约束条件都转为≥0的形式：

根据这个优化问题的形式，可以写出子问题内层优化的拉格朗日函数：

其中，π和θ都是拉格朗日乘子。根据拉格朗日函数进一步写出KKT条件：

这就是初始的KKT条件，接下来看一下文献中是如何对这组式子进行化简的。首先根据式(1)和(4)可以得到：

然后，根据式(1)和(3)可以得到：

这样化简之后，就得到了和原文中相同的形式。根据化简的结果，大家应该能知道化简的目的：消去变量θ，减少优化问题中变量的数目。实际上几乎所有的KKT条件都可以这样进行处理，之后有机会再和大家讲讲，这里不再赘述。

化简后的KKT条件也存在非线性形式，可以使用大M法引入二进制中间变量进行线性化，将其转换为混合整数规划的形式：

对于一般形式的两阶段鲁棒优化问题(如式(1))，Benders对偶切平面算法求解的流程如下：

步骤1：设定目标函数上界UB=＋∞，下界LB=-∞，迭代次数k=0。

步骤2：求解主问题MP：

求出最优解(yk+1*，ηk+1*)，并更新LB=max{LB,cT yk+1*+ηk+1*}；

步骤3：求解子问题SP:

求出子问题的最优目标函数值Qk+1*以及最优解(uk+1*，xk+1*)，并更新UB=min{UB,cT yk*+Q(yk*)}。

步骤4：如果UB-LB ≤ ε(事先设定的运行偏差)，则输出优化结果，并退出循环。否则令k=k+1，将约束添加到主问题MP中并返回步骤2。

从上述步骤中可以看到，算法迭代的过程会不断向主问题添加约束条件，也就是割平面，因此被称为Benders对偶割平面法。

2.2 算例分析

采用文献中给出的运输问题作为算例，使用matlab+yalmip工具箱+gurobi求解器进行求解。

为了求解这个两阶段鲁棒优化问题，我们首先需要把这个优化问题分解成主问题和子问题。而且为了方便理解，重写成符合标准两阶段鲁棒优化问题的形式，其中重写后的优化问题部分变量或系数矩阵和原优化问题中重复，我都加了上标一撇(‘)以示区别，具体步骤如下：

主问题MP_BD：

子问题SP_BD：

步骤1：设定目标函数上界UB=＋∞，下界LB=-∞，迭代次数k=0。

步骤2：求解主问题MP_BD，得到最优解(,)，并更新LB=max{LB,}；

步骤3：求解子问题SP_BD，得到子问题的最优目标函数值Qk*以及最优解(uk*，xk*)，并更新UB=min{UB,}。

步骤4：如果UB-LB ≤ ε(事先设定的允许偏差)，则输出优化结果，并退出循环。否则令k=k+1，将约束添加到主问题MP中并返回步骤2。

采用matlab编程进行求解，结果如下：

与确定性优化的结果对比如下：

变量	确定性优化	两阶段鲁棒优化
最优目标函数值	30566	33680
y	[1 1 0]	[1 0 1]
z	[426 274 0]	[255.2 0 516.8]
x
g	[0 0 0]	[0 1 0.8]
d	[206 274 220]	[206 314 252]

3.列与约束生成算法(C&CG)

3.1 基本原理

Benders对偶割平面法通过将两阶段鲁棒优化分解为主问题和子问题，不断交替求解，并将子问题的求解结果作为主问题增加的约束条件，以此达到迭代收敛的目的。C&CG算法和Benders对偶割平面法原理有一些相似，但也存在明显的差异，其基本原理如下：

在C&CG算法的求解过程中，主问题中首先将不考虑不确定变量的影响，作为一个确定性优化进行求解。但不确定变量的最恶劣场景肯定会对主问题的决策产生影响，所以C&CG算法的本质思想就是在确定性优化求解的基础上，不断添加相对恶劣的场景以及对应的子问题决策变量和约束条件，从而使目标函数上界和下界不断得到改进，直到算法收敛，这种思想也被称为“追索权(recourse)”决策：在主问题不考虑不确定性的基础上，不断根据子问题决策带来的不确定性进行修正决策，来保证最终解的鲁棒性。

根据C&CG算法的基本原理，我们可以想到，如果不确定集是一个离散的有限集，那么也可以通过枚举法来求解两阶段鲁棒优化问题。但实际情况肯定不会这么简单，还是需要通过C&CG算法交替迭代求解，具体步骤如下：

主问题MP_CCG：

其中，xl是第l次迭代添加的决策变量，是第l次迭代子问题的最恶劣场景，因为第l次迭代时子问题中已经求解了的值，所以在主问题中就可以看作已知的参数。

子问题SP_CCG：

步骤1：设定目标函数上界UB=＋∞，下界LB=-∞，迭代次数k=0。

步骤2：求解主问题MP_CCG，得到最优解，并更新LB=max{LB,}；

步骤3：求解子问题SP_CCG，得到子问题的最优目标函数值Qk*以及最优解(uk*，xk*)，并更新UB=min{UB,}。

步骤4：如果UB-LB ≤ ε(事先设定的允许偏差)，则输出优化结果，并退出循环。否则转到步骤5.

步骤5：判断子问题是否存在最优解。

若Q(yk+1*)<+∞，即子问题存在最优解，则创建变量xk+1并给主问题MP_CCG添加以下约束：

其中uk+1*是第k次迭代时子问题求解出来的最恶劣场景。

随后更新k=k+1，并转到步骤2。

若Q(yk+1*)=+∞，即子问题不存在最优解，则创建变量xk+1并给主问题MP_CCG添加以下约束：

随后更新k=k+1，并转到步骤2。

从C&CG算法的步骤可以看到，在迭代的过程中在不断地向主问题添加决策变量以及约束条件，也就是优化问题的行和列都在增加，因此才被称为column-and-constraint generation (C&CG)算法。

原文献中解释了C&CG算法和Benders对偶割平面算法的区别，具体如下：

(1) 在主问题中，C&CG算法通过在每次迭代中引入一组新变量来增加解空间的维数，而Benders对偶割平面算法则使用相同的变量集合。

(2) 在处理可行性问题方面，C&CG算法提供了一种通用的方法，而Benders对偶割平面算法是针对特定问题的。

(3) 在计算复杂度方面，与Benders对偶割平面算法相比，C&CG算法在求解主问题时使用更多变量和约束条件。然而，若第二阶段的决策问题为LP问题，根据原文中的命题1和2，Benders对偶割平面算法的复杂度为O(pq)，C&CG算法的复杂度将为O(p)(p是不确定集U的极值点数，q为满足GTπ≤b和π≥0的集合{π}中的极值点数量，具体见原文献中的描述)。因此C&CG算法的收敛速度要更快。

(4) 在求解问题的能力方面，Benders-dual对偶割平面算法需要将第二阶段的问题转换为线性规划问题，而C&CG算法不关心第二阶段的变量类型。如果第二阶段是混合整数规划，可以采用嵌套C&CG算法进行求解。