Two-stage RO: part 1

作者：刘兴禄，清华大学博士在读

2021.09.28

特别感谢Zeng Bo老师非常耐心地解答我的疑惑，您的解答对我来讲非常重要。

本文非常详细、完整地介绍了求解Two-stage Robust Optimization Model的一种精确算法：Column and Constraint Generation Algorithm。该算法由Zeng Bo, Zhao Long于2013年首次提出，相关论文发表在期刊《Operations Research Letters》上，论文题目为《Solving two-stage robust optimization problems using a column-and-constraint generation method》。本文对该论文进行了非常细致的解读，包括：

• 对鲁棒优化问题进行了分类；
• 对两阶段鲁棒优化模型进行了详细解析；
• 对Bender-dual算法进行了讲解；
• 对Column-and-constraint generation method进行了非常详细的解读、完整的推导和完整的解释；
• 对论文中的案例进行了相当详细的推导；
• 对Column-and-constraint generation method进行了完整的代码复现（Python调用Gurobi实现）。

两阶段鲁棒优化问题（Two-stage Robust Optimization）

鲁棒优化是一类考虑参数不确定的数学规划问题，是运筹学中比较高阶的方法。近几年，关于鲁棒优化的研究越来越火热。常见的鲁棒优化问题包括基本的鲁棒优化、多阶段鲁棒优化、分布式鲁棒优化等。

鲁棒优化旨在处理优化问题中参数不确定的挑战，致力于优化最坏情况(worst case)下的解，或者最坏参数分布下的解(对应分布式鲁棒优化)，从而使最终得到的解具有非常强的抵抗不确定性或者风险的能力。

处理参数不确定的优化问题，还有一类方法，叫作随机规划，关于二者的区别，见往期推文。

两阶段鲁棒优化问题，是多阶段鲁棒优化的一个特例(multi-stage robust optimization)。在两阶段鲁棒优化问题中，一般包含两个不同层级的决策变量，又称为第一阶段的决策(first-stage decision)和第二阶段的决策(second-stage decision)。比如：一个联合优化问题，同时考虑

• 战略层决策和战术层决策；
• 战术层决策和操作层决策
• …

一个典型例子，就是选址-配送规划问题。选址决策属于战略层决策，配送问题，可以视为战术层或者作业层决策。

两阶段鲁棒优化问题描述如下：

• 第一阶段的决策相关的参数是确定的，并且第一阶段的决策需要首先做出；
• 第二阶段的决策的相关参数有一些是不确定的。第二阶段的决策需要在第一阶段的决策确定之后，等到第二阶段参数的揭示后，再做出相应的第二阶段决策。
• 两阶段鲁棒优化的目标为：对两阶段决策进行联合优化，同时考虑第二阶段参数的不确定性。也就是优化在第二阶段参数取到最坏情况下的两个阶段决策对应的总的目标值。

下面我们用标准的数学语言来描述两阶段鲁棒优化问题。两阶段鲁棒优化问题的数学模型如下所示(参考自文献[1])

min ⁡ y c T y + max ⁡ u ∈ U min ⁡ x ∈ F ( y , u ) b T x ( 1 ) s . t . A y ⩾ d , ( 2 ) G x ⩾ h − E y − M u , u ∈ U ( 3 ) S y ⊆ R + n ( 4 ) S x ⊆ R + m ( 5 )

minys.t.cTy+maxu∈Uminx∈F(y,u)bTxAy⩾d,Gx⩾h−Ey−Mu,u∈USy⊆Rn+Sx⊆Rm+(1)(2)(3)(4)(5)$$\begin{array}{rlrl}\underset{\mathbf{y}}{min}& {\mathbf{c}}^{\mathbf{T}}\mathbf{y}+\underset{u\in \mathcal{U}}{max}\underset{\mathbf{x}\in \mathbf{F}(\mathbf{y},u)}{min}{\mathbf{b}}^T\mathbf{x}& & (1)\\ s.t.& \mathbf{A}\mathbf{y}⩾\mathbf{d},& & (2)\\ & \mathbf{G}\mathbf{x}⩾\mathbf{h}\mathbf{-}\mathbf{E}\mathbf{y}\mathbf{-}\mathbf{M}u,u\in \mathcal{U}& & (3)\\ & {\mathbf{S}}_{\mathbf{y}}\subseteq {\mathbb{R}}_{+}^n& & (4)\\ & {\mathbf{S}}_{\mathbf{x}}\subseteq {\mathbb{R}}_{+}^m& & (5)\end{array}$$

