鲁棒随机优化（Robust Stochastic Optimization）和RSOME

记录一下鲁棒优化的学习内容。鲁棒优化和随机优化是不确定性优化中常用的两种方法。其中随机优化的随机变量的概率分布是确定的，对于不确定性观测而言其，确定性的概率分布太乐观了。而对于鲁棒优化是优化不确定下的worst-case的情形，而往往实际的结果没有这么差，也即是其结果太悲观了。因此在随机优化和鲁棒优化之间，产生了概率分布鲁棒优化（Distributionally Robust Optimization,DRO）(不知翻译是否准确),它结合了随机优化的期望目标函数，扩展了不确定性参数的概率分布集合，同时利用了鲁棒优化的worst-case优化原则。相关的DRO基础内容，可以参考斯坦福叶荫宇教授关于DRO的Talk[1]。

文献[2]把树形场景（tree-scenario）的随机线性规划（stochastic linear optimization,SLO）和DRO的不确定性集合，使用event-wise模糊集（event-wise ambiguity set）统一起来，把他们统称为鲁棒随机优化。其每一个离散的场景（或者采样）都对应一个独立完备的不确定性优化问题，也即是一个不确定场景的实例（instance）。 RSO或者DRO相较经典RO的优势在于优化结果大于等于绝对悲观的情景，并且能够结合机器学习方法，利用数据驱动（data-driven）进行优化建模，具有更好的说服力，例如k均值模糊集（K-means Ambiguity Set）。

