本文的主要参考文献:
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表示第二阶段决策变量,表示不确定矢量。在此假设下的两阶段鲁棒优化的一般形式为:
![1d6363b719fa4d598a9234f6874282f6.png](https://img-
blog.csdnimg.cn/1d6363b719fa4d598a9234f6874282f6.png)
其中:
![ce39e4319e764fceb30fbe6c4bf53de4.png](https://img-
blog.csdnimg.cn/ce39e4319e764fceb30fbe6c4bf53de4.png)
向量c,b,d,h和矩阵A , G , E ,
M都是确定性的数值,不确定性体现在向量u上。注意到第二阶段优化的约束条件F(y,u)是关于不确定性u的线性函数。
原文献中提供了以运输问题作为算例,具体如下:
![8c2885d2e5c44711b39dd8de41bc4e89.png](https://img-
blog.csdnimg.cn/8c2885d2e5c44711b39dd8de41bc4e89.png)
其中,yi为0-1变量,表示是否在i地建设仓库,zi表示仓库i储存的商品数量,xij表示从i仓库到j客户运送的商品数量,fi表示建设仓库i的固定成本,ai表示仓库i存储商品的单位成本,cij表示从i仓库到j客户运送单位商品的成本,ki表示仓库i的最大容量,dj表示客户j的需求。
不确定变量为客户的需求,表达方式如下:
![8fcfa7006e894e0bb1dfda86dfbce241.png](https://img-
blog.csdnimg.cn/8fcfa7006e894e0bb1dfda86dfbce241.png)
具体算例参数:
![afef57b843014c028500faf2ab70ffbf.png](https://img-
blog.csdnimg.cn/afef57b843014c028500faf2ab70ffbf.png)
根据上面的公式,我们可以写出各个参数矩阵以及变量的表达式:
![4a02d9e458c44f31bb1e5cbc2bf56e2e.png](https://img-
blog.csdnimg.cn/4a02d9e458c44f31bb1e5cbc2bf56e2e.png)
用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) <= z'];
Constraints=[Constraints ,sum(x,2) >= d];
Constraints=[Constraints ,d == d_ + g*d_wave];
Constraints=[Constraints ,g'*P <= gamma];
%% 设置求解器
ops=sdpsettings('verbose', 3, 'solver', 'gurobi');
sol=optimize(Constraints,objective,ops);
优化结果为:
![4bbd9049d5eb4fc5899642facbd6a8bb.png](https://img-
blog.csdnimg.cn/4bbd9049d5eb4fc5899642facbd6a8bb.png)
进一步把算例写成两阶段鲁棒优化的形式:
![2756e85f1d564b598e5d53c6a169e595.png](https://img-
blog.csdnimg.cn/2756e85f1d564b598e5d53c6a169e595.png)
针对这个两阶段鲁棒优化问题,可以分别采用Benders对偶割平面法和C&CG算法进行求解。
2.Benders对偶割平面法
2.1基本原理
Benders对偶割平面法可以用于解决两阶段鲁棒优化问题,首先将两阶段鲁棒优化问题分解为两部分:主问题(Master
Problem,MP)和子问题(Subproblem,SP)。主问题包含第一阶段的决策变量y以及仅与y有关的约束和子问题返回的割,还包括辅助变量η,用于评估第二阶段目标函数的取值。子问题包含第二阶段的决策变量x和不确定变量u,旨在给出第二阶段目标函数值的一个界限值。针对式(1)中描述的两阶段鲁棒优化问题,其主问题
** MP** 可以写成:
![29a0e85154e448f7b5c3ec915624b018.png](https://img-
blog.csdnimg.cn/29a0e85154e448f7b5c3ec915624b018.png)
主问题是一个线性规划问题。子问题 ** SP** 则为:
![dd7855c30cfe4ea8bde869e002ffd964.png](https://img-
blog.csdnimg.cn/dd7855c30cfe4ea8bde869e002ffd964.png)
而子问题是一个双层线性规划问题(如果不知道双层规划的概念,可以去看看我之前的几篇博客双层优化入门-
CSDN博客),其中上层优化的决策变量是u,下层优化的决策变量是x,而且在下层优化中,变量y和u的值都是确定的,可以视为参数。
对于双层优化形式的子问题的求解,主要有以下几种方式:
1.通过对偶变换将双层优化问题转为单层优化问题,再进行求解,可以使用智能优化算法、等价线性化、二次规划求解器(例如gurobi)等方式进行求解;
2.采用智能优化算法进行求解(可参考博客双层优化入门(3)—基于智能优化算法的求解方法);
3.采用KKT条件进行求解(可参考博客双层优化基本原理与求解方法、基于yalmip的双层优化求解)。
在这里我们使用KKT条件来求解子问题,可以将双层优化的子问题转换为下列单层优化的形式:
![4a4347ddf1f24b3791157143092356ae.png](https://img-
blog.csdnimg.cn/4a4347ddf1f24b3791157143092356ae.png)
(补充说明)
上面给出的式子是经过化简的,并不是最初的KKT条件,因此有读者和我反映,看起来有点懵圈。我在此把详细的转换步骤写一下:
首先由于子问题的内层优化问题是min形式,需要把约束条件都转为≥0的形式:
根据这个优化问题的形式,可以写出子问题内层优化的拉格朗日函数:
其中,π和θ都是拉格朗日乘子。根据拉格朗日函数进一步写出KKT条件:
这就是初始的KKT条件,接下来看一下文献中是如何对这组式子进行化简的。首先根据式(1)和(4)可以得到:
然后,根据式(1)和(3)可以得到:
这样化简之后,就得到了和原文中相同的形式。根据化简的结果,大家应该能知道化简的目的: ** 消去变量** ** θ** **
,减少优化问题中变量的数目** 。实际上几乎所有的KKT条件都可以这样进行处理,之后有机会再和大家讲讲,这里不再赘述。
化简后的KKT条件也存在非线性形式,可以使用大M法引入二进制中间变量进行线性化,将其转换为混合整数规划的形式:
对于一般形式的两阶段鲁棒优化问题(如式(1)),Benders对偶切平面算法求解的流程如下:
步骤1:设定目标函数上界UB=+∞, 下界LB=-∞,迭代次数k=0。
步骤2:求解主问题 ** MP** :
![6ed4d2650ac64246919c6a29e40a8ed0.png](https://img-
blog.csdnimg.cn/6ed4d2650ac64246919c6a29e40a8ed0.png)
求出最优解(yk+1*,ηk+1*),并更新LB=max{LB,c T yk+1*+ηk+1*};
步骤3:求解子问题 ** SP** :
![ac4e316b0d334db3a2b1ff0a7e844c16.png](https://img-
blog.csdnimg.cn/ac4e316b0d334db3a2b1ff0a7e844c16.png)
求出子问题的最优目标函数值Qk+1以及最优解( ** u* k+1*, ** x** k+1*),并更新UB=min{UB,cT yk
+Q(yk)}。
步骤4:如果UB-LB ≤
ε(事先设定的运行偏差),则输出优化结果,并退出循环。否则令k=k+1,将约束![6c73f755451947299a3a949fcdfc3275.png](https://img-
blog.csdnimg.cn/6c73f755451947299a3a949fcdfc3275.png)添加到主问题 ** MP**
中并返回步骤2。
从上述步骤中可以看到,算法迭代的过程会不断向主问题添加约束条件![6c73f755451947299a3a949fcdfc3275.png](https://img-
blog.csdnimg.cn/6c73f755451947299a3a949fcdfc3275.png),也就是割平面,因此被称为Benders对偶割平面法。
2.2 算例分析
采用文献中给出的运输问题作为算例,使用matlab+yalmip工具箱+gurobi求解器进行求解。
![151094cb6ab74e8c8d70da00431dd425.png](https://img-
blog.csdnimg.cn/151094cb6ab74e8c8d70da00431dd425.png)
为了求解这个两阶段鲁棒优化问题,我们首先需要把这个优化问题分解成主问题和子问题。而且为了方便理解,重写成符合标准两阶段鲁棒优化问题的形式,其中重写后的优化问题部分变量或系数矩阵和原优化问题中重复,我都加了上标一撇(‘)以示区别,具体步骤如下:
主问题 ** MP** ** _BD** :
![b3d1cc82fa974cfca8c22c92a19a26b1.png](https://img-
blog.csdnimg.cn/b3d1cc82fa974cfca8c22c92a19a26b1.png)
子问题 ** SP** ** _BD** :
![8c5f695c3bec4048824052dc71f3f32a.png](https://img-
blog.csdnimg.cn/8c5f695c3bec4048824052dc71f3f32a.png)
步骤1:设定目标函数上界UB=+∞, 下界LB=-∞,迭代次数k=0。
步骤2:求解主问题 ** MP** ** _BD** ,得到最优解(,),并更新LB=max{LB,};
步骤3:求解子问题 ** SP** ** _BD** ,得到子问题的最优目标函数值Qk以及最优解( ** u* k*, **
x** k*),并更新UB=min{UB,![83ad75789a1347dcaedd5c30c9b5a663.png](https://img-
blog.csdnimg.cn/83ad75789a1347dcaedd5c30c9b5a663.png)}。
步骤4:如果UB-LB ≤
ε(事先设定的允许偏差),则输出优化结果,并退出循环。否则令k=k+1,将约束![8edacc06710b42d7b258900663e62db6.png](https://img-
blog.csdnimg.cn/8edacc06710b42d7b258900663e62db6.png)添加到主问题 ** MP**
中并返回步骤2。
采用matlab编程进行求解,结果如下:
![8db76cced9354b058ae5d7c0cb5922c2.png](https://img-
blog.csdnimg.cn/8db76cced9354b058ae5d7c0cb5922c2.png)![86e83352b2c643039fc8d38b487a4f6b.png](https://img-
blog.csdnimg.cn/86e83352b2c643039fc8d38b487a4f6b.png)![253fcd100d2d47b5a738f7641dd44338.png](https://img-
blog.csdnimg.cn/253fcd100d2d47b5a738f7641dd44338.png)
与确定性优化的结果对比如下:
变量
|
确定性优化
|
两阶段鲁棒优化
—|—|—
最优目标函数值
|
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** :
![71e7eca16ded453780a6de82eb0fa8fd.png](https://img-
blog.csdnimg.cn/71e7eca16ded453780a6de82eb0fa8fd.png)
其中,x _ l_ 是第 _ l_
次迭代添加的决策变量,![2b90752237e44087ab16455f15390217.png](https://img-
blog.csdnimg.cn/2b90752237e44087ab16455f15390217.png)是第 _ l_
次迭代子问题的最恶劣场景,因为第 _ l_
次迭代时子问题中已经求解了![2b90752237e44087ab16455f15390217.png](https://img-
blog.csdnimg.cn/2b90752237e44087ab16455f15390217.png)的值,所以在主问题中就可以看作已知的参数。
子问题 ** SP** ** _CCG** :
![9e67395329f14aba97be42e47dc56c49.png](https://img-
blog.csdnimg.cn/9e67395329f14aba97be42e47dc56c49.png)
步骤1:设定目标函数上界UB=+∞, 下界LB=-∞,迭代次数k=0。
步骤2:求解主问题 ** MP** ** _CCG**
,得到最优解![9021aef67e764e3f9929f2cf27b53bb3.png](https://img-
blog.csdnimg.cn/9021aef67e764e3f9929f2cf27b53bb3.png),并更新LB=max{LB,![184820b7ccaf4af1971e0f8f2f5d96e8.png](https://img-
blog.csdnimg.cn/184820b7ccaf4af1971e0f8f2f5d96e8.png)};
步骤3:求解子问题 ** SP** ** _CCG** ,得到子问题的最优目标函数值Qk以及最优解( ** u* k*, **
x** k*),并更新UB=min{UB,![59b8c175ae6444a888f634e1c920f897.png](https://img-
blog.csdnimg.cn/59b8c175ae6444a888f634e1c920f897.png)}。
步骤4:如果UB-LB ≤ ε(事先设定的允许偏差),则输出优化结果,并退出循环。否则转到步骤5.
步骤5:判断子问题是否存在最优解。
若 _ Q_ (yk+1*)<+∞,即子问题存在最优解,则创建变量xk+1并给主问题 ** MP** ** _CCG** 添加以下约束:
![eb5cb031bc444eb3af389fc2ca64d88a.png](https://img-
blog.csdnimg.cn/eb5cb031bc444eb3af389fc2ca64d88a.png)
其中uk+1*是第k次迭代时子问题求解出来的最恶劣场景。
随后更新k= k+1,并转到步骤2。
若 _ Q_ (yk+1*)=+∞,即子问题不存在最优解,则创建变量xk+1并给主问题 ** MP** ** _CCG**
添加以下约束:
![fbce73f6355246c99dfc6312e01dab34.png](https://img-
blog.csdnimg.cn/fbce73f6355246c99dfc6312e01dab34.png)
随后更新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是不确定集
** U** 的极值点数,q为满足GTπ≤b和π≥0的集合{π}中的极值点数量,具体见原文献中的描述)。因此C&CG算法的收敛速度要更快。
(4) 在求解问题的能力方面,Benders-
dual对偶割平面算法需要将第二阶段的问题转换为线性规划问题,而C&CG算法不关心第二阶段的变量类型。如果第二阶段是混合整数规划,可以采用嵌套C&CG算法进行求解。
3.2 算例分析
同样采用文献中给出的运输问题作为算例,使用matlab+yalmip工具箱+gurobi求解器进行求解。
![5a88e29763804dd3a671fbf922bd989a.png](https://img-
blog.csdnimg.cn/5a88e29763804dd3a671fbf922bd989a.png)
和Benders-
dual对偶割平面算法一样,我们首先需要把这个优化问题分解成主问题和子问题,并将优化问题重写成符合标准两阶段鲁棒优化问题的形式,具体步骤如下:
主问题 ** MP** ** _CCG** :
![a592c70729fb41f2bf9835871f3d1aed.png](https://img-
blog.csdnimg.cn/a592c70729fb41f2bf9835871f3d1aed.png)
子问题 ** SP** ** _BD** :
![6abb1665a9634cee85b2198458eab7b6.png](https://img-
blog.csdnimg.cn/6abb1665a9634cee85b2198458eab7b6.png)
步骤1:设定目标函数上界UB=+∞, 下界LB=-∞,迭代次数k=0。
步骤2:求解主问题 ** MP** ** _CCG**
,得到最优解![c0006214ed1740bbbf3ba5ec49eb50dc.png](https://img-
blog.csdnimg.cn/c0006214ed1740bbbf3ba5ec49eb50dc.png),并更新LB=max{LB,![c83559d883fa4a708b2d3d05a1410db0.png](https://img-
blog.csdnimg.cn/c83559d883fa4a708b2d3d05a1410db0.png)};
步骤3:求解子问题 ** SP** ** _CCG** ,得到子问题的最优目标函数值Qk以及最优解( ** u* k*, **
x** k*),并更新UB=min{UB,![8e9bdee2a277423d9e2aef636b5a79e7.png](https://img-
blog.csdnimg.cn/8e9bdee2a277423d9e2aef636b5a79e7.png)}。
步骤4:如果UB-LB ≤ ε(事先设定的允许偏差),则输出优化结果,并退出循环。否则转到步骤5.
步骤5:判断子问题是否存在最优解。
若 _ Q_ (yk+1*)<+∞,即子问题存在最优解,则创建变量xk+1并给主问题 ** MP** ** _CCG** 添加以下约束:
![bcdc79cde0ca43a6bbcc39cf8f3797ce.png](https://img-
blog.csdnimg.cn/bcdc79cde0ca43a6bbcc39cf8f3797ce.png)
其中uk+1*是第k次迭代时子问题求解出来的最恶劣场景。随后更新k= k+1,并转到步骤2。
若 _ Q_ (yk+1*)=+∞,即子问题不存在最优解,则创建变量xk+1并给主问题 ** MP** ** _CCG**
添加以下约束:
![2df6c75c072b4385a10435127ff73d33.png](https://img-
blog.csdnimg.cn/2df6c75c072b4385a10435127ff73d33.png)
随后更新k= k+1,并转到步骤2。
运行结果如下,和Benders对偶割平面方法一样,但收敛速度更快:
![dc6bf08586ba42a8a85d56e252b888ce.png](https://img-
blog.csdnimg.cn/dc6bf08586ba42a8a85d56e252b888ce.png)
![70b07236984e44a1bc0779f2f2ebfc86.png](https://img-
blog.csdnimg.cn/70b07236984e44a1bc0779f2f2ebfc86.png)
![52322b88b07e4c638abc15efc96f59e3.png](https://img-
blog.csdnimg.cn/52322b88b07e4c638abc15efc96f59e3.png) 4 ** .
** 完整代码获取链接****
想获取完整代码,可以戳下面这个链接:
两阶段鲁棒优化以及列与约束生成算法(C&CG)的matlab代码
网络安全工程师(白帽子)企业级学习路线
第一阶段:安全基础(入门)
第二阶段:Web渗透(初级网安工程师)
第三阶段:进阶部分(中级网络安全工程师)
如果你对网络安全入门感兴趣,那么你需要的话可以点击这里👉网络安全重磅福利:入门&进阶全套282G学习资源包免费分享!
学习资源分享
![](https://i-blog.csdnimg.cn/blog_migrate/451e72f06f23da2cd3dabe0fa84e00b9.jpeg)