两阶段鲁棒优化及列和约束生成算法_列与约束生成算法

1. 前言

有同学私信我两阶段鲁棒优化的问题，自己之前主要研究单阶段的鲁棒优化，对于两阶段优化不太懂。查了点资料，通过翻译和自己的理解，写下这篇博文，抛砖引玉，以供大家共同学习交流。

2. 两阶段RO

本文主要翻译自2013年Bo Zeng 博士的高被引论文《Solving two-stage robust optimization problems using a column-and-constraint generation method》[1]。

由于传统单阶段的鲁棒优化对于所有不确定性是完全免疫的，求出的结果一般过于保守、比较悲观。两阶段RO也称鲁棒可调优化或者自适应优化，它的引入就是为了应对传统RO存在的上述问题。两阶段RO是指，在作出第一阶段的决策和不确定性部分展现出来之后，再进行第二阶段的决策。 由于改进的建模能力，两阶段RO，广泛应用于网络/运输问题、投资组合优化和电力系统调度等问题。

文献[1]假设一阶段和二阶段决策问题都是线性规划，并且不确定性集合$\mathbf{U}$是离散有限的点集或者多面体。使用$\mathbf{y}$表示第一阶段决策变量，$\mathbf{x}$表示第二阶段决策变量，$u\in \mathbf{U}$表示不确定矢量。在此假设下的两阶段RO的一般形式为：
$$\underset{\mathbf{y}}{\min }{\mathbf{c}}^{\mathrm{T}}\mathbf{y}+\underset{\mathrm{u}\in \mathbf{U}}{\max }\underset{\mathbf{x}\in \mathbf{F}(\mathbf{y},\mathrm{u})}{\min }{\mathbf{b}}^{\mathrm{T}}\mathbf{x}\text{\hspace{1em}(1)}$$ $$\mathrm{s}.\mathrm{t}.{\mathbf{A}}^{\mathrm{T}}\mathbf{y}\ge \mathbf{d},\mathbf{y}\in {\mathbf{S}}_{\mathbf{y}}$$

其中， F ( y , u ) = { x : G x ≥ h − E y − M u , x ∈ S x } \bf F(y,\rm u)=\{\bf x: Gx \geq h-Ey-M \rm u,\bf x \in S_{x}\} F(y,u)={x:Gx≥h−Ey−Mu,x∈Sx​}，${\mathbf{S}}_{\mathbf{y}}\subseteq {\mathbb{R}}_{+}^{\mathbf{n}}$，${\mathbf{S}}_{\mathbf{x}}\subseteq {\mathbb{R}}_{+}^{\mathbf{m}}$，向量$\mathbf{c},\mathbf{b},\mathbf{d},\mathbf{h}$和矩阵$\mathbf{A},\mathbf{G},\mathbf{E},\mathbf{M}$都是确定性的数值，不确定性体现在向量$u$上。注意到第二阶段优化的约束条件$\mathbf{F}(\mathbf{y},\mathrm{u})$是关于不确定性$u$的线性函数。

两阶段RO有优势，当然有缺点。即便是最简单的两阶段RO，也可以是NP-hard问题，计算复杂度很高。为了减轻计算负担，一般有两种方法。第一种是使用近似算法，该种方法假设第二阶段的决策是关于不确定性的简单函数，例如仿射函数。第二种方法试图根据Benders分解得出精确解，即它们使用第二阶段决策问题的对偶解，逐步构造第一阶段决策的值函数（value function）。一般称为Benders对偶割平面算法。

3. Benders对偶割平面法

由于第二阶段的决策是关于$\mathbf{x}$的线性规划问题，可以作以下假设（relatively complete recourse assumption）：该线性规划对于任意给定的$\mathbf{y}$和$u$是可行的，也即是有解 (该假设不是很理解)。假设第二阶段的线性规划的对偶变量为 $\pi$，则将其转化为对偶问题为：
S P 1 : O ( y ) = max ⁡ u , π { ( h − E y − M u ) T π } (2) \rm {SP_{1}}: \mathcal O(\bf y) =\rm \max_{u,\pi} \{ \bf (h-Ey-M \rm u)^{T}\pi\} \tag{2} SP1​:O(y)=u,πmax​{(h−Ey−Mu)Tπ}(2) $$s.t.{\mathbf{G}}^{\mathrm{T}}\pi \le \mathbf{b},\mathrm{u}\in \mathbf{U},\pi \ge 0$$

