隨機優化中的樣本均值近似方法

• 1.SAA原理
• 2.使用Gurobi實操求解
• • 2.1 問題求解
  • 2.2 最壞的情況（worst case）
  • 2.3 機會約束優化（chance constrained optimization）
  • • (1) 25% worst-cast
    • (2) CVaR
• 3.總結

補一下之前隨機優化的學習內容，結合Gurobi的webinar中隨機優化的講座和gurobi求解器進行求解。

1.SAA原理

隨機優化是求解包含隨機變量的一類優化問題，其中隨機變量的分布是確定已知的。而樣本均值近似（sample average approximation,SAA）方法是求解隨機優化的常用方法，其原理部分主要參考翻譯文獻[1]，其中關于SAA的收斂速率的證明，有興趣的同學可以讀讀原文獻，數學太難，在此省略。

考慮以下包含期望的隨機優化問題：
$$\underset{x\in \mathbf{S}}{\min }g(x)={\mathbb{E}}_{W～\mathbf{P}}G(x,W)\text{\hspace{1em}(1)}$$

其中，$W$為隨機變量，其概率分布服從$W～\mathbf{P}$。優化變量屬于有限集合（finite set）$x\in \mathbf{S}$。$w$是隨機變量$W$的樣本，$G(x,w)$是$G(x,W)$在樣本$w$下的一次實現。
式（1）中的期望${\mathbb{E}}_{W～\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].
問題建模如下:
$$\begin{array}{rl}& \max \mathbb{E}[x_p]\\ & s.t.0\le x_s\le \min \{y,d\}\\ & 0\le x_d\le \max \{0,y?x_s\}\\ & x_p=c_s{x}_s?c_{o}y+c_d{x}_d\end{array}\text{\hspace{1em}(3)}$$

其中,$x_p$ 為總體銷售利潤， $x_s$為銷量， $c_s$為售價， $y$為進貨量， $c_o$為成本， $x_d$為打折銷售的數量（滯銷的數量）， $c_d$為滯銷的售價，如果$c_d\ge 0$表示滯銷仍有收益，相反$c_d<0$表示賠本銷售。

需要注意的是$d$表示市場的需求量，為隨機變量，具有不確定性。而進貨量$y$為最終要求解的決策變量。

使用SAA方法求解：

$$\begin{array}{rl}& \max \sum _{i=1}^N\frac{x_{p,i}}{N},?i\in \{1,?,N\}\\ & s.t.0\le x_{s,i}\le \min \{y,d_i\}\\ & 0\le x_{d,i}\le \max \{0,y?x_{s,i}\}\\ & x_{p,i}=c_s{x}_{s,i}?c_{o}y+c_d{x}_{d,i}\end{array}\text{\hspace{1em}(4)}$$

$d_i$為隨機變量第i次的抽樣值，$x_{p,i},x_{s,i},x_{d,i}$為在第i次抽樣對應的優化變量。而最終的求解的變量y在每次抽樣中都保持不變。此時的問題（4）轉化為包含$3?N+1$個變量的確定性優化問題。

2.1 問題求解

進一步，優化問題（4）可以轉化為：
$$\begin{array}{rl}& \max \sum _{i=1}^N\frac{x_{p,i}}{N},?i\in \{1,?,N\}\\ & s.t.0\le x_{s,i}\le 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}\ge 0\end{array}\text{\hspace{1em}(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）

考察一下，最壞的情況。假設用戶希望即便是在利潤最差的情況下也能有較好的結果，這既是魯棒優化問題。此時的目標函數為：
$$\underset{y}{\max }\{\min x_p\}$$
進一步，建立的優化問題為：
$$\begin{array}{rl}& \underset{y}{\max }\omega \\ & s.t.\omega \le x_{p,i},?i\in \{1,?,N\}\\ & 0\le x_{s,i}\le 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}\ge 0\end{array}\text{\hspace{1em}(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$。
寫成數學模型即是：
$$\begin{array}{rl}& \underset{y}{\max }\omega \\ & s.t.\mathrm{P}\mathrm{r}(\omega \ge {\mathrm{x}}_{\mathrm{p},\mathrm{i}})\le 25\%,?i\in \{1,?,N\}\\ & 0\le x_{s,i}\le 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}\ge 0\end{array}\text{\hspace{1em}(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）進一步轉化：
$$\begin{array}{rl}& \underset{y}{\max }\omega \\ & s.t.\omega \le x_{p,i}+B_i(\overline{x_p}?\underline{x_p})\\ & \sum _{i=1}^N\frac{B_i}{N}\le 25\%,?B_i\in \{0,1\},i\in \{1,?,N\}\\ & 0\le x_{s,i}\le 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}\ge 0\end{array}\text{\hspace{1em}(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）。
其數學定義為：
$$\begin{array}{rl}& CVa{R}_{\alpha }(\xi )=\mathbb{E}[\xi |\xi \ge Va{R}_{\alpha }(\xi )]\\ & ={\int }_{\alpha }^{1}Va{R}_{\alpha }(\xi )d\alpha \\ & =\min \{t+\frac{\mathbb{E}[(\xi ?t{)}^{+}]}{1?\alpha }\}\end{array}\text{\hspace{1em}(9)}$$
其中，$(\xi ?t{)}^{+}=\max \{\xi ?t,0\}$，第3個等式沒太弄明白，$t$應該也是指總體利潤。

而需要注意的是，**風險價值VaR(value at risk)**是指在一定置信水平$1?\alpha$下，變量$\xi$面臨的最大損失:

$$Va{R}_{\alpha }(\xi )=\{t:Pr(\xi \le t)=1?\alpha \}$$
注意到在此，$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)

以上代碼對應的優化問題為：
$$\begin{array}{rl}& \max \{?t?\frac{{\sum }_{i=1}^{N}exces{s}_i}{1?\alpha }\}\\ & =\min \{t+\frac{{\sum }_{i=1}^N({\xi }_i?t{)}^{+}}{1?\alpha }\}\\ & =\min CVa{R}_{\alpha }(\xi )\end{array}\text{\hspace{1em}(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.

總結

以上是生活随笔為你收集整理的随机优化中的样本均值近似方法的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。