ymin​s.t.​cTy+u∈Umax​x∈F(y,u)min​bTxAy⩾d,Gx⩾h−Ey−Mu,u∈USy​⊆R+n​Sx​⊆R+m​​​(1)(2)(3)(4)(5)​

其中:

• $\mathbf{y}$为第一阶段的决策变量;
• $\mathbf{x}$为第二阶段的决策变量;
• $\mathbf{c}$为第一阶段的决策变量相关的参数，是确定的;
• $u$为第二阶段的决策变量相关的参数，是不确定的，且其不确定参数的取值范围由不确定集$\mathcal{U}$刻画；具体的例子见往期推文；
• 目标函数中，第一项${\min }_{\mathbf{y}}{\mathbf{c}}^{\mathbf{T}}\mathbf{y}$为第一阶段决策对应的目标函数。该部分决策需要首先确定；
• 目标函数中，第二项${\max }_{u\in \mathcal{U}}{\min }_{\mathbf{x}\in \mathbf{F}(\mathbf{y},u)}{\mathbf{b}}^T\mathbf{x}$为第二阶段决策对应的目标函数。该部分决策需要在第一阶段决策确定之后再做出；该部分目标函数旨在找到worst case;也就是在最坏情况下，第二阶段的目标值；
• 整个目标函数，就是在优化在最坏情况下的总目标函数，使得最坏的情况最好，使得解具有非常好的鲁棒性。
• 约束(2)只跟第一阶段的决策有关；
• 约束(3)中跟第一阶段，第二阶段的决策都有关，且包含不确定参数$u$。

以上就是两阶段鲁棒优化的一个简要介绍。目前还比较晦涩的数学符号，我们将在下文进行更通俗的解读。

• 注意，如果$u$的取值可能性非常少，我们就可以通过穷举$u$的取值，显式地将每个可能的$u$对应的约束都加进模型，然后将两阶段鲁棒优化模型等价为一个可直接求解的一阶段数学规划模型，从而可以达到直接求解上述Two-stage robust optimization模型的目的。因为最坏情况，一定是对应着某一种$u$的取值，而我们已经穷举了所有$u$的取值可能，且将其显式地写进了模型，因此原两阶段鲁棒优化模型可以直接被解决。但是通常，$u$的取值非常多（组合数级别或者无穷多），穷举$u$往往是低效甚至不可行的，因此，需要利用一些其他方法来求解Two-stage robust optimization model。本文将要介绍的Column and Constraint Generation算法就是其中一种（此外，文献[1]中提到的Benders-dual algorithm也是一种）。

两阶段鲁棒优化问题的详细解读

参考文献[1]中提到：

Two-stage RO 的计算非常困难，即使是一个简单的Two-stage RO问题也可能是NP-Hard的。为了克服计算负担，目前主要流行两种求解策略。第一种是使用 approximation algorithms，它假设第二阶段的决策是不确定性的简单函数，如仿射函数。第二类算法是利用Benders分解法求出精确解，即利用第二阶段决策问题的对偶解逐步构造第一阶段决策的价值函数。因此，我们称它们为Benders-dual cutting plane algorithms算法。

文献[1]中关于两阶段鲁棒优化问题模型的描述如下：

该论文探讨的问题，是基于如下的假设：

• 决策变量类型假设：first-stage and second-stage are linear optimization
• 参数不确定性(uncertainty): is either a finite discrete set or a polyhedron. 换句话说，改论文考虑的不确定集为是一个有限的离散集合或者一个多面体集合，比如
  • a finite discrete set：不确定参数的可能取值为一个有限的可选集，例如： u ∈ U = { scenario 1 , scenario 2 , scenario 3 } u \in \mathcal{U} = \{ \text{scenario}_1, \text{scenario}_2, \text{scenario}_3\} u∈U={scenario1​,scenario2​,scenario3​}.
  • a polyhedron: 不确定参数的可能取值的范围为一个多面体，例如 u ∈ U = { 1 ⩽ u 1 ⩽ 2 , 3 ⩽ u 2 ⩽ 5 } u \in \mathcal{U} = \{ 1 \leqslant u_1 \leqslant 2, 3 \leqslant u_2 \leqslant 5\} u∈U={1⩽u1​⩽2,3⩽u2​⩽5}，该不确定集为一个矩形，学术上一般称其为盒子不确定集(box uncertainty set)。

根据上图，两阶段鲁棒优化中，第一阶段的决策需要首先做出，其优化方向为$\min$。