难能可贵的是文献[2]提供了一个对应RSO的Matlab工具包RSOME（Robust Stochastic Optimization Made Easy），其原码和用户手册可以在它们的官网上下载( www.rsomerso.com ）。

工具包中提供了所有的manual examples的源代码，十分方便大家学习。运行环境需要大家安装相应的优化求解器例如cplex，gurobi或者MOSEK。作者论文[2]附录中使用的求解器是MOSEK，但是源代码默认采用cplex作为求解器。RSOME的安装和使用参考其用户手册[3]。

数学公式太复杂，来看一下用户手册中的一个简单例子。

以用户指南16页的3.3.2单产品的报童问题(one-product newsvendor problem)为例[3]。定义 p p p为报纸售价(price)， c c c为报刊成本(cost)， w w w为报纸进货量， u ~ \tilde{u} u~为每天的市场需求量，是随机变量，服从 u ~ ∼ P ,    P ⊆ F \tilde{u} \sim \Bbb{P},\; \Bbb{P} \subseteq \mathscr{F} u~∼P,P⊆F， F \mathscr{F} F是概率分布的 P \Bbb{P} P的不确定集合，即模糊集。

通常情况下，当天的市场需求量 u ~ \tilde{u} u~对于经销商(报童)而言是未知的，其决策进货量 w w w是在随机变量揭示自己之前做出的。报童的目标当然是最大化经销利润。由于存在不确定性的随机变量 u ~ \tilde{u} u~，DRO优化的目标为即便是最差情况下也有最大的平均利润。 即： max ⁡ w    inf ⁡ P ∈ F E P [ p ⋅ min ⁡ { w , u ~ } − c ⋅ w ] (1) \underset{w}{\max} \; \underset{\Bbb{P} \in \mathscr{F}}{\inf}\Bbb{E_{\Bbb{P}}} [p \cdot \min\{w,\tilde{u}\}-c\cdot w]\tag{1} wmax​P∈Finf​EP​[p⋅min{w,u~}−c⋅w](1) s . t . w ≥ 0 s.t. \quad w \ge0 s.t.w≥0 其中， inf ⁡ ( x ) \inf(x) inf(x)为 x x x的下确界， min ⁡ { w , u ~ } \min\{w,\tilde{u}\} min{w,u~}为当天的销售量。

对于最大化和最小化函数 max ⁡ { x , y } 和 min ⁡ { x , y } \max\{x,y\}和\min\{x,y\} max{x,y}和min{x,y}，可以进行如下等效转换：

max ⁡ { x , y } = x + ( y − x ) + \max\{x,y\}=x+(y-x)^{+} max{x,y}=x+(y−x)+ min ⁡ { x , y } = x − ( x − y ) + \min\{x,y\}=x-(x-y)^{+} min{x,y}=x−(x−y)+

其中， ( x ) + = max ⁡ { x , 0 } (x)^{+}=\max\{x,0\} (x)+=max{x,0}即取 x x x和0的最大值。

由于参数 p , c p,c p,c为确定性的值， w w w先于 u ~ \tilde{u} u~确定，把以上等效转换代入公式(1)目标函数可得：

max ⁡ w    inf ⁡ P ∈ F E P [ p ⋅ [ w − ( w − u ~ ) + ] − c ⋅ w ] = max ⁡ w    w ⋅ ( p − c ) + inf ⁡ P ∈ F E P [ − p ⋅ ( w − u ~ ) + ] (2) \underset{w}{\max} \; \underset{\Bbb{P} \in \mathscr{F}}{\inf}\Bbb{E_{\Bbb{P}}} [p \cdot [w-(w-\tilde{u})^{+}]-c\cdot w] \\ \quad =\underset{w}{\max} \; w\cdot(p-c) + \underset{\Bbb{P} \in \mathscr{F}}{\inf}\Bbb{E_{\Bbb{P}}}[-p \cdot (w-\tilde{u})^{+}]\tag{2} wmax​P∈Finf​EP​[p⋅[w−(w−u~)+]−c⋅w]=wmax​w⋅(p−c)+P∈Finf​EP​[−p⋅(w−u~)+](2)

同样，由于 − sup ⁡ ( x ) = inf ⁡ ( − x ) -\sup(x)=\inf(-x) −sup(x)=inf(−x)， sup ⁡ ( x ) \sup(x) sup(x)为 x x x的上确界，公式(2)可以转化为：

max ⁡ w    w ⋅ ( p − c ) + inf ⁡ P ∈ F E P [ − p ⋅ ( w − u ~ ) + ] = max ⁡ w    w ⋅ ( p − c ) − sup ⁡ P ∈ F E P [ p ⋅ ( w − u ~ ) + ] (3) \underset{w}{\max} \; w\cdot(p-c) + \underset{\Bbb{P} \in \mathscr{F}}{\inf}\Bbb{E_{\Bbb{P}}}[-p \cdot (w-\tilde{u})^{+}]\\ \quad=\underset{w}{\max} \; w\cdot(p-c) - \underset{\Bbb{P} \in \mathscr{F}}{\sup}\Bbb{E_{\Bbb{P}}}[p \cdot (w-\tilde{u})^{+}] \tag{3} wmax​w⋅(p−c)+P∈Finf​EP​[−p⋅(w−u~)+]=wmax​w⋅(p−c)−P∈Fsup​EP​[p⋅(w−u~)+](3)

最终优化问题(1)转化为：

max ⁡ w    ( p − c ) w − sup ⁡ p ∈ F    ( E P [ max ⁡ { p ⋅ ( w − u ~ ) , 0 } ] ) (4) \underset{w}{\max}\;(p-c)w-\underset{\Bbb{p}\in \mathscr{F}}{\sup}\;(\Bbb{E}_{\Bbb{P}}[\max\{p \cdot (w-\tilde{u}),0\}]) \tag{4} wmax​(p−c)w−p∈Fsup​(EP​[max{p⋅(w−u~),0}])(4) s . t . w ≥ 0 s.t.\quad w \ge0 s.t.w≥0

在求解以上优化问题时，需要指定随机变量 u ~ \tilde{u} u~的模糊集 F \mathscr{F} F。例3.3.2中使用的是Wasserstein模糊集，其定义为：

F = { ∣ ( u ~ , s ~ ) ∼ P ∣ E p [ ρ ( u ~ , u ^ s ~ ) ∣ s ~ ∈ [ S ]    ] ≤ θ P ∈ P 0 ( R × [ S ] ) ∣ P [ u ~ ∈ [ 0 , U ˉ ] ∣ s ~ = s ] = 1 , ∀ s ∈ [ S ] ∣ P [ s ~ = s ] = 1 / S , ∀ s ∈ [ S ] } \mathscr{F}=\begin{Bmatrix} &|&(\tilde{u},\tilde{s})\sim \Bbb{P}\\ &|&\Bbb{E}_{\Bbb{p}}[\rho(\tilde{u},\hat{u}_{\tilde{s}})|\tilde{s} \in [S]\;] \leq \theta \\ \Bbb{P}\in \mathcal{P}_{0}(\Bbb{R}\times[S])&|&\Bbb{P}[\tilde{u}\in [0,\bar{U}]|\tilde{s}=s]=1,\forall s \in [S]\\ &|&\Bbb{P}[\tilde{s}=s]=1/S, \forall s \in [S] \end{Bmatrix} F=⎩⎪⎪⎨⎪⎪⎧​P∈P0​(R×[S])​∣∣∣∣​(u~,s~)∼PEp​[ρ(u~,u^s~​)∣s~∈[S]]≤θP[u~∈[0,Uˉ]∣s~=s]=1,∀s∈[S]P[s~=s]=1/S,∀s∈[S]​⎭⎪⎪⎬⎪⎪⎫​

其中， [ S ] = { 1 , 2 , ⋯   , S } [S]=\{1,2,\cdots,S\} [S]={1,2,⋯,S}为场景的集合， u ^ \hat{u} u^为 u ~ \tilde{u} u~的观测(采样)值或估计值。 ρ ( u ~ , u ^ s ~ ) \rho(\tilde{u},\hat{u}_{\tilde{s}}) ρ(u~,u^s~​)表示 u ~ \tilde{u} u~和 u ^ \hat{u} u^之间的某种距离测度，右侧第二行表示两者之间的平均距离小于给定的参数 θ \theta θ， θ \theta θ也被称为Wasserstein球半径。

根据文献[1]，通过引入松弛变量 v ~ {\tilde{v}} v~,可以把以上模糊集转化为event-wise模糊集：

F = { ∣ ( ( u ~ , v ~ ) , s ~ ) ∼ P ∣ E p [ v ~ ∣ s ~ ∈ [ S ]    ] ≤ θ P ∈ P 0 ( R 2 × [ S ] ) ∣ P [ u ~ ∈ [ 0 , U ˉ ] ∣ s ~ = s ρ ( u ~ , u ^ s ~ ) ≤ v ~ ∣ ] = 1 , ∀ s ∈ [ S ] ∣ P [ s ~ = s ] = 1 / S , ∀ s ∈ [ S ] } \mathscr{F}=\begin{Bmatrix} &|&((\tilde{u},\tilde{v}),\tilde{s})\sim \Bbb{P}\\ &|&\Bbb{E}_{\Bbb{p}}[\tilde{v}|\tilde{s} \in [S]\;] \leq \theta \\ \Bbb{P}\in \mathcal{P}_{0}(\Bbb{R}^{2}\times[S])&|&\Bbb{P} \begin{bmatrix}\tilde{u}\in [0,\bar{U}]&|&\tilde{s} = s\\ \rho(\tilde{u},\hat{u}_{\tilde{s}}) \leq \tilde{v} &|& \end{bmatrix}=1,\forall s \in [S]\\ &|&\Bbb{P}[\tilde{s}=s]=1/S, \forall s \in [S] \end{Bmatrix} F=⎩⎪⎪⎪⎪⎨⎪⎪⎪⎪⎧​P∈P0​(R2×[S])​∣∣∣∣​((u~,v~),s~)∼PEp​[v~∣s~∈[S]]≤θP[u~∈[0,Uˉ]ρ(u~,u^s~​)≤v~​∣∣​s~=s​]=1,∀s∈[S]P[s~=s]=1/S,∀s∈[S]​⎭⎪⎪⎪⎪⎬⎪⎪⎪⎪⎫​

也即是通过松弛变量 v ~ {\tilde{v}} v~来表示Wasserstein球的半径。

文献[3]的例3.3.2中，参数设置为: U ˉ = 100 , p = 1.5 , c = 1.0 \bar{U}=100,p=1.5,c=1.0 Uˉ=100,p=1.5,c=1.0，采样数量 S = 500 S=500 S=500，并假设每次采样值 u ^ \hat{u} u^服从 [ 0 , U ˉ ] [0,\bar{U}] [0,Uˉ]上的均匀分布， θ = 0.01 U ˉ \theta=0.01\bar{U} θ=0.01Uˉ。以下代码为RSOME用户手册中的代码[3]，非本人代码。