随机优化中的样本均值近似方法

补一下之前随机优化的学习内容，结合Gurobi的webinar中随机优化的讲座和gurobi求解器进行求解。

1.SAA原理

随机优化是求解包含随机变量的一类优化问题，其中随机变量的分布是确定已知的。而样本均值近似（sample average approximation,SAA）方法是求解随机优化的常用方法，其原理部分主要参考翻译文献[1]，其中关于SAA的收敛速率的证明，有兴趣的同学可以读读原文献，数学太难，在此省略。

考虑以下包含期望的随机优化问题：
$$\underset{x\in \mathbf{S}}{\min }g(x)={\mathbb{E}}_{W\sim \mathbf{P}}G(x,W)\text{\hspace{1em}(1)}$$

其中，$W$为随机变量，其概率分布服从$W\sim \mathbf{P}$。优化变量属于有限集合（finite set）$x\in \mathbf{S}$。$w$是随机变量$W$的样本，$G(x,w)$是$G(x,W)$在样本$w$下的一次实现。
式（1）中的期望${\mathbb{E}}_{W\sim \mathbf{P}}G(x,W)=\int G(x,W)f(W)dW$是关于随机变量$W$的积分,$f(W)$是$W$的概率密度函数。

使用SAA的前提[1]
使用SAA需要满足以下三个前提假设：

• （1） 期望函数$g(x)$没有闭合表达式，并且其数值不容易计算
• （2） 函数$G(x,w)$对于给定的$x$和$w$是容易计算的
• （3） 可行解集虽是有限的，但也是数量足够大的，不能使用穷举方法。

定义函数$g(x)$的样本均值近似为：
$${\hat{g}}_N(x)=\frac{1}{N}\sum _{j=1}^{N}G(x,w_j)$$
其中，$w_j$为$W$的第$j$次采样值。

相应的优化问题转化为：
$$\underset{x\in \mathbf{S}}{\min }{\hat{g}}_N(x)\text{\hspace{1em}(2)}$$
把优化问题（1）称为真问题（原始问题），优化问题（2）称为样本均值近似问题，需要注意到${\hat{g}}_N(x)$是$g(x)$的无偏估计量，即$\mathbb{E}[{\hat{g}}_N(x)]=g(x)$ 。

SAA两个重要性质[2]：

• 渐进收敛性：随着样本数量N趋于无穷大，问题（2）的最优解和最优值和收敛于原始问题（1）的最优解和最优值。
• 易处理性：对于大多数函数$G(x,W)$和可行解集$\mathbf{S}$,找到优化问题（2）的最优解和最优值，在计算上是易处理的。

2.使用Gurobi实操求解

在此使用报童问题(the newsvender problem)作为示例，主要参考gurobi在2018年的webinar《Solving Simple Stochastic Optimization Problems with Gurobi》[2].
问题建模如下:
max ⁡ E [ x p ] s . t .    0 ≤ x s ≤ min ⁡ { y , d } 0 ≤ x d ≤ max ⁡ { 0 , y − x s } x p = c s x s − c o y + c d x d (3) \begin{aligned} & \max{\Bbb{E}[{x_{p}}]}\\ & s.t.\; 0 \leq x_{s}\leq \min{\{ y,d\}} \\ & \quad 0 \leq x_{d}\leq \max{\{ 0,y-x_{s}\}} \\ & \quad x_{p}=c_{s}x_{s}-c_{o}y+c_{d}x_{d} \end{aligned}\tag{3} ​maxE[xp​]s.t.0≤xs​≤min{y,d}0≤xd​≤max{0,y−xs​}xp​=cs​xs​−co​y+cd​xd​​(3)

其中,$x_p$ 为总体销售利润， $x_s$为销量， $c_s$为售价， $y$为进货量， $c_o$为成本， $x_d$为打折销售的数量（滞销的数量）， $c_d$为滞销的售价，如果$c_d\ge 0$表示滞销仍有收益，相反$c_d<0$表示赔本销售。

需要注意的是$d$表示市场的需求量，为随机变量，具有不确定性。而进货量$y$为最终要求解的决策变量。

