代码收藏家 优化 | 稳健的优化基础
作者：苏向阳
本文讲解了鲁棒优化的基础内容，内容主要分为3个部分：1.不确定性最优化。2.不确定集。3.对等式转换理论。在末尾附上了一个简单的算例及代码便于读者理解。该文也是【运筹OR帷幄】鲁棒优化电子书系列的文章内容之一。
1. 不确定性最优化（Optimization under uncertainty）
1.1 不确定优化方法
实际生活中很多问题都具有不确定性（Uncertainty），随着最优化理论的不断发展和计算机能力的提高，不确定性优化受到了学界前所未有的重视。早在20世纪50年代，Bellman、Zadeh和Charnes等人便已开始对不确定性优化进行子研究[1,2]。在描述不确定优化问题前，我们先来看一下传统的确定性优化问题：
其中， x是决策向量，f(x)为目标函数，h(x)为约束条件函数。在模型（1）中，无论是约束条件还是目标函数，其对应的参数都是确定的。然而，在实际问题中，我们很难事先确定模型中某些参数。对于一些特定的优化问题，某个参数的细微扰动就可能导致原本所求得的最优解（Optimal solution）变得毫无意义[3]。因此，如何对不确定性条件下的问题进行优化求解就变得十分重要。
随着社会的不断发展，我们所接触到的问题的复杂度不断提高，模型的不确定性也在不断扩大。比如：飞机航班的线路规划、电网的最优调度、物流路径的最优规划等等。在实际生活中，模型参数的不确定性主要来自以下几个方面：
1）数据在统计和采集过程丢失而导致数据偏差过大；
2）天气等不可抗力因素的干扰；
3）认知不全导致现有模型与实际生活中存在偏差；
4）对于一些难以求解的非凸非线性模型，进行简化描述。
我们首先给出不确定性优化数学模型的一般表达：
在模型（2）中，ξ为不确定参数，U表示不确定参数的集合（Uncertainty set）。为求解模型（2），以Bellman等人的工作为开端，相关学者提出了一系列的优化求解方法：随机规划[4]、鲁棒优化[5]、灵敏度分析[5]、模糊规划[6]等等。不确定性优化的理论和方法不断地被开发出来。据分析阶段的不同，不确定优化的理论分为事前分析、和事后分析两大类方法。我们接下来主要对这两大类的不确定理论展开叙述。
1.2 事前分析方法
模糊规划(Fuzzy Programming)
当U是一个模糊逻辑集合时，模型（2）就成为了处理软约束（指约束条件中，等式或不等式两边含有模糊集合，因而不能像精确数学里一样判定等式成立或者不成立。）的规划问题，即模糊规划。这类规划是为了应对实际生活中无法提供准确决策的情况下而产生的一类规划求解理论。对于模糊规划问题的详细求解步骤，可以参照文献[6]。在模糊规划中，需要依据决策者的个人经验来获取不确定参数的模糊隶属函数。因此，模糊规划往往存在较大的主观性，在实际运用过程中因为需要经过反复多次调整而存在诸多限制。
随机规划(Stochastic Optimization)
当U是一个随机不确定集合时，上述模型则变成处理随机性数据的规划求解问题，即随机规划。随机规划根据不同的决策规则，可以分为三类：
1) 期望模型(Expectation model)。首先确定不确定参数的分布模型，然后通过选取离散或连续的概率分布函数对不确定参数进行描述，最终通过求取函数的期望将不确定问题转化为确定性问题并求解。如果目标函数和约束中存在随机参数，只需求取相应函数的期望值，将模型转化为确定性模型进而求解。
2) 机会约束规划模型(Chance constraint)。通俗来讲，机会约束规划是指允许决策不满足约束条件，但是决策满足约束条件的概率不低于事先设定的置信水平的规划求解模型时，目标达到最优的理论。给定置信水平，其一般化的模型描述如下：
3) 相关机会约束规划模型[7]。相关机会约束规划是当决策者面临多个事件时，希望最大化满足这些事件的概率而产生的一种规划方法。无论是期望模型还是机会约束规划模型，最终都是确定性优化求解并得出准确值。相关机会规划虽然求解结果是确定的，但并不代表一定实现，规划的目的是极大化该事件的实现概率。
随机规划的内容会在我们另一个专栏内有详细介绍，请大家多多关注。
1.3 事后分析方法
灵敏度分析（Sensitivity Analysis）是最典型的事后分析方法。灵敏度分析根据需求的不同被划分为局部灵敏度分析(Local sensitivity analysis)和全局灵敏度分析(Global sensitivity analysis)。这里以模型（4）为例对灵敏度分析展开一个简短的描述：
灵敏度分析应对的是不确定优化问题，因此有时会遇到需要添加新约束的情况。这种情况下，如果最优解满足新添加的约束，则原模型的最优解仍是新模型的最优解，需要重新计算模型。但灵敏度分析更多的研究内容是数据变化对最优解产生的影响。即：A、c和b变化导致模型最优值Z发生的改变。其研究内容是，当参数在什么范围内进行波动时，模型的最优解x∗不会发生改变，具体的原理涉及到了基变量、非基变量、对偶单纯形等相关知识，这里不再详细描述，有兴趣的读者可以参考[21]。灵敏度分析方法虽然相对其他不确定优化方法而言比较简单，但灵敏度分析方法仅是一个评价分析工具，大大限制了该方法的使用领域。
鲁棒优化（Robust Optimization）
鲁棒优化也是一类事前分析方法，之所以单独列出来，是因为鲁棒优化是针对传统优化方法的不足，由鲁棒控制理论发展而来的一套方法。在模型（2）中，如果U是一个有界闭集，上述模型则变成处理不确定集合内所有不确定参数的优化问题，即鲁棒优化。相对于传统不确定优化方法，鲁棒优化有如下优点：
1) 鲁棒优化在建模过程中充分考虑了不确定性，并以集合的形式对变量进行描述。相对于随机规划和模糊规划，鲁棒优化不需要不确定参数的分布模型和不确定参数的模糊隶属函数。
2) 鲁棒优化的约束条件是严格成立的，即只要不确定参数ξ属于不确定集合U，所求出的解都能满足约束条件。即优化模型具有较强的鲁棒性，最优解对参数变化的敏感性低。
鲁棒优化虽然有着随机规划和模糊规划没有的优势，但是鲁棒优化模型本身是一个半无限优化问题，很难直接进行求解，鲁棒优化的计算结果受限于不确定集U的不同。我们会在第2节和第3节分别对鲁棒优化中的不确定集U和鲁棒优化对等式及转换理论进行阐述。
1.4 鲁棒优化理论的发展
1973年，Soyster首次用鲁棒优化的思想来解决线性规划中的不确定性[8]。虽然该方法是基于最坏情况的基础上进行考虑，结果过于保守，但是Soyster为不确定优化的发展提供了一个全新的思路，开辟了鲁棒优化发展的道路。Mulvey等人在1995年首次提出鲁棒优化的概念[9]。他们给出了基于情景集鲁棒优化的一般模型框架，提出了鲁棒解（Solution robust）和模型鲁棒（Model robust）的概念，通过将目标函数拆分为聚合函数与罚函数来消除不确定参数对结果的影响。在此之后，不断有学者投入到鲁棒优化的研究中，在这方面的奠基之作由以色列学者Ben-Tal和Nemirovski[10,11]和美国伯克利大学的Ghaoui[12]在20世纪90年代提出。Ben-Tal证明了如果不确定集合U是一个椭球不确定集（后面具体介绍），那么对于一些最重要的一般凸优化问题(线性规划、二次约束规划、半定规划等)，其鲁棒对等式要么是精确的，要么可近似成一个可求解的问题：可以采用诸如内点法的算法在多项式时间内求解。除此之外，Ben-Tal给出了一般不确定半定规划问题的计算可处理的近似鲁棒对等式。在此之后，Ben-Tal等人又提出了可调鲁棒优化概念等概念，并被广泛运用到各行各业中。
21世纪初，Bertsimas和Sim[13]在Soyster、Ben-Tal和Nemirovski的研究基础上提出了全新的鲁棒优化框架。Bertsimas和Sim的鲁棒优化涵盖了离散优化，最主要的特点是所建立的鲁棒对等式不增加问题求解的复杂度。另一方面，Bertsimas和Sim的鲁棒优化允许出现约束违背(Constraint Violation)的情况，在这种情况下得到的鲁棒解大概率具有可行性。Bertsimas和Sim的理论由于其易处理性及实用性，受到了学界的广泛认可。
1.5 鲁棒优化问题研究路线
鲁棒优化近年来受到广泛关注，不断地被各个领域的学者应用到各行各业中，鲁棒优化问题的研究路线大致如下：
图3.1:鲁棒优化问题研究路线
前面我们已经提到，当模型（2）中的不确定集合为闭集合时，模型（2）可以视为一个鲁棒优化模型。此时，目标函数和约束条件中均含有不确定参数，为了更通俗地描述，我们对模型（2）进行转换：
转换后可以明显看出(应注意转换是否等价)，无论原鲁棒优化模型的目标函数是线性还是非线性，是否含不确定参数，都可以由模型（5）表示。但是模型（5）通常很难直接求解。为了方便求解我们需要通过数学优化理论将模型（5）转换为一个多项式时间内可求解的凸优化问题，即鲁棒对等问题(Robust Counterpart)。
目前，对相关不确定性问题用鲁棒优化方法进行研究主要体现在不确定集合的选取及鲁棒对等转换理论上：
1.不确定集合的选取。如何选取合适的不确定集合对不确定参数进行准确的描述，直接影响了模型的优化结果，而且不同的不确定集合所对应的鲁棒对等问题也不同。
2.鲁棒对等转换理论。如何把已经构建好的鲁棒优化模型转化成一个可以在多项式时间内求解的模型，直接影响了优化时间和优化结果。
2. 不确定集（Uncertainty set)
鲁棒优化中，不同的不确定集对结果影响十分明显，当不确定集越精细时，模型复杂度越高，求解越困难。当不确定集越宽泛时，所求出的最优解越保守，越不符合实际。为了权衡二者的关系，如何选择一个适合的不确定集来求解鲁棒优化一直是相关学者的研究热点。常见的不确定集主要有如下几类:
1) 盒式不确定集（Box Uncertainty Set）
盒式不确定集是最简单的不确定集，也被称作区间集。由于鲁棒优化是考虑最差情况下的优化求解方法，针对于一些模型可能会出现所有不确定参数都在区间集上下界进行优化的情况，然而实际中该情况发生的概率可能极低或不会发生，因此，很容易出现过度保守的情况。
2) 椭球不确定集（Ellipsoidal Uncertainty Set）[5,14]：椭球集/椭球交集
在Ben-tal的经典著作《Robust Optimization》称集合(7)为椭球不确定集，集合(8)为椭球交集不确定集。在Boyd的经典著作《Convex Optimization》中称(8)为椭球集，(7)为退化的椭球。为了方便描述，本文以Ben-tal的描述为准。上述公式中，ξ为不确定参数向量，U为不确定参数的期望或预测值向量，R为协方差矩阵，Ω为不确定度，用以刻画不确定参数扰动范围。相对于椭球集，椭球交集能更准确地对不确定参数进行描述，但是求解二次规划（quadratic programming）、锥二次优化（Second order cone programming）、半定规划（Semidefinite programming）问题时，如果不确定集是椭球交集，模型将难以直接求解。Ben-tal已经证明了这些优化问题中使用椭球交集时是NP-hard问题[Citation]。如果采用椭球集，在线性规划、二次优化问题和锥二次优化时可以转化为可处理问题，但是在半定规划中，仍需满足诸多限制才能求解。椭球不确定集虽然可以很好地表示很多类型集合，方便数据输入，在一定程度上可以体现不确定参数之间的关联性。但是椭球不确定集会增加问题求解的复杂度，因此未能得到广泛运用。
3) 多面体不确定集（Polyhedral Uncertainty Set）
多面体集合[15]可以看作是椭球集合的一种特殊表现形式[11]。尽管多面体不确定集难以刻画不确定参数间的相关性，但其具有线性结构、易于控制不确定度，在实际工程问题中广受青睐。
4) 基数/预算不确定集（Budget uncertainty Set）
5) 数据驱动不确定集（Data-driven Uncertainty Set）:无论是采用盒式还是椭球式不确定集合，都会出现所得到的解过于保守的情况。为了解决解的过度保守，一些学者根据历史数据进行不确定集的构造，也被称为数据驱动不确定集。
数据驱动不确定集的构建，是使用统计假设检验的置信区间来精确描述不确定参数的分布。在04年Bertsimas、Sim和09年Ben-Tal等人的研究中，假定不确定参数为，其分布不能精确获得。最初始的研究是基于数据结构特性做出的先验假设。这些方法假设是没有依赖的，但是不会认为的边界分布是确定的。在13年Bertsimas对数据驱动不确定集的构造进行了进一步的改进。Bertsimas假设数据S有独立分布，这些S可以为不确定参数的分布添加更多的细节。并且通过这些细节信息，设计一个概率保证的集合，相对于传统的集合，新的集合更小，所求出的结果也不是那么的精确。关于数据驱动的先验假设和假设检验可以参照Bertsimas的研究[19]。
除去以上常见的不确定集，一些学者为了适应不同的情况以及更精确地对不确定参数进行描述，还衍生出了很多种组合不确定集，具体如：盒式+椭球式不确定集、盒式+多面体不确定集、盒式+椭球式+多面体不确定集等等。
3. 对等式转换理论（Robust counterpart）
前文我们已经提到，鲁棒优化是一个半无限优化问题，通常情况下很难直接求解。为了对鲁棒优化问题进行求解，我们需要对原模型做出一定程度的转化，转化后模型需要在多项式时间内可以求解。在这方面做出重要突破的学者有Soyster、Ben-tal、Bertsimas和Sim等。Soyster提出用线性优化模型求解问题，并使用凸集内所有数据，但该方法存在过于保守的弊端。Ben-tal和Nemirovski、Bertsimas和Sim在Soyster的基础上进行了改进，使得鲁棒优化适用性更广。接下来我们分别对这些理论进行描述。
3.1 Soyster的鲁棒对等模型
3.2 Ben-tal和Nemirovski鲁棒对等模型
因为随机变量对称分布，所以Ben-tal和Nemirovski采用随机变量的期望和方差来进行替换相关约束问题。（14）中Ω大于0，需要事先给定。我们可以明显看出，模型（12）的可行解是模型（14）可行解的子集。因此，Ben-tal和Nemirovski的方法保守度更低一定，同时可以通过对Ω的调整来改变问题的保守度。但是模型（14）不再是线性规划，模型求解困难度上升，这也是Ben-tal和Nemirovski鲁棒优化的不足之处。
Ben-tal和Nemirovski的方法虽然能降低模型的保守程度，但是并不能求解离散优化问题，针对这一痛点，Bertsimas和Sim提出了一种全新的鲁棒架构[13]
3.3 Bertsimas和Sim鲁棒对等模型
模型转化的理论依据是对偶理论，具体的证明过程在作者稿件中有详细描述，有兴趣的读者可以自行翻阅，这里不再赘述。
在此基础上，Bertsimas和Sim提出了鲁棒离散优化的模型[20]。
模型（11）的鲁棒对等模型如下所示：
4. 算例
为了方便大家理解，附上了一个小例子进行说明。我们以(13)中的仿真算例进行说明：
该问题在Python3.6环境下调用gurobi的求解代码如下
(https://github.com/Feeling-well/robust-optimization)：
"""
learn from xprog and bertsimas's paper(price of robustness)
http://xprog.weebly.com/
the initial model cann't be solved by gurobi or cplex
the problem has an equivalent linear formulation as follows:
max z
s.t. z<=sum(p[i]x[i])-(sum(Q[i])+Γ[i]m[i])
Q[i]+m[i]>=σ[i]x[i]
Q[i]>=0
m[i]>=0
sum(x[i])=1
x[i]>=0
Γis the protection level of the actual portfolio return in the following sense
"""
from gurobipy import *
try:
m = Model("RO")
#establishment of constant
σ = []
p = []
Γ= []
for n in range(1,151):
σ.append(0.05/450*(2*n*150*151)**0.5)
p.append(1.15+n*0.05/150)
Γ.append(5)
#Add variables
x = m.addVars(150,lb=0,name='x')
z = m.addVar(name='z')
Q = m.addVars(150,name='Q')
mm = m.addVars(150,name='m')
px = sum(p[i]*x[i] for i in range(150))
QC =sum(Q[i] for i in range(150))
# Add Constrs
m.addConstrs((z<=px-mm[i]*Γ[i]-QC for i in range(150)),name='first')
m.addConstrs((mm[i]+Q[i]>=σ[i]*x[i] for i in range(150)),name='second')
m.addConstr(sum(x[i] for i in range(150))==1,name='third')
#set obj
m.setObjective(z,GRB.MAXIMIZE)
#print model
m.write ("RO.lp")
#solve model
m.optimize()
#print variables
# for v in m.getVars ():
# print ('%s %g' % (v.varName, v.x))
except GurobiError as e:
print('Error code ' + str(e.errno) + ": " + str(e))
except AttributeError:
print('Encountered an attribute error')参考文献：
[1].Charnes,A., & Cooper, W. W., Chance-constrained programming. Managementscience, 1959. 6(1), p. 73-79.
[2].Bellman,R. E., & Zadeh, L. A., Decision-making in a fuzzy environment. Managementscience, 1970. 17(4), p. B-141.
[3].ElGhaoui, L., Oustry, F., & Lebret, H., Robust solutions to uncertainsemidefinite programs. SIAM Journal on Optimization, 1998. 9(1),p. 33-52.
[4].Birge,J. R., & Louveaux, F., Introduction to stochastic programming.Springer Science & Business Media. 2011.
[5].Ben-Tal,A., El Ghaoui, L., & Nemirovski, A., Robust optimization (Vol.28). Princeton University Press. 2009.
[6].刘宝碇, & 赵瑞清. (1998). 随机规划与模糊规划. 清华大学出版社.
[7]. Liu, B., Dependent-chance programming: A class of stochastic optimization. Computers & Mathematics with Applications, 1997. 34(12): p. 89-104.
[8]. Soyster, A.L., Technical Note – Convex Programming with Set-Inclusive Constraints and Applications to Inexact Linear Programming. Operations Research, 1973. 21(5): p. 1154-1157.
[9]. Mulvey,J.M., Vanderbei, R.J. and Zenios, S.A.
[10]. Ben-Tal, A. and A. Nemirovski, Robust Convex Optimization. Mathematics of Operations Research, 1998. 23(4): p. 769-805.
[11]. Ben-Tal, A. and A. Nemirovski, Robust solutions of uncertain linear programs. Operations Research Letters, 1999: p. 1-13.
[12].ElGhaoui L, Lebret H. Robust solutions to least-squares problems with uncertaindata. SIAM Journal on matrix analysis and applications, 1997.18(4): p.1035-64.
[13]. Bertsimas, D. and M. Sim, The Price of Robustness. Operations Research, 2004. 52(1): p. 35-53.
[14].BoydS, Vandenberghe L. Convex optimization. Cambridge university press. 2004.
[15].BertsimasD, Thiele A. Robust and data-driven optimization: modern decision making underuncertainty. InModels, methods, and applications for innovative decisionmaking, 2006. p. 95-122. INFORMS.
[16]. Minoux, M., On robust maximum flow with polyhedral uncertainty sets. Optimization Letters, 2009. 3(3): p. 367-376.
[17].JiangR, Zhang M, Li G, Guan Y. Two-stage robust power grid optimization problem. OptimizationOnline. 2010.
[18]. Baringo, L. and A. Baringo, A Stochastic Adaptive Robust Optimization Approach for the Generation and Transmission Expansion Planning. IEEE Transactions on Power Systems, 2017. PP(99): p. 1-1.
[19]. Bertsimas, D., V. Gupta and N. Kallus, Data-driven robust optimization. Mathematical Programming, 2018. 167(2): p. 235-292.
[20]. Bertsimas, D. and M. Sim, Robust discrete optimization and network flows. Mathematical Programming, 2003. 98(1-3): p. 49-71.
[21] 陈宝林, 最优化理论与算法, 清华大学出版社, 2005.
