鲁棒优化｜利用rome求解鲁棒优化简单案例

super1 · 发表于 2022-5-25 15:35

实际问题中往往存在着许多不确定性，常见的处理不确定性问题的方法是随机规划、鲁棒优化等。一般我们都认为随机规划问题针对的是一个给定的概率分布，但实际中这个概率分布往往不能得到，通常只能通过历史数据去估计，但是当估计误差与真实分布相差较大时，随机规划的方案可能会带来极大的结果偏差。鲁棒优化不对不确定性变量做任何分布假设，而是通过构建不确定集，在其中寻找“最坏的情况最好”。
一、简单的股票投资案例

考虑一个由 $N$ 个股票构成的投资组合问题，我们认为，投资股票$i$带来的未来回报为 $r_i$ ，现在希望寻找一种投资组合策略，最大限度得提高总回报。这可以通过求解以下线性规划来实现： $\begin{aligned} \underset{x \in \mathbb{R}^{n}}{\operatorname{max}} \quad&\sum_{i=1}^{n} r_{i} x_{i} \\ \text { s.t. } &\sum_{i=1}^{n} x_{i}=100 \% \\ & 0 \leq x_i \leq 1 \quad \forall i \in\{1, \ldots, n\} \end{aligned}$
其中 $x_i$ 表示投资股票 $i$ 的比例，第一个约束意味着我们希望所有的预算都被投资，而 $x_i\ge 0$ 表示我们希望避免卖空股票。但实际上这个问题刻画是不符合实际的，因为真正的投资回报是不确定的，我们称之为不确定收益率，表述为：

$\bm{r}=\bm{\mu}+\bm{P}\bm{z}$
其中 $\mu$ 为收益的均值， $P$ 为收益的上下波动, 为不确定性集，  属于0到1，表示实际需求的波动程度，用不确定性集刻画。
二、CVaR不确定性集

风险价值（VaR）

VaR可以定义为在给定置信区间下可能发生的最大损失金额。在特定场景下，它也可以被视为最严重的损失。例：假设以95%的置信度每月VaR为1亿美元，其含义为在正常情况下，在95%的几个月内，我们预计该基金的利润或损失不超过1亿美元。换句话说，在任何给定的月份损失1亿美元或以上的概率为5%。
VaR的不足在于：
1）它没有描述可能最严重的损失。
2）VaR没有描述左尾的损失。它指示了值发生的概率，但未能描述左尾损失的分布。
3）VaR估计既受模型风险也受实施风险的影响。模型风险源于错误的假设，而实施风险是实施过程中出错的风险。
条件风险价值(CVaR)

条件风险价值（CVaR），也称为 expected shortfall (ES)，是由概率定义的损失期望值。换句话说，是预期损失投资回报率低于预先指定的最坏情况的回报率。例：基金5%的VaR是-25%，其含义为，在5%的时间里，该基金的回报率低于-25%。
简单来说 0.95 VaR衡量的是95% 里面的损失，而0.95CVaR衡量了96%-100%里面的损失，因此，CVaR的损失大于VaR。

图源https://analystprep.com/cfa-level-1-exam/portfolio-management/measuring-modifying-risks/

CVaR被认为是比VaR更好的风险衡量标准，因为与VaR不同，CVaR估计了不利事件的损失程度。我们构造CVaR不确定性集如下：