但是第二阶段的决策，却是一个双层问题(bi-level)。具体来讲，目标函数第二项的符号为
$$\underset{u\in \mathcal{U}}{\max }\underset{\mathbf{x}\in \mathbf{F}(\mathbf{y},u)}{\min }{\mathbf{b}}^T\mathbf{x}$$
这个符号需要详细解释一下：

• $\mathbf{x}\in \mathbf{F}(\mathbf{y},u)$的意思: 表示given 不确定参数$u$的一个取值，以及given first-stage决策变量$\mathbf{y}$的取值，我们决策$\mathbf{x}$，使得${\mathbf{b}}^T\mathbf{x}$最小。并且满足如下约束
  F ( y , u ) = { x ∈ S x : G x ⩾ h − E y − M u }
  F(y,u)={x∈Sx:Gx⩾h−Ey−Mu}$$\begin{array}{r}\mathbf{F}(\mathbf{y},u)=\{\mathbf{x}\in {\mathbf{S}}_{\mathbf{x}}:\mathbf{G}\mathbf{x}⩾\mathbf{h}\mathbf{-}\mathbf{E}\mathbf{y}\mathbf{-}\mathbf{M}u\}\end{array}$$
  F(y,u)={x∈Sx​:Gx⩾h−Ey−Mu}​
  with ${\mathbf{S}}_{\mathbf{y}}\subseteq {\mathbb{R}}_{+}^n$ and ${\mathbf{S}}_{\mathbf{x}}\subseteq {\mathbb{R}}_{+}^m$。

其实${\mathbf{S}}_{\mathbf{y}}\subseteq {\mathbb{R}}_{+}^n$ and ${\mathbf{S}}_{\mathbf{x}}\subseteq {\mathbb{R}}_{+}^m$就是说，$\mathbf{x}$和$\mathbf{y}$分别是$n\times 1$维和$m\times 1$维的非负连续决策变量而已。

我们将上图中的模型展开，写成一个比较容易看明白的形式：

min ⁡ y c T y + max ⁡ u ∈ U min ⁡ x ∈ F ( y , u ) b T x ( 1 ) s . t . A y ⩾ d , ( 2 ) G x ⩾ h − E y − M u , u ∈ U ( 3 ) S y ⊆ R + n ( 4 ) S x ⊆ R + m ( 5 )

minys.t.cTy+maxu∈Uminx∈F(y,u)bTxAy⩾d,Gx⩾h−Ey−Mu,u∈USy⊆Rn+Sx⊆Rm+(1)(2)(3)(4)(5)$$\begin{array}{rlrl}\underset{\mathbf{y}}{min}& {\mathbf{c}}^{\mathbf{T}}\mathbf{y}+\underset{u\in \mathcal{U}}{max}\underset{\mathbf{x}\in \mathbf{F}(\mathbf{y},u)}{min}{\mathbf{b}}^T\mathbf{x}& & (1)\\ s.t.& \mathbf{A}\mathbf{y}⩾\mathbf{d},& & (2)\\ & \mathbf{G}\mathbf{x}⩾\mathbf{h}\mathbf{-}\mathbf{E}\mathbf{y}\mathbf{-}\mathbf{M}u,u\in \mathcal{U}& & (3)\\ & {\mathbf{S}}_{\mathbf{y}}\subseteq {\mathbb{R}}_{+}^n& & (4)\\ & {\mathbf{S}}_{\mathbf{x}}\subseteq {\mathbb{R}}_{+}^m& & (5)\end{array}$$

ymin​s.t.​cTy+u∈Umax​x∈F(y,u)min​bTxAy⩾d,Gx⩾h−Ey−Mu,u∈USy​⊆R+n​Sx​⊆R+m​​​(1)(2)(3)(4)(5)​

只不过，论文截图中$\mathbf{y}\in {\mathbf{S}}_{\mathbf{y}}$，就等价于$\mathbf{y}$是下面模型的可行解

min ⁡ y c T y + max ⁡ u ∈ U min ⁡ x ∈ F ( y , u ) b T x s . t . A y ⩾ d , S y ⊆ R + n

minys.t.cTy+maxu∈Uminx∈F(y,u)bTxAy⩾d,Sy⊆Rn+$$\begin{array}{rl}\underset{\mathbf{y}}{min}& {\mathbf{c}}^{\mathbf{T}}\mathbf{y}+\underset{u\in \mathcal{U}}{max}\underset{\mathbf{x}\in \mathbf{F}(\mathbf{y},u)}{min}{\mathbf{b}}^T\mathbf{x}\\ s.t.& \mathbf{A}\mathbf{y}⩾\mathbf{d},\\ & {\mathbf{S}}_{\mathbf{y}}\subseteq {\mathbb{R}}_{+}^n\end{array}$$