使用SAA方法求解：

max ⁡ ∑ i = 1 N x p , i N , ∀ i ∈ { 1 , ⋯   , N } s . t .    0 ≤ x s , i ≤ min ⁡ { y , d i } 0 ≤ x d , i ≤ max ⁡ { 0 , y − x s , i } x p , i = c s x s , i − c o y + c d x d , i (4) \begin{aligned} &\max \sum_{i=1}^{N}\frac{x_{p,i}}{N}, \forall i \in \{1,\cdots,N\} \\ & s.t.\; 0 \leq x_{s,i}\leq \min{\{ y,d_{i}\}} \\ & \quad 0 \leq x_{d,i}\leq \max{\{ 0,y-x_{s,i}\}} \\ & \quad x_{p,i}=c_{s}x_{s,i}-c_{o}y+c_{d}x_{d,i} \end{aligned}\tag{4} ​maxi=1∑N​Nxp,i​​,∀i∈{1,⋯,N}s.t.0≤xs,i​≤min{y,di​}0≤xd,i​≤max{0,y−xs,i​}xp,i​=cs​xs,i​−co​y+cd​xd,i​​(4)

$d_i$为随机变量第i次的抽样值，$x_{p,i},x_{s,i},x_{d,i}$为在第i次抽样对应的优化变量。而最终的求解的变量y在每次抽样中都保持不变。此时的问题（4）转化为包含$3\ast N+1$个变量的确定性优化问题。

2.1 问题求解

进一步，优化问题（4）可以转化为：
max ⁡ ∑ i = 1 N x p , i N , ∀ i ∈ { 1 , ⋯   , N } s . t .    0 ≤ x s , i ≤ d i x s , i + x d , i = y x p , i = c s x s , i − c o y + c d x d , i x d , i ≥ 0 (5) \begin{aligned} & \max \sum_{i=1}^{N}\frac{x_{p,i}}{N}, \forall i \in \{1,\cdots,N \}\\ & s.t.\; 0 \leq x_{s,i}\leq d_{i} \\ & \quad x_{s,i}+x_{d,i}=y \\ & \quad x_{p,i}=c_{s}x_{s,i}-c_{o}y+c_{d}x_{d,i} \\ & \quad x_{d,i}\geq 0 \end{aligned}\tag{5} ​maxi=1∑N​Nxp,i​​,∀i∈{1,⋯,N}s.t.0≤xs,i​≤di​xs,i​+xd,i​=yxp,i​=cs​xs,i​−co​y+cd​xd,i​xd,i​≥0​(5)

使用gurobi求解优化问题(5)，代码如下[3]：

#导包
from gurobipy import *
import random
random.seed(a=100) #设置随机生成器的种子，保证可重复性
import matplotlib.pyplot as plt

#参数设置
cost    = 2  #成本
retail  = 15 #售价
recover = -3 #滞销的价格
samples = 10000 #样本数量
# 对于市场需求量，使用截断正态分布模拟
sigma   = 100
mu      = 400 #均值
demand  = [max(random.normalvariate(mu,sigma),0) for i in range(samples)]

# 总体利润 x_{p}的最大最小值，约束x_{p}
maxrev  = max(demand)*(retail-cost) #最大值（无滞销）
minrev  = max(demand)*(recover-cost)+min(demand)*retail #最小值

#建模求解
m = Model()
# Set to maximize
m.ModelSense = -1
# 添加变量 
#gurobi中添加的变量默认情况下是 >=0 的。
order    = m.addVar(name='order') #对于优化变量y
profit   = m.addVars(samples,obj=1.0/samples,lb=minrev,ub=maxrev,name='profit') #obj为目标函数中变量的系数，此处为1.0/samples
sales    = m.addVars(samples,ub=demand,name='sales') #ub是变量的上界，lb是下界
discount = m.addVars(samples,name='discount')
# 设置约束条件
m.addConstrs((profit[i] == sales[i] * retail - order * cost + recover * discount[i] for i in range(samples)),name='profit')
m.addConstrs((sales[i]+discount[i] == order for i in range(samples)),name='demand')
m.update()
m.optimize() #求解
print("进货量:",order.x)
print("总体利润:",m.objVal)

