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

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

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

• 1. 前言
• 2. 两阶段RO
• 3. Benders对偶割平面法
• 4. 列和约束生成(C&CG)算法
• 5 存在的难点

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]假设一阶段和二阶段决策问题都是线性规划，并且不确定性集合 U bf U U是离散有限的点集或者多面体。使用 y bf y y表示第一阶段决策变量， x bf x x表示第二阶段决策变量， u ∈ U u in bf U u∈U表示不确定矢量。在此假设下的两阶段RO的一般形式为：
min ⁡ y c T y + max ⁡ u ∈ U min ⁡ x ∈ F ( y , u ) b T x (1) min_{bf y } bf c^{rm T}y+rm max_{u in bf U} min_{bf x in F(y,rm u)} bf b^{rm T}x tag{1} ymin​cTy+u∈Umax​x∈F(y,u)min​bTx(1) s . t . A T y ≥ d , y ∈ S y quad bf A^{rm T}y geq d,y in S_{y} s.t.ATy≥d,y∈Sy​

其中， 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​}， S y ⊆ R + n bf S_{y}subseteq Bbb R_{+}^{n} Sy​⊆R+n​， S x ⊆ R + m bf S_{x}subseteq Bbb R_{+}^{m} Sx​⊆R+m​，向量 c , b , d , h bf c,b,d,h c,b,d,h和矩阵 A , G , E , M bf A,G,E,M A,G,E,M都是确定性的数值，不确定性体现在向量 u u u上。注意到第二阶段优化的约束条件 F ( y , u ) bf F(y,rm u) F(y,u)是关于不确定性 u u u的线性函数。

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

3. Benders对偶割平面法

由于第二阶段的决策是关于 x bf x x的线性规划问题，可以作以下假设（relatively complete recourse assumption）：该线性规划对于任意给定的 y bf y y和 u u 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 . G T π ≤ b , u ∈ U , π ≥ ; bf G^{rm T}pi leq b, rm u in bf U, pi geq bf GTπ≤b,u∈U,π≥0

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

假设对于给定的 y k ∗ bf y_{mathit k}^{*} yk∗​，子问题(2) 的最优解为 ( u k ∗ , π k ∗ ) (u_{ k}^{*},pi_{ k}^{*}) (uk∗​,πk∗​)，根据对偶定理，则可以构建以下割平面：
η ≥ ( h − E y − M u k ∗ ) T π k ∗ (3) eta geq bf (h-Ey-M rm u_{mathit k}^{*})^{T}pi_{ mathit k}^{*} tag{3} η≥(h−Ey−Muk∗​)Tπk∗​(3)
其中， η = max ⁡ u ∈ U min ⁡ x ∈ F ( y , u ) b T x eta = max_{u in bf U} min_{bf x in F(y,rm u)} bf b^{rm T}x η=maxu∈U​minx∈F(y,u)​bTx为一维标量。注意该割平面是关于第一阶段决策变量 y bf y y的约束。

将割平面约束添加第一阶段优化中，可以得到：
M P 1 : min ⁡ y c T y + η (4) rm MP_{1}: min_{bf y } bf c^{rm T}y+ eta tag{4} MP1​:ymin​cTy+η(4) s . t . A T y ≥ d , y ∈ S y quad bf A^{rm T}y geq d,y in S_{y} s.t.ATy≥d,y∈Sy​ η ≥ ( h − E y − M u l ∗ ) T π l ∗ , ∀ l ≤ k eta geq bf (h-Ey-M rm u_{mathit l}^*)^{T}pi_{ mathit l}^* ,forall ; mathit {l leq k} η≥(h−Ey−Mul∗​)Tπl∗​,∀l≤k y ∈ S y , η ∈ R bf y in S_{y}, eta in Bbb R y∈Sy​,η∈R

其中 ( u l ∗ , π l ∗ ) (u_{ l}^{*},pi_{ l}^{*}) (ul∗​,πl∗​)已知。求解主问题(4)可以获得其最优解 ( y k + 1 ∗ , η k + 1 ∗ ) (bf y_{mathit k+1}^{*},eta_{mathit k+1}^{*}) (yk+1∗​,ηk+1∗​)。

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