$\mathcal{Z}:=\left\{z \in \mathbb{R}^{n} \mid \exists \theta \in \mathbb{R}^{K}, z=\sum_{k=1}^{K} \theta_{k} \bar{z}_{k}, \theta \geq 0, \sum_{k=1}^{K} \theta_{k}=1, \theta \leq \frac{1}{K \alpha}\right\}$
注意其中 $z=\sum_{i=1}^K{\theta_k}\tilde{z_k}$ 是通过不同场景的安全平均求解z的值； ${\theta}\geq0$ 表示每个场景的出现概率大于等于0； $\sum_{i=1}^K{\theta_k = 1}$ 表示所有场景的概率和等于1； ${\theta}\leq\frac{1}{K_{\alpha}}$ 表示每个出现概率小于等于一个上界。是不确定性集大小的参数，可以简单地把 $\bar{z}_{k}$ 做为所有观察到的实现的集合[1]。如果我们想要找到最大的  ，就是需要解决下面这个线性规划问题：
$\begin{aligned} \underset{\theta, \gamma}{\operatorname{min}} &\quad \gamma \\ \text { s.t. } & z=\sum_{i=1}^{K} \theta_{k} \bar{z}_{k} \\ & 0 \leq \theta \leq \gamma / K \\ & \sum_{k=1}^{K} \theta_{k}=1 \end{aligned}$
并令 $\alpha=1/\gamma$ 。
三、鲁棒优化的模型的reformulation

原问题无法用solver直接求解因此需要进行reformulation。
step1: 将原问题写成epigraph form
$\begin{aligned} \max _{x,t}&\quad t \\ \text { s.t. } &t \leq(\mu+Pz)^Tx&&（1）\\ &0 \leq x_i \leq 1 \quad\forall i \in\{1, \ldots, n\}&&（2）\\ &\sum_{i=1}^{n} x_{i}=100 \%&&（3） \end{aligned}$

step2: 对于step1中的约束（1），求证 $\forall z \in Z,-P z^{T} x\leq-t+\mu^{T} x$ ，即求解以下问题： $\begin{aligned} \max _{z} &\quad-P z^{T} x \\ \text { s.t. } &z=\sum_{k=1}^{K} \theta_{k} \bar{z}_{k} \\ &\sum_{k=1}^{K} \theta_{k}=1 \\ &\theta_{k} \leq \frac{1}{K \alpha} \quad \forall k \in\{1, \ldots, K\} \\ &\theta_{k} \geq 0 \quad \forall k \in\{1, \ldots, K\} \end{aligned}$
将 $z=\sum_{k=1}^{K} \theta_{k} \bar{z}_{k}$ 代入
$\begin{aligned} \max _{\theta} &\quad-Px^{T}\left(\sum_{k=1}^{K} \theta_{k} \bar{z}_{k}\right) \\ \text { s.t. }& \sum_{k=1}^{K} \theta_{k}=1 \\ &\theta_{k} \leq \frac{1}{K \alpha} \quad \forall k \in\{1, \ldots, n\} \\ &\theta_{k} \geq 0 \quad \forall k \in\{1, \ldots, n\} \end{aligned}$

step3: 对step2中的问题求对偶：
$\begin{aligned} \min _{\lambda, \beta} &\quad\lambda+\frac{1}{K \alpha} \sum_{k=1}^{K} \beta_{k} \\ \text { s.t. }& \lambda+\beta_{k} \geq -Px^{T} \bar{z}_{k} \quad \forall k \in\{1, \ldots, K\} \\ &\beta_{k} \geq 0 \quad \forall k \in\{1, \ldots, K\} \end{aligned}$

step4: 原问题的reformulation为：
$\begin{aligned} \max _{x,t, \lambda, \beta}&\quad t \\ \text { s.t. } &t-\mu^Tx+\lambda+\frac{1}{K \alpha} \sum_{k=1}^{K} \beta_{k} \leq 0 \\ &\lambda+\beta_{k} \geq -Px^{T} \bar{z}_{k}\quad \forall k \in\{1, \ldots, K\}\\ &\beta_{k} \geq 0 \quad \forall i \in\{1, \ldots, K\}\\ &0 \leq x_i \leq 1 \quad\forall i \in\{1, \ldots, n\}\\ &\sum_{i=1}^{n} x_{i}=100 \% \end{aligned}$
至此原问题就转化为了一个可直接用solver求解的线性规划问题。
四、基于rome的鲁棒优化模型实现

rome包简介

rome是Joel Goh和Melvyn Sim开发的Matlab库，以促进鲁棒优化在线性规划建模的决策问题上的应用，常用于鲁棒模型的验证。该库在[2]中首次引入，可以用于以代数形式描述线性规划模型的细节，可在网站： http://www.robustopt.com 上下载并查阅相关文档。
代码实现