结果如下：

考察样本数量$N$变化对SAA方法的影响，如下图所示。设置最大的样本数量为10000个，每次迭代增加500个样本，可见随着样本数量增加，优化变量$y$和目标函数都在逐步下降，趋于收敛。

2.2 最坏的情况（worst case）

考察一下，最坏的情况。假设用户希望即便是在利润最差的情况下也能有较好的结果，这既是鲁棒优化问题。此时的目标函数为：
max ⁡ y { min ⁡ x p } \max_{y} \{\min x_{p}\} ymax​{minxp​}
进一步，建立的优化问题为：
max ⁡ y ω s . t .    ω ≤ x p , i , ∀ i ∈ { 1 , ⋯   , N } 0 ≤ x s , i ≤ d i x s , i + x d , i = y x p , i = c s x s , i − c o y + c d x d , i x d , i ≥ 0 (6) \begin{aligned} & \max_{y} \omega \\ & s.t. \; \omega \leq x_{p,i},\forall i \in \{1,\cdots,N \} \\ &\quad 0 \leq x_{s,i}\leq d_{i} \\ & \quad x_{s,i}+x_{d,i}=y \\ & \quad x_{p,i}=c_{s}x_{s,i}-c_{o}y+c_{d}x_{d,i} \\ & \quad x_{d,i}\geq 0 \end{aligned} \tag{6} ​ymax​ωs.t.ω≤xp,i​,∀i∈{1,⋯,N}0≤xs,i​≤di​xs,i​+xd,i​=yxp,i​=cs​xs,i​−co​y+cd​xd,i​xd,i​≥0​(6)
其中，第一个约束条件保证$\omega$为$x_{p,i}$的最小值。
求解代码如下：

#主要在模型建立方面发生变化，其他参数与2.1节保持一致
m = Model()
# Set to maximize
m.ModelSense = -1
# Add variables
worst    = m.addVar(lb=minrev,ub=maxrev,obj=1,name='worst') #obj为变量在目标函数中的系数
order    = m.addVar(name='order')
profit   = m.addVars(samples,lb=minrev,ub=maxrev,name='profit')
sales    = m.addVars(samples,ub=demand,name='sales')
discount = m.addVars(samples,name='discount')
# Set constraints
m.addConstrs((profit[i] == -order * cost + sales[i] * retail + recover * discount[i] for i in range(samples)),name='profit')
m.addConstrs((sales[i]+discount[i] == order for i in range(samples)),name='demand')
m.addConstrs((worst <= profit[i] for i in range(samples)),name='worst')# worst <= min(profit[i])
m.update()
m.optimize()
print("worst-case进货量:",order.x)
print("worst-case总体利润:",m.objVal)

结果如下图，此时的目标函数为504，性能较差，可见worst-case情况下求出的最优解比较保守。

2.3 机会约束优化（chance constrained optimization）

(1) 25% worst-cast

在一般的鲁棒优化worst-case情况下，求得的最优解和最优值结果太悲观，实际的情况往往好于worst-case情况。因此可以使用机会约束（chance constrained）去松弛（中和）worst-case情况下的强约束（也称不可违背约束），使得优化结果比worst-case具有更好的结果。

例如，对于优化问题（6），优化目标是最差的情况worst-case，在此可以将目标函数松弛25%，也即是优化最差的25%部分对应的分位点。不是很好理解，使用下图说明一下。