ymin​s.t.​cTy+u∈Umax​x∈F(y,u)min​bTxAy⩾d,Sy​⊆R+n​​
其中$u,\mathbf{x}$是一个给定的值。

上述Two-stage Robust Optimization模型，一般是无法直接调用优化求解器进行求解的，需要设计对应的算法进行求解。下面就来介绍求解该模型的算法。

Two-stage RO和Benders-dual cutting plane method

类似于使用标准的Benders Decomposition算法求解确定性问题，Two-stage RO也是可以使用Benders Decomposition算法来求解的。文献中将该算法称之为Benders-dual cutting plane method。

本文仅对Benders-dual cutting plane method做简要介绍。小编其实早就已经写好了Benders-dual cutting plane method的完整笔记和代码，只不过想留着后面再分享，尽请期待！

Benders-dual cutting plane method求解Two-stage RO问题可以类比Benders Decomposition求解确定性问题的情形。我们可以将Two-stage RO分解为：

• 一个主问题(Master Problem)，主要包含第一阶段的决策$\mathbf{y}$及只跟$\mathbf{y}$有关的约束、子问题返回的割，评估第二阶段目标函数取值的辅助变量$\eta$（类似于Benders Decomposition算法中的$q$）;
• 一个子问题(Subproblem)，主要包含第二阶段的决策$\mathbf{x}$和不确定变量$u$，旨在给出第二阶段目标函数值的一个界限（类似于标准的Benders Decomposition算法）。

Benders-dual cutting plane method中的初始主问题和Benders Decomposition算法类似，为

min ⁡ y c T y + η ( 1 ) s . t . A y ⩾ d , ( 2 ) S y ⊆ R + n , ( 3 ) η ⩾ 0. ( 4 )

minys.t.cTy+ηAy⩾d,Sy⊆Rn+,η⩾0.(1)(2)(3)(4)$$\begin{array}{rlrl}\underset{\mathbf{y}}{min}& {\mathbf{c}}^{\mathbf{T}}\mathbf{y}+\eta & & (1)\\ s.t.& \mathbf{A}\mathbf{y}⩾\mathbf{d},& & (2)\\ & {\mathbf{S}}_{\mathbf{y}}\subseteq {\mathbb{R}}_{+}^n,& & (3)\\ & \eta ⩾0.& & (4)\end{array}$$

ymin​s.t.​cTy+ηAy⩾d,Sy​⊆R+n​,η⩾0.​​(1)(2)(3)(4)​
初始主问题并没有任何Benders cut的加入。

下面我们来大致推导Benders-dual cutting plane method中的割平面。首先来看second-stage的部分，也就是

max ⁡ u ∈ U min ⁡ x ∈ F ( y , u ) b T x s . t . G x ⩾ h − E y − M u , u ∈ U S x ⊆ R + m

maxu∈Uminx∈F(y,u)s.t.bTxGx⩾h−Ey−Mu,u∈USx⊆Rm+$$\begin{array}{rl}\underset{u\in \mathcal{U}}{max}\underset{\mathbf{x}\in \mathbf{F}(\mathbf{y},u)}{min}& {\mathbf{b}}^T\mathbf{x}\\ s.t.& \mathbf{G}\mathbf{x}⩾\mathbf{h}\mathbf{-}\mathbf{E}\mathbf{y}\mathbf{-}\mathbf{M}u,u\in \mathcal{U}\\ & {\mathbf{S}}_{\mathbf{x}}\subseteq {\mathbb{R}}_{+}^m\end{array}$$

u∈Umax​x∈F(y,u)min​s.t.​bTxGx⩾h−Ey−Mu,u∈USx​⊆R+m​​

这部分原文是这么描述的

Consider the case where the second-stage decision problem is a linear programming (LP) problem in $\mathbf{x}$. We first take the relatively complete recourse assumption that this LP is feasible for any given $\mathbf{y}$ and $u$. Let $\pi$ be its dual variables. Then, we obtain its dual problem, which is a maximization problem and can be merged with the maximization over $u$. As a result, we have the following problem, which yields the subproblem in the Benders-dual method.

由于上式中，外层是${\max }_{u\in \mathcal{U}}$，外层的决策变量是不确定变量$u$，因此对于内层，${\min }_{\mathbf{x}\in \mathbf{F}(\mathbf{y},u)}$而言，$\mathbf{y}$和$u$都是已知(given)的，也就是，inner level的模型其实是 given $u$, given $\mathbf{y}$,

