前言

记录一下鲁棒优化的学习内容。鲁棒优化和随机优化是不确定性优化中常用的两种方法。其中随机优化的随机变量的概率分布是确定的，对于不确定性观测而言其，确定性的概率分布太乐观了。而对于鲁棒优化是优化不确定下的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 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} wmaxP∈FinfEP[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} wmaxP∈FinfEP[p⋅[w−(w−u~)+]−c⋅w]=wmaxw⋅(p−c)+P∈FinfEP[−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} wmaxw⋅(p−c)+P∈FinfEP[−p⋅(w−u~)+]=wmaxw⋅(p−c)−P∈FsupEP[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

2.2 Wasserstein 模糊集

在求解以上优化问题时，需要指定随机变量 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球的半径。

2.3 测试求解

文献[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]，非本人代码。

Ubar = 100;                             % Upper bounds of random demands
S = 500;                                % Number of samples
Uhat = Ubar * rand(1, S);               % Empirical data of demandsp = 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:SP(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 w=49.5177 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

鲁棒随机优化（Robust Stochastic Optimization）和RSOME相关推荐

Robust Stochastic Optimization Made Easy with RSOME
摘要我们提出了一种新的分布式鲁棒优化模型,称为鲁棒随机优化 (RSO),它将基于场景树的随机线性优化和分布式鲁棒优化统一在一个可行的框架中,该框架可以使用最先进的商业优化求解器来解决 . 我们还开发 ...
拜占庭鲁棒随机聚合的分布式学习方法
文章目录拜占庭鲁棒随机聚合的分布式学习方法介绍分布式SGD 鲁棒分布式学习的RSA算法 l1l_1l1-范数 RSA p范数 RSA 结论原文链接:https://arxiv.org/abs ...
RoBERTa：一种鲁棒地优化BERT预训练的方法
RoBERTa:一种鲁棒地优化BERT预训练的方法文章目录 RoBERTa:一种鲁棒地优化BERT预训练的方法前言背景实验静态 VS 动态 Masking 输入形式与NSP任务更大的bat ...
语音合成（speech synthesis）方向二：鲁棒TTS（Robust TTS)
声明:工作以来主要从事TTS工作,工程算法都有涉及,平时看些文章做些笔记.文章中难免存在错误的地方,还望大家海涵.平时搜集一些资料,方便查阅学习:TTS 论文列表 http://yqli.tech/p ...
计及需求侧响应日前、日内两阶段鲁棒备用优化(Matlab代码实现）
计及需求侧响应日前、日内两阶段鲁棒备用优化【IEEE6节点】(Matlab代码实现）
鲁棒优化入门（一）——工具箱Xprog和RSOME的安装与使用
工具箱可以从这里下载:鲁棒优化工具箱Xprog和RSOME 一.Xprog 1.工具箱简介 Xprog是由新加坡国立大学的Peng Xiong于2016年发布的一款matlab工具箱,可以用于求解确定 ...
【鲁棒优化笔记】以Coding入门鲁棒优化：以一个例子引入（二）-错误版
[鲁棒优化笔记]以Coding入门鲁棒优化:以一个例子引入(二) 投资组合的例子鲁棒优化模型的reformulation: 利用对偶进行reformulation 利用对偶进行reformulati ...
电网两阶段鲁棒优化调度模型（含matlab程序）
目录一两阶段鲁棒优化理论二两阶段鲁棒优化程序实现 2.1 主/子问题变量要分清 2.2 对偶问题 2.3 线性化处理 2.4 编程小技巧 2.5 迭代问题三程序运行效果视频讲解两阶段鲁 ...

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