如图所示，使用N=10作为一个示意，对$x_{p,i}$进行升序排序。在（a）图中worst-case情况下，$\omega$应当小于所有的$x_{p,i}$，此时的$\omega$为最左边的$x_{p,i}$。现在希望针对最差的25%部分进行优化，也即是$w$在原来的基础上向右平移了25%*N，变成了(b)图箭头指示部分，此时即是优化该位置的$\omega$。
写成数学模型即是：
max ⁡ y ω s . t .    P r ( ω ≥ x p , i ) ≤ 25 % , ∀ i ∈ { 1 , ⋯   , N } 0 ≤ x s , i ≤ d i x s , i + x d , i = y x p , i = c s x s , i − c o y + c d x d , i x d , i ≥ 0 (7) \begin{aligned} & \max_{y} \omega \\ & s.t. \; \rm{Pr}( \omega \geq x_{p,i})\leq 25\%,\forall \it i \in \{1,\cdots,N \} \\ & \quad 0 \leq x_{s,i}\leq d_{i} \\ & \quad x_{s,i}+x_{d,i}=y \\ & \quad x_{p,i}=c_{s}x_{s,i}-c_{o}y+c_{d}x_{d,i} \\ & \quad x_{d,i}\geq 0 \end{aligned}\tag{7} ​ymax​ωs.t.Pr(ω≥xp,i​)≤25%,∀i∈{1,⋯,N}0≤xs,i​≤di​xs,i​+xd,i​=yxp,i​=cs​xs,i​−co​y+cd​xd,i​xd,i​≥0​(7)

优化问题（7）中的第一个约束条件是机会约束条件或者叫概率约束条件，
其中，$\mathrm{P}\mathrm{r}(\omega \ge {\mathrm{x}}_{\mathrm{p},\mathrm{i}})\le 25\%=\mathrm{P}\mathrm{r}(\omega \le {\mathrm{x}}_{\mathrm{p},\mathrm{i}})\ge 1-25\%$。

实际编程求解时需要对问题（7）进一步转化：
max ⁡ y ω s . t .    ω ≤ x p , i + B i ( x p ‾ − x p ‾ ) ∑ i = 1 N B i N ≤ 25 % , ∀ B i ∈ { 0 , 1 } , i ∈ { 1 , ⋯   , N } 0 ≤ x s , i ≤ d i x s , i + x d , i = y x p , i = c s x s , i − c o y + c d x d , i x d , i ≥ 0 (8) \begin{aligned} & \max_{y} \omega \\ & s.t. \; \omega \leq x_{p,i}+ B_{i}(\overline{x_{p}}-\underline{x_{p}})\\ & \quad \sum_{i=1}^{N}\frac{B_{i}}{N}\leq 25\%,\forall B_{i} \in \{0,1\}, i \in \{1,\cdots,N \} \\ & \quad 0 \leq x_{s,i}\leq d_{i} \\ & \quad x_{s,i}+x_{d,i}=y \\ & \quad x_{p,i}=c_{s}x_{s,i}-c_{o}y+c_{d}x_{d,i} \\ & \quad x_{d,i}\geq 0 \end{aligned}\tag{8} ​ymax​ωs.t.ω≤xp,i​+Bi​(xp​​−xp​​)i=1∑N​NBi​​≤25%,∀Bi​∈{0,1},i∈{1,⋯,N}0≤xs,i​≤di​xs,i​+xd,i​=yxp,i​=cs​xs,i​−co​y+cd​xd,i​xd,i​≥0​(8)
问题（8）的第一二行约束等价于问题（7）的第一个约束，其中，$\overline{x_p}=\max x_{p,i},\underline{x_p}=\min x_{p,i}$,$B_i$为引入的0-1变量。对于$\overline{x_p}$可以使用程序中的maxrev近似，而$\underline{x_p}$可以使用minrev近似。问题（8）为混合整数线性规划（mixed integer linear programming, MILP）问题，求解十分复杂，耗时很长。

对于问题（8）中的第一二约束的理解可以使用上图中的(b)图。在第一个约束中有大于75%的约束满足$\omega \le x_{p,i}$，有小于25%的约束满足$\omega \le x_{p,i}+(\overline{x_p}-\underline{x_p})$。为使$\omega$具有最大值，应当在最小的25%$x_{p,i}$上加上$\overline{x_p}-\underline{x_p}$，也即是（b）图中红蓝相间的线条。此时的$\omega$就近似为$x_{p,i}$的25%分为点的数值。这种转化方法很巧妙，值得学习。

求解代码如下：