min ⁡ x ∈ F ( y ˉ , u ˉ ) b T x s . t . G x ⩾ h − E y ˉ − M u ˉ , → π S x ⊆ R + m

minx∈F(y¯,u¯)s.t.bTxGx⩾h−Ey¯−Mu¯,→πSx⊆Rm+$$\begin{array}{rl}\underset{\mathbf{x}\in \mathbf{F}(\overline{\mathbf{y}},\overline{u})}{min}& {\mathbf{b}}^T\mathbf{x}\\ s.t.& \mathbf{G}\mathbf{x}⩾\mathbf{h}\mathbf{-}\mathbf{E}\overline{\mathbf{y}}\mathbf{-}\mathbf{M}\overline{u},\to \pi \\ & {\mathbf{S}}_{\mathbf{x}}\subseteq {\mathbb{R}}_{+}^m\end{array}$$

x∈F(yˉ​,uˉ)min​s.t.​bTxGx⩾h−Eyˉ​−Muˉ,→πSx​⊆R+m​​
其中$\bar{\mathbf{y}}$是求解first-stage的模型得到的解，而$\bar{u}$是在second-stage的outer-level中已经fix的不确定的量的取值。另外，$\pi$是一个列向量，并不是标量。因此，$h-E\bar{y}-\mathbf{M}\bar{u}$是已知常数列向量，是约束的右端常数项。

由于上述inner level的模型是线性规划，所以我们直接写出上述模型的对偶问题，为

max ⁡ π ( h − E y ˉ − M u ˉ ) T π s . t . G T π ⩽ b , π ⩾ 0.

maxπs.t.(h−Ey¯−Mu¯)TπGTπ⩽b,π⩾0.$$\begin{array}{rl}\underset{\pi }{max}& (\mathbf{h}\mathbf{-}\mathbf{E}\overline{\mathbf{y}}\mathbf{-}\mathbf{M}\overline{u}{)}^T\pi \\ s.t.& {\mathbf{G}}^{\mathbf{T}}\pi ⩽\mathbf{b},\\ & \pi ⩾\mathbf{0}.\end{array}$$

πmax​s.t.​(h−Eyˉ​−Muˉ)TπGTπ⩽b,π⩾0.​
注意，$\pi$是一个列向量，并不是标量。以及$0$是一个元素全部为0的列向量。

如果我们把second-stage的outer-level和inner level的模型写在一起，$\max \max$就等价于$\max$。 于是我们就将第二阶段bi-level（双层）优化模型转化成了一个单层优化模型，也就是把

max ⁡ u max ⁡ π ( h − E y ˉ − M u ) T π s . t . G T π ⩽ b , π ⩾ 0 u ∈ U .

maxumaxπs.t.(h−Ey¯−Mu)TπGTπ⩽b,π⩾0u∈U.$$\begin{array}{rl}\underset{u}{max}\underset{\pi }{max}& (\mathbf{h}\mathbf{-}\mathbf{E}\overline{\mathbf{y}}\mathbf{-}\mathbf{M}u{)}^T\pi \\ s.t.& {\mathbf{G}}^{\mathbf{T}}\pi ⩽\mathbf{b},\\ & \pi ⩾\mathbf{0}\\ & u\in \mathcal{U}.\end{array}$$

umax​πmax​s.t.​(h−Eyˉ​−Mu)TπGTπ⩽b,π⩾0u∈U.​

等价为

max ⁡ π , u ( h − E y ˉ − M u ) T π s . t . G T π ⩽ b , π ⩾ 0 u ∈ U .

maxπ,us.t.(h−Ey¯−Mu)TπGTπ⩽b,π⩾0u∈U.$$\begin{array}{rl}\underset{\pi ,u}{max}& (\mathbf{h}\mathbf{-}\mathbf{E}\overline{\mathbf{y}}\mathbf{-}\mathbf{M}u{)}^T\pi \\ s.t.& {\mathbf{G}}^{\mathbf{T}}\pi ⩽\mathbf{b},\\ & \pi ⩾\mathbf{0}\\ & u\in \mathcal{U}.\end{array}$$

π,umax​s.t.​(h−Eyˉ​−Mu)TπGTπ⩽b,π⩾0u∈U.​

• 注意，由于是将$\max \max$等价成了$\max$，因此outer-level的不确定变量$u$也就纳入到等价后的一阶段模型中来了。其中$\bar{\mathbf{y}}$是从第一阶段得到的fix的值。

