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

前言

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

1. RSO

文献[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]。

2. 简单示例

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

2.1 单产品的报童问题(one-product newsvendor problem)

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

通常情况下，当天的市场需求量$\tilde{u}$对于经销商(报童)而言是未知的，其决策进货量$w$是在随机变量揭示自己之前做出的。报童的目标当然是最大化经销利润。由于存在不确定性的随机变量$\tilde{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\ge 0$$
其中，$\inf (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$和0的最大值。

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

$$\begin{gathered} \underset{w}{\max }\underset{\mathbb{P}\in \mathcal{F}}{\inf }{\mathbb{E}}_{\mathbb{P}}[p\cdot [w-(w-\tilde{u}{)}^{+}]-c\cdot w] \\ =\underset{w}{\max }w\cdot (p-c)+\underset{\mathbb{P}\in \mathcal{F}}{\inf }{\mathbb{E}}_{\mathbb{P}}[-p\cdot (w-\tilde{u}{)}^{+}]\text{\hspace{1em}(2)} \end{gathered}$$

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

$$\begin{gathered} \underset{w}{\max }w\cdot (p-c)+\underset{\mathbb{P}\in \mathcal{F}}{\inf }{\mathbb{E}}_{\mathbb{P}}[-p\cdot (w-\tilde{u}{)}^{+}] \\ =\underset{w}{\max }w\cdot (p-c)-\underset{\mathbb{P}\in \mathcal{F}}{\sup }{\mathbb{E}}_{\mathbb{P}}[p\cdot (w-\tilde{u}{)}^{+}]\text{\hspace{1em}(3)} \end{gathered}$$

最终优化问题(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\ge 0$$

2.2 Wasserstein 模糊集

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

$$\mathcal{F}=\left\{\begin{array}{ccc}& |& (\tilde{u},\tilde{s})\sim \mathbb{P}\\ & |& {\mathbb{E}}_p[\rho (\tilde{u},{\hat{u}}_{\tilde{s}})|\tilde{s}\in [S]]\le \theta \\ \mathbb{P}\in {\mathcal{P}}_0(\mathbb{R}\times [S])& |& \mathbb{P}[\tilde{u}\in [0,\bar{U}]|\tilde{s}=s]=1,\forall s\in [S]\\ & |& \mathbb{P}[\tilde{s}=s]=1/S,\forall s\in [S]\end{array}\right\}$$

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

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

$$\mathcal{F}=\left\{\begin{array}{ccc}& |& ((\tilde{u},\tilde{v}),\tilde{s})\sim \mathbb{P}\\ & |& {\mathbb{E}}_p[\tilde{v}|\tilde{s}\in [S]]\le \theta \\ \mathbb{P}\in {\mathcal{P}}_0({\mathbb{R}}^2\times [S])& |& \mathbb{P}\left[\begin{array}{ccc}\tilde{u}\in [0,\bar{U}]& |& \tilde{s}=s\\ \rho (\tilde{u},{\hat{u}}_{\tilde{s}})\le \tilde{v}& |& \end{array}\right]=1,\forall s\in [S]\\ & |& \mathbb{P}[\tilde{s}=s]=1/S,\forall s\in [S]\end{array}\right\}$$

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

2.3 测试求解

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

Ubar = 100;                             % Upper bounds of random demands
S = 500;                                % Number of samples 
Uhat = Ubar * rand(1, S);               % Empirical data of demands

p = 1.5;                                % Selling price of the product
c = 1.0;                                % Production cost of the product
theta = Ubar*0.01;                      % The Wasserstein distance

%% Create a RSOME model
model = rsome('newsvendor');            % Create a model object named "newsvendor"

%% Random variables and a type-1 Wasserstein ambiguity set
u = model.random;                       % Random demand
v = model.random;                       % Auxiliary random variable 
P = model.ambiguity(S);                 % Create an ambiguity set with S scenarios
for n = 1:S
    P(n).suppset(0 <= u, u <= Ubar, ...   
                 norm(u-Uhat(n)) <= v );% Define the support set for each scenario
end
P.exptset(expect(v) <= theta);
prob = P.prob;                          % pr for all scenario probabilities
P.probset(prob == 1/S);                 % The probability set of the ambiguity set
model.with(P);                          % The ambiguity set of the model is P

%% Decision
w = model.decision;                     % Define a decision variable w

%% Objective function
loss = maxfun({p*(w-u), 0});            % Profit loss due to unsold products
model.max((p-c)*w - expect(loss));      % Objective as the worst-case expectation

%% Constraints
model.append(w >= 0);                   % The variable w is non-negative

%% Solution
model.solve;                            % Solve the model

我的运行环境是MATLAB2018a和cplex12.8教育版。由于cplex的教育版最大支持1000个优化变量的优化问题，在运行example_3_3_2_exact_nv.m文件时，rsome类的solve方法在调用Cplex求解时，报Cplex/subsref错误。

在这个例子中，样本数量S=500，导致变量超出了cplex限制，报了下标(subsref)的错误。我把S改为S=70运行正确。如果大家求解问题规模过大，还得需要完全、无限制版的优化求解器。
最后结果为$w=49.5177$。

总结

RSOME在使用体验上同yalmip十分相似，但是需要对DRO具有一定的理解才能真正入手。这只是其中的入门的例子，本人数学较差，后面还得补下RO和DRO的基础知识。在此十分感谢RSOME作者香港城市大学副教授Chen Zhi的答疑和其提供的国内网盘RSOME下载链接。

参考文献
[1]Math for DS 直播 NO.4 | 斯坦福大学管理科学与工程叶荫宇. https://www.bilibili.com/video/BV1Kt4y1y7ZU
[2]Chen, Zhi & Sim, Melvyn & Xiong, Peng. (2020). Robust Stochastic Optimization Made Easy with RSOME. Management Science. DOI: 10.1287/mnsc.2020.3603.
[3]Zhi Chen, Melvyn Sim, Peng Xiong.Users Guide for RSOME, version 1.2.
[4]链接: https://pan.baidu.com/s/1ejWU-lRwsmmfo5YgrG5hyA 提取码: k37r