注意到 c T y k ∗ + O ( y k ∗ ) bf c^{rm T}y_{mathit k}^{*}+mathcal O(bf y_{mathit k}^{*}) cTyk∗​+O(yk∗​)是原始问题(1)的上界。如何理解？ y k ∗ bf y_{mathit k}^{*} yk∗​是第一阶段的一个可行解， O ( y k ∗ ) mathcal O(bf y_{mathit k}^{*}) O(yk∗​)是第二阶段优化在 y k ∗ bf y_{mathit k}^{*} yk∗​下的最优目标值，也即是说存在二阶段的最优解 x k ∗ bf x_{mathit k}^* xk∗​。此时的 y k ∗ , x k ∗ bf y_{mathit k}^{*},bf x_{mathit k}^{*} yk∗​,xk∗​只是问题(1)的可行解，而不一定是最优解。由于是求最小化问题， c T y k ∗ + O ( y k ∗ ) bf c^{rm T}y_{mathit k}^{*}+mathcal O(bf y_{mathit k}^{*}) cTyk∗​+O(yk∗​)是原始问题(1)的上界。

把 y k + 1 ∗ bf y_{mathit k+1}^{*} yk+1∗​带入到子问题(2)中，求解问题(2)得到最优解 ( u k + 1 ∗ , π k + 1 ∗ ) (u_{ k+1}^{*},pi_{ k+1}^{*}) (uk+1∗​,πk+1∗​) 和最优值 O ( y k + 1 ∗ ) mathcal O(bf y_{mathit k+1}^{*}) O(yk+1∗​)，进而可以构造新的割平面。因此，迭代的引入割平面(3)，计算主问题(1)，上届和下界逐渐收敛到问题(1)的最优解。

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

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

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

列和约束生成(column-and-constraint generation) 算法是基于以下事实：
假设 x l bf x^{mathit l} xl是不确定性 u = u l u=u_{l} u=ul​(一个实例)下 η = max ⁡ u ∈ U min ⁡ x ∈ F ( y , u ) b T x eta = max_{u in bf U} min_{bf x in F(y,rm u)} bf b^{rm T}x η=maxu∈U​minx∈F(y,u)​bTx的最优解，则一定存在如下约束成立：
η ≥ b T x l eta geq bf b^{rm T}x^{mathit l} η≥bTxl G x l ≥ h − E y − M u l bf Gx^{mathit l} geq h-Ey-M rm u_{mathit l} Gxl≥h−Ey−Mul​
可以根据以上约束构造割平面，形成C&CG算法。

– Step 1 ：设置 L B = − ∞ rm LB=-infty LB=−∞， U B = + ∞ rm UB=+infty UB=+∞，索引 k = 0 k=0 k=0，集合 O = ∅ bf O=emptyset O=∅
– Step 2： 求解如下主问题：
M P 2 : min ⁡ y c T y + η (5) rm MP_{2}: min_{bf y } bf c^{rm T}y+ eta tag{5} MP2​:ymin​cTy+η(5) s . t . A T y ≥ d , y ∈ S y quad bf A^{rm T}y geq d,y in S_{y} s.t.ATy≥d,y∈Sy​ η ≥ b T x l , ∀ l ∈ O eta geq bf b^{rm T}x^{mathit l}, forall mathit {l in bf O} η≥bTxl,∀l∈O G x l ≥ h − E y − M u l ∗ ∀ l ≤ k bf Gx^{mathit l} geq h-Ey-M rm u_{mathit l}^* forall mathit {l leq k} Gxl≥h−Ey−Mul∗​∀l≤k y ∈ S y , η ∈ R , x l ∈ S x ∀ l ≤ k bf y in S_{y}, eta in Bbb R, x^{mathit l} in S_{x} ; forall mathit {l leq k} y∈Sy​,η∈R,xl∈Sx​∀l≤k
求出最优解 ( y k + 1 ∗ , η k + 1 ∗ , x 1 ∗ , ⋯ , x k ∗ ) (bf y_{mathit k+1}^*,eta_{mathit k+1}^*,x^{mathit 1*},cdots,x^{mathit k*}) (yk+1∗​,ηk+1∗​,x1∗,⋯,xk∗)，更新下界 L B = c T y k + 1 ∗ + η k + 1 ∗ rm LB=bf c^{rm T}y_{mathit k+1}^{*}+eta_{mathit k+1}^{*} LB=cTyk+1∗​+ηk+1∗​