在对偶问题中，目标函数从原始的最小化转换为关于对偶变量$\pi$的最大化，同时与(1)式中的最大化$u$合并，得到上述子问题(2)。此时不确定性向量转化为对偶问题(2)的决策变量，需要注意到子问题(2)是存在$u^T\pi$的双线性项。此时问题(2)被称为关于$u,\pi$的双线性规划。对于双线性规划的求解，放在后面再说。

假设对于给定的${\mathbf{y}}_k^{\ast }$，子问题(2) 的最优解为 $(u_k^{\ast },{\pi }_k^{\ast })$，根据对偶定理，则可以构建以下割平面：
$$\eta \ge (\mathbf{h}-\mathbf{E}\mathbf{y}-\mathbf{M}{\mathrm{u}}_k^{\ast }{)}^{\mathrm{T}}{\pi }_k^{\ast }\text{\hspace{1em}(3)}$$
其中，$\eta ={\max }_{u\in \mathbf{U}}{\min }_{\mathbf{x}\in \mathbf{F}(\mathbf{y},\mathrm{u})}{\mathbf{b}}^{\mathrm{T}}\mathbf{x}$为一维标量。注意该割平面是关于第一阶段决策变量$\mathbf{y}$的约束。

将割平面约束添加第一阶段优化中，可以得到：
$$\mathrm{M}{\mathrm{P}}_1:\underset{\mathbf{y}}{\min }{\mathbf{c}}^{\mathrm{T}}\mathbf{y}+\eta \text{\hspace{1em}(4)}$$ $$\mathrm{s}.\mathrm{t}.{\mathbf{A}}^{\mathrm{T}}\mathbf{y}\ge \mathbf{d},\mathbf{y}\in {\mathbf{S}}_{\mathbf{y}}$$ $$\eta \ge (\mathbf{h}-\mathbf{E}\mathbf{y}-\mathbf{M}{\mathrm{u}}_l^{\ast }{)}^{\mathrm{T}}{\pi }_l^{\ast },\forall l\le k$$ $$\mathbf{y}\in {\mathbf{S}}_{\mathbf{y}},\eta \in \mathbb{R}$$

其中$(u_l^{\ast },{\pi }_l^{\ast })$已知。求解主问题(4)可以获得其最优解$({\mathbf{y}}_{k+1}^{\ast },{\eta }_{k+1}^{\ast })$。

此时问题(4)的目标函数值${\mathbf{c}}^{\mathrm{T}}{\mathbf{y}}_{k+1}^{\ast }+{\eta }_{k+1}^{\ast }$是原始两阶段问题(1)的下界(lower bound, LB)。如何理解？原始的两阶段RO是在worse-case情况下的目标值，也即是对于所有不确定性都具有包容性，可以理解成其最优决策在所有的不确定性下都成立。通过把不确定性转化为割平面，主问题(4)中添加的割平面只是在部分不确定性$u_l$下的优化，也即是问题(4)是关于问题(1)的松弛问题，由于是求最小化，问题(4)的最优目标值一定是原问题(1)的一个下界。

注意到${\mathbf{c}}^{\mathrm{T}}{\mathbf{y}}_k^{\ast }+\mathcal{O}({\mathbf{y}}_k^{\ast })$是原始问题(1)的上界。如何理解？${\mathbf{y}}_k^{\ast }$是第一阶段的一个可行解，$\mathcal{O}({\mathbf{y}}_k^{\ast })$是第二阶段优化在${\mathbf{y}}_k^{\ast }$下的最优目标值，也即是说存在二阶段的最优解${\mathbf{x}}_k^{\ast }$。此时的${\mathbf{y}}_k^{\ast },{\mathbf{x}}_k^{\ast }$只是问题(1)的可行解，而不一定是最优解。由于是求最小化问题，${\mathbf{c}}^{\mathrm{T}}{\mathbf{y}}_k^{\ast }+\mathcal{O}({\mathbf{y}}_k^{\ast })$是原始问题(1)的上界。