接下来实现一个股票投资问题的CVaR不确定性集构造，并利用rome分别求解原模型和reformulation后的模型
VaRProb = 0.05;
n = 10; %number of stocks, 10
k = 120; %number of months, 120
%Generate parameters randomly
temp=rand(n,1);
returns=zeros(n,k);
for i=1:n
returns(i,:)=[rand(1,k)-repelem(0.5,k)+temp(i,1)];
end
Zs = returns-mu*ones(1,k); %
rhat = max(abs(Zs),[],2); % a column vector containing the maximum derivation of each row
Zs = diag(1./rhat)*Zs; % each row of Zs divide its maximum
P = diag(rhat);

% Calibrating the CVaR uncertainty set
h = rome_begin;
newvar gamma_CvaR;
newvar Theta(m,m);
for j=1:k       %for Zs_j, use the jth column of Theta to do the combination
rome_constraint(Zs(:,j) == Zs*Theta(:,j));
rome_constraint(sum(Theta(:,j))==1);
end
rome_constraint(Theta <= gamma_CvaR/k);
rome_constraint(Theta>=0);
rome_minimize(gamma_CvaR);
h.solve;
obj_gamma_CvaR = h.objective;
Theta_Value = h.eval(Theta);
Theta_colum_max = sort(max(Theta_Value, [], 1)); % a row vector containing the maximum value of each column
Theta_CVaR =Theta_colum_max(k+1-ceil((1-VaRProb)*end));
alpha = 1/(Theta_CVaR*k);
rome_end;