由于上述模型中，$\bar{\mathbf{y}}$已经被固定，因此$(\mathbf{h}-\mathbf{E}\bar{\mathbf{y}}-\mathbf{M}u{)}^T\pi$就为最优解中第二阶段的目标函数提供了一个下界。

为什么$(\mathbf{h}-\mathbf{E}\bar{\mathbf{y}}-\mathbf{M}u{)}^T\pi$就为最优解中第二阶段的目标函数提供了一个下界呢？因为，如果$\mathbf{y}$没有被固定，则${\max }_{\pi ,u}(\mathbf{h}-\mathbf{E}\mathbf{y}-\mathbf{M}u{)}^T\pi ⩾{\max }_{\pi ,u}(\mathbf{h}-\mathbf{E}\bar{\mathbf{y}}-\mathbf{M}u{)}^T\pi$一定成立，因此$(\mathbf{h}-\mathbf{E}\bar{\mathbf{y}}-\mathbf{M}u{)}^T\pi$就为第二阶段的目标函数提供了一个下界。

我们将上述模型简略地写为

S P 1 : Q ( y ) = max ⁡ π , u { ( h − E y ˉ − M u ) T π : G T π ⩽ b , u ∈ U , π ⩾ 0. }

SP1:Q(y)=maxπ,u{(h−Ey¯−Mu)Tπ:GTπ⩽b,u∈U,π⩾0.}$$\begin{array}{r}\mathbf{S}{\mathbf{P}}_{\mathbf{1}}:\mathcal{Q}(\mathbf{y})=\underset{\pi ,u}{max}\{(\mathbf{h}\mathbf{-}\mathbf{E}\overline{\mathbf{y}}\mathbf{-}\mathbf{M}u{)}^T\pi :{\mathbf{G}}^{\mathbf{T}}\pi ⩽\mathbf{b},u\in \mathcal{U},\pi ⩾\mathbf{0}.\}\end{array}$$

SP1​:Q(y)=π,umax​{(h−Eyˉ​−Mu)Tπ:GTπ⩽b,u∈U,π⩾0.}​
这相当于把上面的数学规划模型，写成了一个函数形式，实际上是等价的，只是另一个简便，节省地方的写法。这个部分，也被称为是recourse problem，可以理解为，第一阶段决策确定了，第二阶段等到不确定参数揭示了，我们做第二阶段的决策，来对整个问题进行追索，秋后算账。

写好了子问题的模型$\mathbf{S}{\mathbf{P}}_1$，那接下来的重要工作当然是：求解$\mathbf{S}{\mathbf{P}}_1$。

$\mathbf{S}{\mathbf{P}}_1$如何求解呢？：

注意，$\mathbf{S}{\mathbf{P}}_1$是一个 bilinear optimization problem （因为$u\cdot \pi$是二次项，两个都是决策变量），根据文献[1]中的参考文献，可知：

• 可以用一些启发式算法求解$\mathbf{S}{\mathbf{P}}_1$; 这种方法得到的不一定是$\mathbf{S}{\mathbf{P}}_1$的最优解，因此算法收敛速度也许会慢一些；
• 另外，如果$\mathcal{U}$具有特殊的结构，可以构造精确性算法进行精确求解; 比如，如果$u$是0-1变量，可以进行等价线性化，变成MIP进行求解；
• 当然，根据笔者的经验，还可以使用Gurobi，COPT等求解器来求解$\mathbf{S}{\mathbf{P}}_1$;
• 最后，还可以使用文献[1]介绍的方法：使用KKT条件来求解$\mathbf{S}{\mathbf{P}}_1$。本推文也是使用了KKT条件来求解$\mathbf{S}{\mathbf{P}}_1$的，见下文的超级详细的介绍；

综上，大概可以用四种可以来求解$\mathbf{S}{\mathbf{P}}_1$：

1. 启发式算法
2. $\mathcal{U}$具有特殊的结构，可以构造精确性算法进行精确求解;
3. 使用Gurobi，COPT等求解器来求解；
4. 使用KKT条件来求解。

这部分原文是这么描述的

Note that the resulting problem in (2) is a bilinear optimization problem. Several solution strategies have been developed, either in a heuristic fashion or for instances with specially-structured $\mathcal{U}$.

求解了$\mathbf{S}{\mathbf{P}}_1$之后，又能有什么用呢？

答案是：为主问题提供割平面，从而改进全局界限，改进当前解。

假设在Benders-dual算法中，对于给定的第$k$步迭代中第一阶段（Master Problem）得到的${\mathbf{y}}_k^{\ast }$，我们都能得到一个相应$\mathcal{Q}({\mathbf{y}}_k^{\ast })$的第二阶段的解$(u_k^{\ast },{\pi }_k^{\ast })$，我们就可以生成一个下面形式的割平面：