把${\mathbf{y}}_{k+1}^{\ast }$带入到子问题(2)中，求解问题(2)得到最优解$(u_{k+1}^{\ast },{\pi }_{k+1}^{\ast })$ 和最优值$\mathcal{O}({\mathbf{y}}_{k+1}^{\ast })$，进而可以构造新的割平面。因此，迭代的引入割平面(3)，计算主问题(1)，上届和下界逐渐收敛到问题(1)的最优解。

计算复杂度[1]
命题1：如果不确定集合$\mathbf{U}$是多面体，用$p$表示极点的数量，如果$\mathbf{U}$是离散点集，用$p$表示集合的势（集合中元素的数量）。用$q$表示多面体 { π : G T π ≤ b , π ≥ 0 } \{\pi: \bf G^{\rm T}\pi \leq b, \pi \geq 0 \} {π:GTπ≤b,π≥0} 的极点数量。Benders对偶算法需要经过$o(pq)$次迭代才能求出问题(1)的最优解。

Benders 对偶割平面法主要的问题是求解问题(2)的双线性规划问题。

4. 列和约束生成(C&CG)算法

列和约束生成(column-and-constraint generation) 算法是基于以下事实：
假设${\mathbf{x}}^l$是不确定性$u=u_l$(一个实例)下$\eta ={\max }_{u\in \mathbf{U}}{\min }_{\mathbf{x}\in \mathbf{F}(\mathbf{y},\mathrm{u})}{\mathbf{b}}^{\mathrm{T}}\mathbf{x}$的最优解，则一定存在如下约束成立：
$$\eta \ge {\mathbf{b}}^{\mathrm{T}}{\mathbf{x}}^l$$ $$\mathbf{G}{\mathbf{x}}^l\ge \mathbf{h}-\mathbf{E}\mathbf{y}-\mathbf{M}{\mathrm{u}}_l$$
可以根据以上约束构造割平面，形成C&CG算法。

– Step 1 ：设置$\mathrm{L}\mathrm{B}=-\infty$，$\mathrm{U}\mathrm{B}=+\infty$，索引$k=0$，集合$\mathbf{O}=\varnothing$
– Step 2： 求解如下主问题：
$$\mathrm{M}{\mathrm{P}}_2:\underset{\mathbf{y}}{\min }{\mathbf{c}}^{\mathrm{T}}\mathbf{y}+\eta \text{\hspace{1em}(5)}$$ $$\mathrm{s}.\mathrm{t}.{\mathbf{A}}^{\mathrm{T}}\mathbf{y}\ge \mathbf{d},\mathbf{y}\in {\mathbf{S}}_{\mathbf{y}}$$ $$\eta \ge {\mathbf{b}}^{\mathrm{T}}{\mathbf{x}}^l,\forall l\in \mathbf{O}$$ $$\mathbf{G}{\mathbf{x}}^l\ge \mathbf{h}-\mathbf{E}\mathbf{y}-\mathbf{M}{\mathrm{u}}_l^{\ast }\forall l\le k$$ $$\mathbf{y}\in {\mathbf{S}}_{\mathbf{y}},\eta \in \mathbb{R},{\mathbf{x}}^l\in {\mathbf{S}}_{\mathbf{x}}\forall l\le k$$
求出最优解$({\mathbf{y}}_{k+1}^{\ast },{\eta }_{k+1}^{\ast },{\mathbf{x}}^{1\ast },\cdots ,{\mathbf{x}}^{k\ast })$，更新下界$\mathrm{L}\mathrm{B}={\mathbf{c}}^{\mathrm{T}}{\mathbf{y}}_{k+1}^{\ast }+{\eta }_{k+1}^{\ast }$