# Maximize the 25% 25% worst-case profit
mo = Model()
# Set to maximize
mo.ModelSense = -1
# Add variables
worst    = mo.addVar(lb=minrev,ub=maxrev,obj=1,name='worst') #obj目标函数中的系数
order    = mo.addVar(name='order')
chance   = mo.addVars(samples,vtype='B',name='chance') #添加的0-1二进制变量
profit   = mo.addVars(samples,lb=minrev,ub=maxrev,name='profit')
sales    = mo.addVars(samples,ub=demand,name='sales')
discount = mo.addVars(samples,name='discount')
# Set constraints
mo.addConstrs((profit[i] == -order * cost + sales[i] * retail + recover * discount[i] for i in range(samples)),name='profit')
mo.addConstrs((sales[i]+discount[i] == order for i in range(samples)),name='demand')
mo.addConstrs((worst - (maxrev-minrev)*chance[i] <= profit[i] for i in range(samples)),name='worst') #最少75%的概率小于原始约束
mo.addConstr(chance.sum() <= samples*0.25)
mo.update()
mo.params.TimeLimit=1200 #设置求解时间门限1200秒
mo.optimize()
print("25% worst-case进货量:",order.x)
print("25% worst-case总体利润:",mo.objVal)

结果如下图，25%worst-case 情况下，目标函数为3307.8，但是求解用时为20分钟。：

(2) CVaR

由于以上机会约束引入了0-1二进制变量，以上问题（8）为包含整数的混合整数线性规划问题，求解十分复杂，耗时很长。在此使用简单的线性规划进行近似估计，根据以上gurubi的webinar，使用风险度量的指标进行衡量。

条件风险价值（Conditional value at risk，CVaR）[4]: 是指损失序列的分布中，升序排序后,排序最高的alpha%部分的均值。英文解释更容易理解（Conditional value at risk (for losses) is the expected value of the worst alpha% tail of the realizations in the random variable）。
其数学定义为：
C V a R α ( ξ ) = E [ ξ ∣ ξ ≥ V a R α ( ξ ) ] = ∫ α 1 V a R α ( ξ ) d α = min ⁡ { t + E [ ( ξ − t ) + ] 1 − α } (9) \begin{aligned} & CVaR_{\alpha}(\xi)=\Bbb{E}[\xi|\xi \geq VaR_{\alpha}(\xi)] \\ & \qquad = \int_{\alpha}^{1}VaR_{\alpha}(\xi)d\alpha \\ & \qquad = \min\{t + \frac{\Bbb{E}[(\xi -t)^{+}]} {1-\alpha}\} \end{aligned} \tag{9} ​CVaRα​(ξ)=E[ξ∣ξ≥VaRα​(ξ)]=∫α1​VaRα​(ξ)dα=min{t+1−αE[(ξ−t)+]​}​(9)
其中， ( ξ − t ) + = max ⁡ { ξ − t , 0 } (\xi -t)^{+}=\max \{ \xi-t,0\} (ξ−t)+=max{ξ−t,0}，第3个等式没太弄明白，$t$应该也是指总体利润。

而需要注意的是，**风险价值VaR(value at risk)**是指在一定置信水平$1-\alpha$下，变量$\xi$面临的最大损失:

V a R α ( ξ ) = { t : P r ( ξ ≤ t ) = 1 − α } VaR_{\alpha}(\xi)=\{t:Pr(\xi \leq t)=1-\alpha\} VaRα​(ξ)={t:Pr(ξ≤t)=1−α}
注意到在此，$CVa{R}_{\alpha }(\xi )$和$Va{R}_{、alpha}(\xi )$都是针对损失（Loss）而言的。

CVaR的计算很简单，先对损失数据进行排序，找到分位点，再分位点之后的数据计算均值就是对应的CVaR值。

#CVaR函数
def CVaR(data,alpha):
    data.sort()
    n = len(data)
    m = int(alpha*n)
    return sum(data[m:-1])/len(data[m:-1])