• 上面已经介绍过，${\max }_{\pi ,u}(\mathbf{h}-\mathbf{E}\bar{\mathbf{y}}-\mathbf{M}u{)}^T\pi$为最优解中第二阶段的目标函数的取值提供了一个下界）

因此，割平面的形式如下：

η ⩾ ( h − E y − M u k ∗ ) T π k ∗

η⩾(h−Ey−Mu∗k)Tπ∗k$$\begin{array}{r}\eta ⩾\mathbf{(}\mathbf{h}\mathbf{-}\mathbf{E}\mathbf{y}\mathbf{-}\mathbf{M}{u}_k^{\ast }{)}^T{\pi }_k^{\ast }\end{array}$$

η⩾(h−Ey−Muk∗​)Tπk∗​​
其中，$\eta$表示最优解中第二阶段的目标函数的取值。这个割平面，其实就可以看做是Benders decomposition中的optimality cut。（比较直观的理解请参考我们往期讲解Benders decomposition的推文中最后的Cut 5和 Cut 7）

这个割平面可以被加入到主问题中，主问题就变化为

M P 1 : min ⁡ y , η c T y + η s . t . A y ⩾ d η ⩾ ( h − E y − M u l ∗ ) T π l ∗ , ∀ l ⩽ k y ∈ S y , η ∈ R .

MP1:s.t.miny,ηcTy+ηAy⩾dη⩾(h−Ey−Mu∗l)Tπ∗l,∀l⩽ky∈Sy,η∈R.$$\begin{array}{rl}\mathbf{M}{\mathbf{P}}_{\mathbf{1}}:& \underset{\mathbf{y},\eta }{min}{\mathbf{c}}^{\mathbf{T}}\mathbf{y}\mathbf{+}\eta \\ s.t.& \mathbf{A}\mathbf{y}⩾\mathbf{d}\\ & \eta ⩾\mathbf{(}\mathbf{h}\mathbf{-}\mathbf{E}\mathbf{y}\mathbf{-}\mathbf{M}{u}_l^{\ast }{)}^T{\pi }_l^{\ast },\mathrm{\forall }l⩽k\\ & \mathbf{y}\in {\mathbf{S}}_{\mathbf{y}},\\ & \eta \in \mathbb{R}.\end{array}$$

MP1​:s.t.​y,ηmin​cTy+ηAy⩾dη⩾(h−Ey−Mul∗​)Tπl∗​,∀l⩽ky∈Sy​,η∈R.​
注意，第2条约束中，$\forall l⩽k$指的是在前$k$次的迭代中，每一次产生的割平面都被加入到了$\mathbf{M}{\mathbf{P}}_1$中。并且$\eta$相当于一个辅助变量 ，是为了识别worst-case下的second-stage产生的成本的(也叫作recourse)。

上述模型，我们可以计算到目前为止得到的当前最优解$({\mathbf{y}}_{\mathbf{k}+1}^{\ast },{\eta }_{k+1}^{\ast })$， 这也就是第$k+1$次迭代的时候二者的取值。

注意到，随着割平面 $\eta ⩾(\mathbf{h}-\mathbf{E}\mathbf{y}-\mathbf{M}{u}_l^{\ast }{)}^T{\pi }_l^{\ast },\forall l⩽k$的不断添加，$\mathbf{M}{\mathbf{P}}_1$中$\eta$的值会不断增加(至少是不减)，因此$\mathbf{M}{\mathbf{P}}_1$的目标函数会不断上升。下面会讲到，$\mathbf{M}{\mathbf{P}}_1$的目标函数是全局最优解的一个下界。

与传统的Benders Decomposition算法中相同，Benders-dual算法迭代过程中，也可以计算任意一步迭代过程中，全局的上界和下界的变化，以此来判断当前解的质量，以及是否达到了最优解。

文献[1]中原文对Benders-dual method下的主问题的表述如下：

注意：文中提到了一些非常关键的信息，即Benders-dual method迭代过程中，全局上界和下界的计算。根据原文：