– Step 3 ：代入$\mathbf{y}={\mathbf{y}}_{k+1}^{\ast }$，求解如下子问题：
S P 2 : O ( y ) = max ⁡ u ∈ U min ⁡ x ∈ F ( y , u ) { b T x : G x ≥ h − E y − M u , x ∈ S x } (6) \rm SP2: \mathcal O(\bf y)= \rm \max_{u \in \bf U} \min_{\bf x \in F(y,\rm u)} \{ \bf b^{\rm T}x :\bf Gx \geq h-Ey-M \rm u, \bf x \in S_{x} \} \tag{6} SP2:O(y)=u∈Umax​x∈F(y,u)min​{bTx:Gx≥h−Ey−Mu,x∈Sx​}(6)
更新上界 U B = min ⁡ { U B , c T y k + 1 ∗ + O ( y k + 1 ∗ ) } \rm UB=\min\{UB,\bf c^{\rm T}y_{\mathit k+1}^{*}+\mathcal O(\bf y_{\mathit k+1}^{*})\} UB=min{UB,cTyk+1∗​+O(yk+1∗​)}。

– Step 4 ： 如果$\mathrm{U}\mathrm{B}-\mathrm{L}\mathrm{B}\le ϵ$, 返回${\mathbf{y}}_{k+1}^{\ast }$，程序终止。否则：
(a) 如果$\mathcal{O}({\mathbf{y}}_{k+1}^{\ast })<+\infty$，添加变量${\mathbf{x}}^{k+1}$，添加如下约束：
$$\eta \ge {\mathbf{b}}^{\mathrm{T}}{\mathbf{x}}^{k+1}\text{\hspace{1em}(7)}$$ $$\mathbf{G}{\mathbf{x}}^{k+1}\ge \mathbf{h}-\mathbf{E}\mathbf{y}-\mathbf{M}{\mathrm{u}}_{k+1}^{\ast }\text{\hspace{1em}(8)}$$
到问题(5)中。其中$u_{k+1}^{\ast }$是问题(6)在$\mathbf{y}={\mathbf{y}}_{k+1}^{\ast }$下的最优场景(不确定性)，对问题(6)可以利用数据库（Call the oracle to solve subproblem SP2）进行枚举求得。
然后，更新$k=k+1$，集合$\mathbf{O}=\mathbf{O}\bigcup k+1$，跳转至步骤Step 2。

(b) 如果$\mathcal{O}({\mathbf{y}}_{k+1}^{\ast })=+\infty$ (对于某些$u^{\ast }\in \mathbf{U}$，如果第二阶段决策$\mathcal{O}({\mathbf{y}}_{k+1}^{\ast })$不可行(infeasible)，则把$\mathcal{O}({\mathbf{y}}_{k+1}^{\ast })$记为$+\infty$)，添加变量${\mathbf{x}}^{k+1}$，添加如下约束：
$$\mathbf{G}{\mathbf{x}}^{k+1}\ge \mathbf{h}-\mathbf{E}\mathbf{y}-\mathbf{M}{\mathrm{u}}_{k+1}^{\ast }\text{\hspace{1em}(9)}$$ 到问题(5)中。其中$u_{k+1}^{\ast }$是问题(6)在$\mathbf{y}={\mathbf{y}}_{k+1}^{\ast }$下不可行的不确定性$u$的值。
然后，更新$k=k+1$，跳转至步骤Step 2。

对于步骤Step 4(a)的约束(7)和(8)被称为最优割(optimality cuts)，而Step 4(b)的约束(9)被称为可性割(feasibility cuts)。对于可行割不是很明白，虽然在某些$\mathbf{y}={\mathbf{y}}_{k+1}^{\ast }$下不可行，但是仍然是其中一个不确定性$u_{k+1}^{\ast }$对二阶段决策的一个反映，因为二阶段是包含所有的不确定性。 不知以上理解是否正确。 由于约束(7)对于不可行的场景仍然是成立的，此时的步骤4(a)和4(b)可以合并，因此可以用统一的方法和程序去处理最优性和可行性。