% CvaR
disp(&#39;Solving the Origin model&#39;);
h = rome_begin;
newvar x(n) nonneg;
newvar y(1) ;
newvar z(n) uncertain;
newvar Theta(m) uncertain;
rome_constraint(z==returns*Theta);
rome_constraint(Theta>=0);
rome_constraint(sum(Theta)==1);
rome_constraint(Theta<=1/(k*alpha));
rome_maximize((mu+P*z)&#39;*x);
rome_constraint(x<=1);
rome_constraint(sum(x)==1);
h.solve;
obj_CVaR = h.objective;
xx_CVaR = h.eval(x); %get optimal portfolio
rome_end;

string=[&#39;worstObjCvaR=&#39;,num2str(obj_CVaR)];
%print answer
disp(string);

%reformulation式
disp(&#39;Solving the Reformulation&#39;);
h = rome_begin;
newvar x(n) nonneg;
newvar lambda;
newvar y;
newvar beta(m) nonneg;
rome_maximize(y);
rome_constraint(-mu&#39;*x+lambda+1/(k*alpha)*sum(beta)+y<=0);
for i=1:k
rome_constraint(lambda+beta(i)>=-((P*returns(:,i))&#39;*x));
end
rome_constraint(x<=1);
rome_constraint(sum(x)==1);
h.solve;
obj_ref_CVaR = h.objective;
xx_ref_CVaR = h.eval(x);
rome_end;

string2=[&#39;worstObjCvaR2=&#39;,num2str(obj_ref_CVaR)];
%print answer
disp(string2);求解结果示例：

求解结果

注意：由于案例数据随机生成，每次求解的目标值数值可能不一样，但是每次原模型和reformulation之后的模型结果应当一致。
参考文献

[1] J. Goh and M. Sim. Robust optimization made easy with rome. Operations Research, 59(4):973–985, 2011.
[2] Delage 2019 LectureNotes_v8 Robust Optimization.
[3] https://analystprep.com/cfa-level-1-exam/portfolio-management/measuring-modifying-risks/

<hr/>往期推文合集

第85篇：非线性优化 | 非线性问题线性化以及求解的详细案例及Python+Gurobi求解
第84篇：ORers&#39; Bling Chat | 【高光聊天记录集锦-02】：运小筹读者群里那些热烈的讨论
第83篇：Machine-Learning–Based Column Selection for Column Generation
第82篇：最新！205页运小筹优化理论学习笔记发布(2021年9月--2022年3月)！
第81篇：【鲁棒优化】| 补充证明:为什么最优解一定满足取值等于绝对值(论文笔记：The Price of Robustness)
第80篇：【鲁棒优化】| 论文笔记：The Price of Robustness - 列不确定性模型部分的推导和一些思考
第79篇：ORers&#39; Bling Chat | 【高光聊天记录集锦-01】：运小筹读者群里那些热烈的讨论
第78篇：优化| 手把手教你学会杉数求解器(COPT)的安装、配置与测试
第77篇：【教学视频】优化 | 线性化(2)：连续变量 * 0-1变量的线性化
第76篇：【教学视频】优化 | 线性化：两个0-1变量相乘的线性化
第75篇：强化学习实战 | DQN和Double DQN保姆级教程：以Cart-Pole为例
第74篇：强化学习| 概念梳理：强化学习、马尔科夫决策过程与动态规划
第73篇：强化学习实战 | Q-learning求解最短路(SPP)问题
第72篇：鲁棒优化 | 以Coding入门鲁棒优化：以一个例子引入（二）
第71篇：鲁棒优化|基于ROME编程入门鲁棒优化：以一个例子引入(一)
第70篇：优化|含分式的非线性规划求解: 基于Outer Approximation的Branch-and-cut 算法及其Java实现
第69篇：科研工具 | 手把手教你玩转文献管理神器：Endnote
第68篇：相约深圳 | 2022 INFORMS 服务科学国际会议·征稿通知
第67篇：鲁棒优化｜利用rome求解鲁棒优化简单案例：以CVaR不确定性集为例
第66篇：机器学习 | 交通流特征工程小技巧和思考
第65篇：最新！145页运小筹优化理论学习笔记发布(2021年4月--9月)！
第64篇：优化 | 手把手教你用Python调用SCIP求解最优化模型
第63篇：优化 | 随机规划案例：The value of the stochastic solution
第62篇：工具 | draw.io: 科研流程示意图必备大杀器
第61篇：优化 | 开源求解器SCIP的安装与入门教程（Windows+Python）
第60篇：优化|Gurobi处理非线性目标函数和约束的详细案例
第59篇：让出租车更“懂”您的出行
第58篇：测试算例下载：《运筹优化常用模型、算法及案例实战：代码手册》
第57篇：鲁棒优化|分布式鲁棒优化转化为SOCP案例及代码实现(Python+RSOME)
第56篇：鲁棒优化 | 分布式鲁棒优化简单案例及实战(RSOME+Gurobi)
第55篇：【尾款+发货】|【代码手册】《运筹优化常用模型、算法及案例实战：Python+Java
第54篇：深度强化学习之：PPO训练红白机1942
第53篇：简单装配线平衡问题
第52篇：【视频讲解】CPLEX的OPL语言：可视化的优化模型求解解决方案
第51篇：算法 | 基于英雄联盟寻路背景的A星算法及python实现
第50篇：【转发】清华大学深圳国际研究生院2021年物流工程与管理项目优秀大学生夏令营报名通知
第49篇：优化 | 精确算法之分支定界介绍和实现（附Python代码）
第48篇：【顶刊论文速递】综述：服务科学视角下的金融科技
第47篇：【重新发布】|《运筹优化常用模型、算法及案例实战：Python+Java实现》【代码手册】开始预购啦！！！
第46篇：智慧交通可视化：手把手教你用Python做出行数据可视化-osmnx 包入门教程
第45篇：优化 | Pick and delivery problem的介绍与建模实现（二）
第44篇：推一个知乎学弱猹的公众号
第43篇：元启发式算法 | 遗传算法(GA)解决TSP问题(Python实现)
第42篇：优化|视频详解Python调用Gurobi的应用案例：TSP
第41篇：最新！213页运小筹优化理论系列笔记发布！
第40篇：运小筹第四期月刊发布！
第39篇：开源交通仿真平台SimMobility的安装教程
第38篇：浅谈 | P问题、NP问题、NPC问题、NP-Hard问题
第37篇：一份掏心掏肺的信息素养笔记分享
第36篇：强化学习｜Q-learning (王者荣耀视角)
第35篇：优化｜高级建模方法（Gurobi）：线性化表达小技巧
第34篇：深度强化学习介绍 | Deep Q Network——理论与代码实现
第33篇：优化 | 列生成算法及Java调用cplex实现
第32篇：优化 | Pick and delivery problem的简介与建模实现（一）
第31篇：最新！运小筹月刊-1月份发布！
第30篇：“Learn to Improve”（L2I）：ICLR文章分享 | 运用RL解决VRP问题
第29篇：线性规划求解小程序 | 基于Python 的Cplex 求解器图形化界面
第28篇：运筹学与管理科学TOP期刊揭秘 —TR系列
第27篇：Julia安装与配置Jupyter Notebook
第26篇：期刊征文| IEEE TRANSACTIONS应对COVID-19的特刊征文消息速递
第25篇：两阶段随机规划(Two-stage Stochastic Programming)：一个详细的例子
第24篇：最新！运小筹月刊-12月份发布！
第23篇：Python调用Gurobi：Assignment Problem简单案例
第22篇：初识随机规划：一个小小例子
第21篇：机器学习运用到VRP的若干小知识
第20篇：运筹学与管理科学TOP期刊揭秘 —Service Science
第19篇：手把手教你用Python实现Dijkstra算法（伪代码+Python实现）
第18篇：运筹学与管理科学大揭秘—TOP期刊主编及研究方向一览
第17篇：优化 | 手把手教你用Python实现动态规划Labeling算法求解SPPRC问题
第16篇：代码 | 运小筹GitHub项目上线啦，快来标星收藏吧！！
第15篇：最新！运小筹月刊首次发布！
第14篇：优化| 手把手教你用Python实现Dijkstra算法（附Python代码和可视化）
第13篇：优化|手把手教你用Python调用Gurobi求解最短路问题
第12篇：优化| 手把手教你用Java实现单纯形法
第11篇：优化| 手把手教你用Python调用Gurobi求解VRPTW
第10篇：优化 | 两阶段启发式算法求解带时间窗车辆路径问题（附Java代码）
第9篇：Java调用cplex求解运输问题
第8篇：优化 | 如何优雅地写出大规模线性规划的对偶
第7篇：优化 | TSP中两种不同消除子环路的方法及Callback实现（Python调用Gurobi求解）
第6篇：元启发式算法 | 禁忌搜索(Tabu Search)解决TSP问题(Python代码实现)
第5篇：论文代码复现 | 无人机与卡车联合配送(Python+Gurobi)
第4篇：优化|Shortest Path Problem及其对偶问题的一些探讨(附Python调用Gurobi实现)
第3篇：可视化|Python+OpenStreetMap的交通数据可视化(二)：OpenStreetMap+Python画城市交通路网图
第2篇：优化|最速下降法：详细案例+Python实现
第1篇：Python+networkX+OpenStreetMap实现交通数据可视化(一)：用OpenStreetMap下载地图数据
<hr/>运小筹公众号致力于分享运筹优化(LP、MIP、NLP、随机规划、鲁棒优化)、凸优化、强化学习等研究领域的内容以及涉及到的算法的代码实现。编程语言和工具包括Java、Python、Matlab、CPLEX、Gurobi、SCIP 等。欢迎广大同行投稿！欢迎大家关注我们的公众号“运小筹”以及粉丝群！
可以私信编辑 MunSum3 拉入读者群
<hr/>作者：杨影，清华大学，深圳国际研究生院物流工程与管理
编辑：张瑞三，四川大学，硕士在读，邮箱：zrssum3@stu.scu.edu.cn、sum3rebirth@gmail.com

		自动登录	找回密码
密码			立即注册

鲁棒优化｜ 利用rome求解鲁棒优化简单案例

本帖子中包含更多资源

鲁棒优化｜利用rome求解鲁棒优化简单案例