在此使用利润$x_{p,i}的负值，即$$Loss=-x_{p,i}$，作为损失。因而损失越大，收益越小。令$\alpha =75\%$，则根据公式（9），损失的75%CVaR计算的是损失大于75%VaR部分的均值。而对应的大于75%的部分，即是损失$Loss$最大的25%部分，也即是收益$x_{p,i}$最差的25%。

因此最小化$Loss$的75%CVaR,也即是最大化收益$x_{p,i}$最差的25%部分的均值。也即是用均值近似问题（8）中的机会约束。

具体求解代码如下：

# Maximize CVaR_0.75 profit
alpha = 0.75
m = Model()
# Set to maximize
m.ModelSense = -1
# Add variables
t        = m.addVar(lb=minrev,ub=maxrev,obj=-1,name='worst') # t应该也是指总利润
order    = m.addVar(name='order')
excess   = m.addVars(samples,obj=-1.0/((1-alpha)*samples),ub=maxrev-minrev,name='chance') #
# excess[i] 为新引入的变量，代表公式(9)期望中的 (\xi-t)^+部分
profit   = m.addVars(samples,lb=minrev,ub=maxrev,name='profit')
sales    = m.addVars(samples,ub=demand,name='sales')
discount = m.addVars(samples,name='discount')
# Set constraints
m.addConstrs((profit[i] == -order * cost + sales[i] * retail + recover * discount[i] for i in range(samples)),name='profit')
m.addConstrs((sales[i]+discount[i] == order for i in range(samples)),name='demand')
m.addConstrs((-profit[i]-t <= excess[i] for i in range(samples)),name='excess')
# -profit[i]即是损失Loss=-x_{p,i} 这一个约束表达的是 \xi - t <= excess[i]
m.update()
m.optimize()
print("CVaR 进货量:",order.x)
print("CVaR 总体利润:",m.objVal)

以上代码对应的优化问题为：
max ⁡ { − t − ∑ i = 1 N e x c e s s i 1 − α } = min ⁡ { t + ∑ i = 1 N ( ξ i − t ) + 1 − α } = min ⁡ C V a R α ( ξ ) (10) \begin{aligned} & \max \{ -t - \frac{\sum_{i=1}^{N}excess_{i}}{1-\alpha} \} \\ & =\min \{ t+\frac{\sum_{i=1}^{N}{(\xi_{i}-t)^+}}{1-\alpha} \} \\ & =\min CVaR_{\alpha}(\xi) \end{aligned} \tag{10} ​max{−t−1−α∑i=1N​excessi​​}=min{t+1−α∑i=1N​(ξi​−t)+​}=minCVaRα​(ξ)​(10)

其为线性规划问题，容易求解，结果如下，目标函数为3078.7，比优化问题（8）目标函数近小7%，而最优解十分接近。

3.总结

对于随机优化，要求随机变量具有确定概率分布，使用SAA方法其求解复杂度较低，结果可以接受。但是对于包含概率约束的场景，需要进行转化利用CVaR进行近似求解。需要注意的是，对于VaR的定义比较混乱，可以从数据统计值直方图的左尾定义，也可以从右尾定义，容易搞混，理解很麻烦。对于gurobi的软件许可和相关使用教程可以访问gurobi中国网站（http://www.gurobi.cn/）。

参考文献

[1] Kleywegt Anton, et al. The sample average approximation method for stochastic discrete optimization [J], society for industrial applied mathematics, 2001.

[2] Lauren A. Hannah. Stochastic Optimization. April 4, 2014. http://www.stat.columbia.edu/~liam/teaching/compstat-spr14/lauren-notes.pdf

[3] Gurobi webinar的视频. Stochastic Programing Part I: Why should we care about uncertainty_哔哩哔哩 (゜-゜)つロ 干杯~-bilibili[EB/OL]. [2020-9-4]. https://www.bilibili.com/video/BV1Lp4y1Q7RN/?spm_id_from=333.788.videocard.0

[4] 条件风险价值CVaR_zte10096334的博客-CSDN博客[EB/OL]. [2020-9-4]. https://blog.csdn.net/zte10096334/article/details/94761411.