• ${\mathbf{c}}^{\mathbf{T}}{\mathbf{y}}_{\mathbf{k}}^{\ast }+\mathcal{Q}({\mathbf{y}}_{\mathbf{k}}^{\ast })$ 为原始问题目标值 Obj = { min ⁡ y c T y + max ⁡ u ∈ U min ⁡ x ∈ F ( y , u ) b T x ∣ A y ⩾ d , y ∈ S y } \text{Obj} = \{\min_{\mathbf{y}} \mathbf{c^Ty} +\max_{u \in \mathcal{U}} \min_{ \mathbf{x}\in {\mathbf{F}(\mathbf{y}, u )}} \mathbf{b}^T \mathbf{x} | \mathbf{Ay} \geqslant \mathbf{d}, \mathbf{y} \in \mathbf{S}_{\mathbf{y}}\} Obj={miny​cTy+maxu∈U​minx∈F(y,u)​bTx∣Ay⩾d,y∈Sy​} （即原始问题）提供了一个上界；
• 而$\mathbf{M}{\mathbf{P}}_1:{\mathbf{c}}^{\mathbf{T}}{\mathbf{y}}_{\mathbf{k}+1}^{\ast }+{\eta }_{k+1}^{\ast }$为原始问题提供了一个下界。
• 在Benders-dual method算法迭代的过程中，上界和下界会得到不断改进。并且，循环地求解主问题$\mathbf{M}{\mathbf{P}}_1$，得到固定的$\bar{\mathbf{y}}$，并将$\bar{\mathbf{y}}$传给子问题$\mathbf{S}{\mathbf{P}}_1$，更新子问题$\mathbf{S}{\mathbf{P}}_1$，然后求解$\mathbf{S}{\mathbf{P}}_1$，根据$\mathbf{S}{\mathbf{P}}_1$的解，构造割平面$\eta ⩾(\mathbf{h}-\mathbf{E}\mathbf{y}-\mathbf{M}{u}_k^{\ast }{)}^T{\pi }_k^{\ast }$，并且将其添加到$\mathbf{M}{\mathbf{P}}_1$中。
• 同时，更新全局的上界和下界。循环这个过程，最终上下界就可以在有限步之内收敛。

总结一句，Benders-dual算法迭代过程中，任意一步迭代中的全局的上界和下界可以用下面的公式计算得出：

L B = max ⁡ { c T y k + 1 ∗ + η k + 1 ∗ , L B } U B = min ⁡ { c T y k ∗ + Q ( y k ∗ ) , U B }

LB=max{cTy∗k+1+η∗k+1,LB}UB=min{cTy∗k+Q(y∗k),UB}$$\begin{array}{rl}& LB=max\{{\mathbf{c}}^{\mathbf{T}}{\mathbf{y}}_{\mathbf{k}\mathbf{+}\mathbf{1}}^{\mathbf{\ast }}+{\eta }_{k+1}^{\ast },LB\}\\ & UB=min\{{\mathbf{c}}^{\mathbf{T}}{\mathbf{y}}_{\mathbf{k}}^{\mathbf{\ast }}\mathbf{+}\mathcal{Q}\mathbf{(}{\mathbf{y}}_{\mathbf{k}}^{\mathbf{\ast }}\mathbf{)},UB\}\end{array}$$

​LB=max{cTyk+1∗​+ηk+1∗​,LB}UB=min{cTyk∗​+Q(yk∗​),UB}​

这个很好理解：

• $LB={\mathbf{c}}^{\mathbf{T}}{\mathbf{y}}_{\mathbf{k}+1}^{\ast }+{\eta }_{k+1}^{\ast }$，是只考虑了一部分随机场景，松弛了一部分约束，因此对其求$\min$，一定是更小(或者等于最小)，因此它提供了全局最优解的一个下界；直观理解就是少考虑了约束，又是min问题，因此是个下界；
• $UB={\mathbf{c}}^{\mathbf{T}}{\mathbf{y}}_{\mathbf{k}}^{\ast }+\mathcal{Q}({\mathbf{y}}_{\mathbf{k}}^{\ast })$，是相当于将原问题拆分成两个相互独立的问题进行分别求解，然后将各自的目标函数相加，分开优化的可行解，一定也是联合优化的可行解；而分开优化的解，对于联合优化而言，只是一个可行解而已，并不一定是最优解，联合优化的解，一定好于或者等于分开优化的解。因此${\mathbf{c}}^{\mathbf{T}}{\mathbf{y}}_{\mathbf{k}}^{\ast }+\mathcal{Q}({\mathbf{y}}_{\mathbf{k}}^{\ast })$为Two-stage RO原模型的最优解提供了一个上界。
• 当算法迭代过程中，达到$UB=LB$，则算法终止，得到了最优解。

注意到，上述算法描述中，${\pi }_k^{\ast }$和$u_k^{\ast }$都是在他们各自的可行域中的极点(或者离散的点)。这样一来，Benders-dual算法的复杂度如下：