– Step 3 ：代入 y = y k + 1 ∗ bf y=y_{mathit k+1}^* y=yk+1∗​，求解如下子问题：
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 ： 如果 U B − L B ≤ ϵ rm UB-LBleq epsilon UB−LB≤ϵ, 返回 y k + 1 ∗ bf y_{mathit k+1}^* yk+1∗​，程序终止。否则：
(a) 如果 O ( y k + 1 ∗ ) < + ∞ mathcal O(bf y_{mathit k+1}^{*}) < +infty O(yk+1∗​)<+∞，添加变量 x k + 1 bf x^mathit {k+1} xk+1，添加如下约束：
η ≥ b T x k + 1 (7) eta geq bf b^{rm T}x^{mathit {k+1}} tag{7} η≥bTxk+1(7) G x k + 1 ≥ h − E y − M u k + 1 ∗ (8) bf Gx^{mathit {k+1}} geq h-Ey-M rm u_{mathit {k+1}}^* tag{8} Gxk+1≥h−Ey−Muk+1∗​(8)
到问题(5)中。其中 u k + 1 ∗ u_{mathit {k+1}}^* uk+1∗​是问题(6)在 y = y k + 1 ∗ bf y=y_{mathit k+1}^* y=yk+1∗​下的最优场景(不确定性)，对问题(6)可以利用数据库（Call the oracle to solve subproblem SP2）进行枚举求得。
然后，更新 k = k + 1 k=k+1 k=k+1，集合 O = O ⋃ k + 1 bf O= O bigcup {mathit {k+1}} O=O⋃k+1，跳转至步骤Step 2。

(b) 如果 O ( y k + 1 ∗ ) = + ∞ mathcal O(bf y_{mathit k+1}^{*}) =+infty O(yk+1∗​)=+∞ (对于某些 u ∗ ∈ U u^*in bf U u∗∈U，如果第二阶段决策 O ( y k + 1 ∗ ) mathcal O(bf y_{mathit k+1}^{*}) O(yk+1∗​)不可行(infeasible)，则把 O ( y k + 1 ∗ ) mathcal O(bf y_{mathit k+1}^{*}) O(yk+1∗​)记为 + ∞ +infty +∞)，添加变量 x k + 1 bf x^mathit {k+1} xk+1，添加如下约束：
G x k + 1 ≥ h − E y − M u k + 1 ∗ (9) bf Gx^{mathit {k+1}} geq h-Ey-M rm u_{mathit {k+1}}^* tag{9} Gxk+1≥h−Ey−Muk+1∗​(9) 到问题(5)中。其中 u k + 1 ∗ u_{mathit {k+1}}^* uk+1∗​是问题(6)在 y = y k + 1 ∗ bf y=y_{mathit k+1}^* y=yk+1∗​下不可行的不确定性 u u u的值。
然后，更新 k = k + 1 k=k+1 k=k+1，跳转至步骤Step 2。

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

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

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

C&CG算法与Benders对偶的区别：
– (1) 主问题中决策的决策变量。C&CG算法每次迭代都添加决策变量 x l bf x^mathit l xl，而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的最优解为不确定集合 U U U的极点。对于多面体的不确定集合，其极点是有限的，需要找出所有的极点，然后通过枚举就可以求出最优的 u u u。同样根据文献[4]，在问题(2)达到最优时，变量 u u u 和对偶变量 π pi π取值总是其各自可行集的极点，因此也需要求出 π pi π的极点。如果变量 u u u 和 π pi π的维数比较大的话，枚举的复杂度仍然比较高。

对于问题(6)本实质上是单阶段的RO问题，可以使用求解一般RO的求解器进行求解。但是它们与传统的RO仍有区别，不论是问题(2)还是(6)都需要求出具体的不确定性 u u u的值。可以认为，传统的线性RO其最大的不确定性，在不确定性集合的极点处取得，因此还是需要计算不确定集合 U U U的极点。先使用单阶段RO求解出决策变量 x l bf x^{mathit l} xl，然后在枚举所有的 U U U的极点，找到满足临界约束条件的点，应该就是具体的 u 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. .pdf