由于生成的割平面是由第二阶段决策变量（wait-and-see decision variables/recourse decision variables）以二阶段决策约束的形式定义的，整个过程其实就是列和约束生成，其中列生产是指在主问题中添加第二阶段决策变量，约束生成是指添加割平面约束。

C&CG算法的复杂性
命题2：如果不确定集合$\mathbf{U}$是多面体，用$p$表示极点的数量，如果$\mathbf{U}$是离散点集，用$p$表示集合的势（集合中元素的数量）。C&CG算法需要经过$o(p)$次迭代才能求出问题(1)的最优解。

C&CG算法与Benders对偶的区别：
– (1) 主问题中决策的决策变量。C&CG算法每次迭代都添加决策变量${\mathbf{x}}^l$，而Benders对偶每次迭代决策变量不变。

– (2)计算复杂度。由以上两个命题可以看出，在( the relatively complete recourse assumption)假设下，C&CG算法具有更低的复杂度。

– (3)求解问题的能力。Benders分解要求第二阶段的决策问题为行性规划问题，C&CG算法对于第二阶段优化问题的变量的类型不敏感，可以是混合整数规划问题。

5 存在的难点

以上两种算法存在的难点在于：Benders对偶算法求解子问题(2)的双线性规划；C&CG算法求解子问题(6)。

针对相对简单的多面体不确定性集，一般使用外部近似算法(outer approximation algorithm)和混合整数线性重构(mixed integer linear reformulation)两种方法求解。 第一种是求解SP1的启发式算法。 第二种是使用不确定性集的特殊结构，将双线性规划SP1转换为等效的混合整数线性规划。 文献[1] 使用经典的Karush–Kuhn–Tucker（KKT）条件来处理一般的多面体不确定性集，只要( the relatively complete recourse assumption)假设成立即可。文献[1]对于问题(6)SP2，通过引入对偶变量、利用KKT、结合大M方法，将鲁棒优化问题(6)转化为确定性的整数线性规划问题，再利用现成的求解器进行求解。具体转化就不翻译了，自己没看懂。

文献[2]中表述：根据文献[3]的结论，问题(2)的双线性规划，变量$u$的最优解为不确定集合$U$的极点。对于多面体的不确定集合，其极点是有限的，需要找出所有的极点，然后通过枚举就可以求出最优的$u$。同样根据文献[4]，在问题(2)达到最优时，变量$u$ 和对偶变量$\pi$取值总是其各自可行集的极点，因此也需要求出$\pi$的极点。如果变量$u$ 和$\pi$的维数比较大的话，枚举的复杂度仍然比较高。

对于问题(6)本实质上是单阶段的RO问题，可以使用求解一般RO的求解器进行求解。但是它们与传统的RO仍有区别，不论是问题(2)还是(6)都需要求出具体的不确定性$u$的值。可以认为，传统的线性RO其最大的不确定性，在不确定性集合的极点处取得，因此还是需要计算不确定集合$U$的极点。先使用单阶段RO求解出决策变量${\mathbf{x}}^l$，然后在枚举所有的$U$的极点，找到满足临界约束条件的点，应该就是具体的$u$。不知理解是否正确，后面有时间再写代码吧。

参考文献
[1] 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.
[2]刘一欣, 郭力, 王成山. 微电网两阶段鲁棒优化经济调度方法[J]. 中国电机工程学报, 2018, 038(014):4013-4022.
[3] D. Bertsimas, E. Litvinov, X.A. Sun, Jinye Zhao, Tongxin Zheng, Adaptive robust optimization for the security constrained unit commitment problem, IEEE Transactions on Power Systems 28 (1) (2013) 52–63.
[4] Bo Zeng. Solving Two-stage Robust Optimization Problems by A Constraint-and-Column Generation Method. 2011. http://www.optimization-online.org/DB_FILE/2011/06/3065.pdf