随机优化中的样本均值近似方法
补一下之前随机优化的学习内容,结合Gurobi的webinar中随机优化的讲座和gurobi求解器进行求解。
1.SAA原理
随机优化是求解包含随机变量的一类优化问题,其中随机变量的分布是确定已知的。而样本均值近似(sample average approximation,SAA)方法是求解随机优化的常用方法,其原理部分主要参考翻译文献[1],其中关于SAA的收敛速率的证明,有兴趣的同学可以读读原文献,数学太难,在此省略。
考虑以下包含期望的随机优化问题:
min
x
∈
S
g
(
x
)
=
E
W
∼
P
G
(
x
,
W
)
(1)
\min_{x \in \bf{S}}{g(x)=\Bbb{E}_{W \sim \bf{P}}G(x,W)} \tag{1}
x∈Sming(x)=EW∼PG(x,W)(1)
其中,
W
W
W为随机变量,其概率分布服从
W
∼
P
W \sim \bf{P}
W∼P。优化变量属于有限集合(finite set)
x
∈
S
x \in \bf{S}
x∈S。
w
w
w是随机变量
W
W
W的样本,
G
(
x
,
w
)
G(x,w)
G(x,w)是
G
(
x
,
W
)
G(x,W)
G(x,W)在样本
w
w
w下的一次实现。
式(1)中的期望
E
W
∼
P
G
(
x
,
W
)
=
∫
G
(
x
,
W
)
f
(
W
)
d
W
\Bbb{E}_{W \sim \bf{P}}G(x,W)=\int G(x,W)f(W)\ dW
EW∼PG(x,W)=∫G(x,W)f(W) dW是关于随机变量
W
W
W的积分,
f
(
W
)
f(W)
f(W)是
W
W
W的概率密度函数。
使用SAA的前提[1]
使用SAA需要满足以下三个前提假设:
- (1) 期望函数 g ( x ) g(x) g(x)没有闭合表达式,并且其数值不容易计算
- (2) 函数 G ( x , w ) G(x,w) G(x,w)对于给定的 x x x和 w w w是容易计算的
- (3) 可行解集虽是有限的,但也是数量足够大的,不能使用穷举方法。
定义函数
g
(
x
)
g(x)
g(x)的样本均值近似为:
g
^
N
(
x
)
=
1
N
∑
j
=
1
N
G
(
x
,
w
j
)
\hat{g}_{N}(x)=\frac{1}{N}\sum_{j=1}^{N}G(x,w_{j})
g^N(x)=N1j=1∑NG(x,wj)
其中,
w
j
w_{j}
wj为
W
W
W的第
j
j
j次采样值。
相应的优化问题转化为:
min
x
∈
S
g
^
N
(
x
)
(2)
\min_{x \in \bf{S}}{\hat{g}_{N}(x)} \tag{2}
x∈Sming^N(x)(2)
把优化问题(1)称为真问题(原始问题),优化问题(2)称为样本均值近似问题,需要注意到
g
^
N
(
x
)
\hat{g}_{N}(x)
g^N(x)是
g
(
x
)
g(x)
g(x)的无偏估计量,即
E
[
g
^
N
(
x
)
]
=
g
(
x
)
\Bbb{E}[\hat{g}_{N}(x)]=g(x)
E[g^N(x)]=g(x) 。
SAA两个重要性质[2]:
- 渐进收敛性:随着样本数量N趋于无穷大,问题(2)的最优解和最优值和收敛于原始问题(1)的最优解和最优值。
- 易处理性:对于大多数函数 G ( x , W ) G(x,W) G(x,W)和可行解集 S \bf{S} 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=csxs−coy+cdxd(3)
其中, x p x_{p} xp 为总体销售利润, x s x_{s} xs为销量, c s c_{s} cs为售价, y y y为进货量, c o c_{o} co为成本, x d x_{d} xd为打折销售的数量(滞销的数量), c d c_{d} cd为滞销的售价,如果 c d ≥ 0 c_{d}\geq 0 cd≥0表示滞销仍有收益,相反 c d < 0 c_{d}<0 cd<0表示赔本销售。
需要注意的是 d d d表示市场的需求量,为随机变量,具有不确定性。而进货量 y y 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∑NNxp,i,∀i∈{1,⋯,N}s.t.0≤xs,i≤min{y,di}0≤xd,i≤max{0,y−xs,i}xp,i=csxs,i−coy+cdxd,i(4)
d i d_{i} di为随机变量第i次的抽样值, x p , i , x s , i , x d , i x_{p,i},x_{s,i},x_{d,i} xp,i,xs,i,xd,i为在第i次抽样对应的优化变量。而最终的求解的变量y在每次抽样中都保持不变。此时的问题(4)转化为包含 3 ∗ N + 1 3*N+1 3∗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∑NNxp,i,∀i∈{1,⋯,N}s.t.0≤xs,i≤dixs,i+xd,i=yxp,i=csxs,i−coy+cdxd,ixd,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
N
N变化对SAA方法的影响,如下图所示。设置最大的样本数量为10000个,每次迭代增加500个样本,可见随着样本数量增加,优化变量
y
y
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≤dixs,i+xd,i=yxp,i=csxs,i−coy+cdxd,ixd,i≥0(6)
其中,第一个约束条件保证
ω
\omega
ω为
x
p
,
i
x_{p,i}
xp,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
x_{p,i}
xp,i进行升序排序。在(a)图中worst-case情况下,
ω
\omega
ω应当小于所有的
x
p
,
i
x_{p,i}
xp,i,此时的
ω
\omega
ω为最左边的
x
p
,
i
x_{p,i}
xp,i。现在希望针对最差的25%部分进行优化,也即是
w
w
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≤dixs,i+xd,i=yxp,i=csxs,i−coy+cdxd,ixd,i≥0(7)
优化问题(7)中的第一个约束条件是机会约束条件或者叫概率约束条件,
其中,
P
r
(
ω
≥
x
p
,
i
)
≤
25
%
=
P
r
(
ω
≤
x
p
,
i
)
≥
1
−
25
%
\rm{Pr}( \omega \geq x_{p,i})\leq 25\%=\rm{Pr}( \omega \leq x_{p,i})\geq 1-25\%
Pr(ω≥xp,i)≤25%=Pr(ω≤xp,i)≥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∑NNBi≤25%,∀Bi∈{0,1},i∈{1,⋯,N}0≤xs,i≤dixs,i+xd,i=yxp,i=csxs,i−coy+cdxd,ixd,i≥0(8)
问题(8)的第一二行约束等价于问题(7)的第一个约束,其中,
x
p
‾
=
max
x
p
,
i
,
x
p
‾
=
min
x
p
,
i
\overline{x_{p}}=\max x_{p,i},\underline{x_{p}}=\min x_{p,i}
xp=maxxp,i,xp=minxp,i,
B
i
B_{i}
Bi为引入的0-1变量。对于
x
p
‾
\overline{x_{p}}
xp可以使用程序中的maxrev近似,而
x
p
‾
\underline{x_{p}}
xp可以使用minrev近似。问题(8)为混合整数线性规划(mixed integer linear programming, MILP)问题,求解十分复杂,耗时很长。
对于问题(8)中的第一二约束的理解可以使用上图中的(b)图。在第一个约束中有大于75%的约束满足 ω ≤ x p , i \omega \leq x_{p,i} ω≤xp,i,有小于25%的约束满足 ω ≤ x p , i + ( x p ‾ − x p ‾ ) \omega \leq x_{p,i}+ (\overline{x_{p}}-\underline{x_{p}}) ω≤xp,i+(xp−xp)。为使 ω \omega ω具有最大值,应当在最小的25% x p , i x_{p,i} xp,i上加上 x p ‾ − x p ‾ \overline{x_{p}}-\underline{x_{p}} xp−xp,也即是(b)图中红蓝相间的线条。此时的 ω \omega ω就近似为 x p , i x_{p,i} xp,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α(ξ)]=∫α1VaRα(ξ)dα=min{t+1−αE[(ξ−t)+]}(9)
其中,
(
ξ
−
t
)
+
=
max
{
ξ
−
t
,
0
}
(\xi -t)^{+}=\max \{ \xi-t,0\}
(ξ−t)+=max{ξ−t,0},第3个等式没太弄明白,
t
t
t应该也是指总体利润。
而需要注意的是,**风险价值VaR(value at risk)**是指在一定置信水平 1 − α 1-\alpha 1−α下,变量 ξ \xi ξ面临的最大损失:
V
a
R
α
(
ξ
)
=
{
t
:
P
r
(
ξ
≤
t
)
=
1
−
α
}
VaR_{\alpha}(\xi)=\{t:Pr(\xi \leq t)=1-\alpha\}
VaRα(ξ)={t:Pr(ξ≤t)=1−α}
注意到在此,
C
V
a
R
α
(
ξ
)
CVaR_{\alpha}(\xi)
CVaRα(ξ)和
V
a
R
、
a
l
p
h
a
(
ξ
)
VaR_{、alpha}(\xi)
VaR、alpha(ξ)都是针对损失(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 的 负 值 , 即 x_{p,i}的负值,即 xp,i的负值,即 L o s s = − x p , i Loss=-x_{p,i} Loss=−xp,i,作为损失。因而损失越大,收益越小。令 α = 75 % \alpha=75\% α=75%,则根据公式(9),损失的75%CVaR计算的是损失大于75%VaR部分的均值。而对应的大于75%的部分,即是损失 L o s s Loss Loss最大的25%部分,也即是收益 x p , i x_{p,i} xp,i最差的25%。
因此最小化 L o s s Loss Loss的75%CVaR,也即是最大化收益 x p , i x_{p,i} xp,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=1Nexcessi}=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.