鲁棒优化是处理数学规划中处理不确定性的重要方法论。这里我们以几个简单的例子来介绍鲁棒优化。同时，我们结合鲁棒优化工具包Robust Optimization Made Easy (ROME)来编程实现，方便读者理解和复现。

线性规划中的不确定性

在生产实际中，一些生产计划、运输问题可以建模为线性规划进行求解，从而达到利润的最大化或者成本的最小化。但是这些建模是基于一个重要假设：

• 生产成本/生产要素相关的参数是确定的。比如药物生产中原材料的含量，这是基于检测得到的。

这些线性规划模型，一般可以表示为下面的形式：

$$\begin{array}{rl}\max & c^{T}x\\ s.t.& Ax⩽b\\ & x⩾0\end{array}$$
其中，$c$为$n\times 1$的列向量，$b$为$m\times 1$的列向量，$A$为$m\times n$的矩阵，$x$为决策变量，纬度为$n\times 1$。$c$为价值系数，$b$为资源向量，$A$为生产要素矩阵。

目标函数$\max c^{T}x$ 表示要最大化收益，约束$Ax⩽b$ 表示生产所用的原材料要小于等于可用的生产资源。

但是在实际生产生活中，价值系数$c$，资源向量$b$，生产要素矩阵$A$都有可能是不确定的，也就是具有不确定性。例如，股票投资中，股票的预期收益是不确定的，无法提前预知的。再比如实际生产中，生产实际消耗的资源(对应$A$)也是不确定的，会由于施工精度，施工条件和工人的操作产生偏差。

如果我们在建模中考虑这些不确定的因素，则价值系数$c$，资源向量$b$和生产要素矩阵$A$则可以分别表示为$\tilde{c},\tilde{b},\tilde{A}$. 因此，考虑不确定因素的线性规划模型可以写为

$$\begin{array}{rl}\max & {\tilde{c}}^{T}x\\ s.t.& \tilde{A}x⩽\tilde{b}\\ & x⩾0\end{array}$$

当然啦，实际中，上述不确定也许不会同时出现。或者，决策者不想同时考虑。因此有时决策者只会考虑其中一种(例如只考虑目标函数系数$c$不确定)，或者多种。

处理这种带有不确定性的线性规划模型，一般有两大类方法论：

1. 鲁棒优化(Robust Optimization): 这种优化是考虑在不确定性的影响下，追求最优化最坏的情形下(worst case)的解。即：求大佬把我最点背时候的解帮我优化到最好，好让我最点背的时候也能尽量损失少。这种方法论一般比较保守，因为现实生活中，最点背(worst case)的情况发生的概率比较小，为了避免worst case而完全忽略眼前的利益，确实有点太保守。这也正像曾国藩的话，未来不迎，考虑未来太过度了也是不好的。
2. 随机规划(Stochastic Programming): 随机规划是倾向于追求在不确定性下，最大化目标函数的期望，这种方法论会综合考虑可能发生的场景，将不确定性以扩大问题规模的方式，纳入模型，最终得到一个可以一定程度上抵御不确定性的解。随机规划的解相较鲁棒优化而言，保守性显著降低，但是一个主要的缺点就是：模型规模被显著扩大，给求解带来了极大的挑战。关于随机规划，我们已经有几期推文进行了详细的介绍，并附有完整代码，感兴趣的读者们请关注《运小筹》公众号查看，也可以点击下方链接查看学习。

我们这期推文主要来讲解鲁棒优化。

一个生产的例子

本例子来自文献[1]。

某公司生产两种包含特定成分A 的药物：分别为药物1 和药物2。成分A 是从市场上购买的原材料提取而来。相关的生产成本以及原料数据见下表。我们的目标是指定一个生产计划，使得公司能够获得的利润最大。

我们将其建模为一个线性规划模型，如下

$$\begin{array}{rl}\underset{x,y}{\max }& 6200y_1+6900y_2-\left(100x_1+199.90x_2+700y_1+800y_2\right)\\ & x_1+x_2⩽1000\text{(Storage)}\\ & 90y_1+100y_2⩽2000\text{(Manpower)}\\ & 40y_1+50y_2⩽800\text{(Equipment)}\\ & 100x_1+199.90x_2+700y_1+800y_2⩽100000\text{(Budget)}\\ & 0.01x_1+0.02x_2-0.5y_1-0.6y_2⩾0\text{(Agent A)}\\ & x_1,x_2,y_1,y_2⩾0\end{array}$$

Matlab+ROME求解

网址：https://robustopt.com/

ROME包简介:
ROME全称：Robust Optimization Made Easy(ROME) (调皮一下，我给这个包来个翻译：轻松拿下鲁棒优化)，是两位鲁棒优化大佬Joel Goh和Melvyn Sim团队开发的基于Matlab的开源鲁棒优化工具包，可以实现不进行reformulation就直接求解鲁棒优化模型的功能，也可以用于求解一般的优化模型。改开源工具包目前已经有了python版本，是Zhi Chen，Melvyn Sim，Peng Xiong联合开发的，叫做Robust Stochastic Optimization Made Easy, RSOME (见熊鹏老师的github)(同样给个翻译：轻松拿下鲁棒随机优化)， 见网址https://xiongpengnus.github.io/rsome/examples。这些软件也都伴随着Operations Research和Management Science论文的发表，论文如下：

• ROME: Joel Goh and Melvyn Sim. Robust Optimization Made Easy with ROME. Operations Research, 2011, 59(4), pp.973-985.
• Joel Goh and Melvyn Sim. Distributionally Robust Optimization and its Tractable Approximations. Operations Research, 2010, 58(4), pp. 902-917.

以及RSOME对应的论文

• Chen, Zhi, and Peng Xiong. 2021. RSOME in Python: an open-source package for robust stochastic optimization made easy. Optimization Online.
• Chen, Zhi, Melvyn Sim, Peng Xiong. 2020. Robust stochastic optimization made easy with RSOME. Management Science 66(8) 3329–3339.

在这里向前辈大佬们致以崇高的敬意！！！

我们用Matlab调用ROME求解上述模型，代码如下(代码部分参考自官网代码和参考文献)：

clc;
clear;
% begin rome
h = rome_begin;
an1 = 0.01;
an2 = 0.02;
% decision variables
newvar Drug1 nonneg;     % nonneg:非负
newvar Drug2 nonneg;
newvar Raw1 nonneg;
newvar Raw2 nonneg;% Objective to maximize the profit
rome_maximize(6200 * Drug1 + 6900*Drug2 - (100 * Raw1 + 199 * Raw2 + 700 * Drug1 + 800 * Drug2));% Storage constraint
rome_constraint(Raw1 + Raw2 <= 1000);% manpower constraint
rome_constraint(90 * Drug1 + 100 * Drug2 <= 2000);% Equipment constraint
rome_constraint(40 * Drug1 + 50 * Drug2 <= 800);%Budget constraint
rome_constraint(100 * Raw1 + 199.9 * Raw2 + 700 * Drug1 + 800 * Drug2 <= 100000);%Constraint to have enough active agent A
rome_constraint(an1 * Raw1 + an2 * Raw2 - 0.5 * Drug1 - 0.6 * Drug2 >= 0);%solve the model
h.solve;
optobj_det = h.objective;               % get optimal value
det_sol.Drug1 = h.eval(Drug1);      % get optimal solution
det_sol.Drug2 = h.eval(Drug2);     % get optimal solution
det_sol.Raw1 = h.eval(Raw1);        % get optimal solution
det_sol.Raw2 = h.eval(Raw2);       % get optimal solutionfprintf('Profit:\t  %10.3f\n', optobj_det);
fprintf('Drug1:\t  %10.3f\n', det_sol.Drug1);
fprintf('Drug2:\t  %10.3f\n', det_sol.Drug2);
fprintf('Raw1:\t  %10.3f\n', det_sol.Raw1);
fprintf('Raw2:\t  %10.3f\n', det_sol.Raw2);
rome_end;

求解得到结果如下：

ROME is using CPLEX class because it is available
Status: integer optimal solution
Profit:	    8819.658
Drug1:	      17.552
Drug2:	       0.000
Raw1:	       0.000
Raw2:	     438.789

根据结果知，生产Drug 1的量为17552，不生产Drug 2，最终利润为8819.658。

但是由于原材料所含成分A的含量实际上是不确定的，上面模型中个给出每kg原材料1包含物质A的量为0.01g，每kg原材料2包含物质A的量为0.02g。如果这个数字由于一些测量和机器操作等误差导致偏差，比如实际每kg原材料1包含物质A的量减少了0.005g，仅为0.0095g等，那么生产计划和利润如何呢？

我们假设，A物质的含量变化为0.0095和0.0196 （分别减少5%和2%），则模型的参数变化为

an1 = 0.0095;
an2 = 0.0196;

求解得到的结果变化为

Status: Optimal solution found.
Profit:	    6905.956
Drug1:	      17.243
Drug2:	       0.000
Raw1:	       0.000
Raw2:	     439.870

我们发现，利润只有8294.567，相比原来的8819.658减少了1913.7，居然减少了21.7%。原来能生产17552份Drug1，现在只能生产17243。这样的变化实在是太大了。这说明，在实际生产中，如果由于参数测量或者预估不准确，有可能会导致企业产生巨大的经济损失。

当然了，有些时候，还会导致无可行解，即生产任务无法完成。

通过上述例子，相信读者们都对优化中考虑参数不确定性的重要性有了直观的认识，接下来我们以一个具体的例子来阐述如何通过鲁棒优化来应对这种不确定性，提高解的抵抗不确定性(风险)的能力。

鲁棒优化简单案例：一个投资组合的例子

在上面的例子中，我们发现，当模型的一些参数与我们初始预估的参数发生轻微
的变化的时候，我们得到的解就会变得不可行，也就是不可能被实施。

下面我们来看一个股票投资组合的例。见文献Goh. 2011。本例说明了，除了约束以外，目标函数中有不确定性因素，也会影响到问题的解。

• M. Sim J. Goh. Robust optimization made easy software. 2011.

有$n$种股票，在确定性问题下，每种股票的投资回报率为$r_i$ ，我们需要决策每种股票的投资占比，从而使得收益最大。确定型模型如下。

$$\begin{array}{rlrl}\max & \sum _{i=1}^n{r}_i{x}_i& & \\ & \sum _{i=1}^n{x}_i=100\%& & \\ & x_i⩾0,& & \forall i=1,\cdots ,n\end{array}$$

但是投资回报率是有波动的，确定型问题中的回报率$r_i$只是一个期望值。

不确定回报率

为了建立robust optimization模型，我们考虑投资回报率$r_i$是不确定的变量，即

$$\begin{array}{r}{\tilde{r}}_i={\mu }_i+{\sigma }_i{z}_i\end{array}$$其中${\mu }_i,{\sigma }_i$是已知量，是第$i$种股票投资回报率的均值和方差，$z_i\in \left[-1,1\right]$ 是具有不确定性的决策变量, 主要是为了控制${\tilde{r}}_i$在$\left[r_i-{\sigma }_i,r_i+{\sigma }_i\right]$之间变动。

我们考虑以下算例：有150种股票，$n=150$，每种股票的均值和方差都是随着股票的ID依次递增。并且
$$\begin{array}{rl}{\mu }_i& =0.15+i\frac{0.05}{150}\\ {\sigma }_i& =\frac{0.05}{450}\sqrt{2ni\left(n+1\right)}\end{array}$$
因此，可以得到，150种股票的投资回报率的${\mu }_i,{\sigma }_i$以及投资回报率${\tilde{r}}_i$的波动范围。

• 150种股票的投资回报率的${\mu }_i$

-----------------------------------------以下是每支股票的收益均值----------------------------
#  1 -- 0.1503      #  2 -- 0.1507      #  3 -- 0.1510      #  4 -- 0.1513      #  5 -- 0.1517      #  6 -- 0.1520      
#  7 -- 0.1523      #  8 -- 0.1527      #  9 -- 0.1530      # 10 -- 0.1533      # 11 -- 0.1537      # 12 -- 0.1540      
# 13 -- 0.1543      # 14 -- 0.1547      # 15 -- 0.1550      # 16 -- 0.1553      # 17 -- 0.1557      # 18 -- 0.1560      
# 19 -- 0.1563      # 20 -- 0.1567      # 21 -- 0.1570      # 22 -- 0.1573      # 23 -- 0.1577      # 24 -- 0.1580      
# 25 -- 0.1583      # 26 -- 0.1587      # 27 -- 0.1590      # 28 -- 0.1593      # 29 -- 0.1597      # 30 -- 0.1600      
# 31 -- 0.1603      # 32 -- 0.1607      # 33 -- 0.1610      # 34 -- 0.1613      # 35 -- 0.1617      # 36 -- 0.1620      
# 37 -- 0.1623      # 38 -- 0.1627      # 39 -- 0.1630      # 40 -- 0.1633      # 41 -- 0.1637      # 42 -- 0.1640      
# 43 -- 0.1643      # 44 -- 0.1647      # 45 -- 0.1650      # 46 -- 0.1653      # 47 -- 0.1657      # 48 -- 0.1660      
# 49 -- 0.1663      # 50 -- 0.1667      # 51 -- 0.1670      # 52 -- 0.1673      # 53 -- 0.1677      # 54 -- 0.1680      
# 55 -- 0.1683      # 56 -- 0.1687      # 57 -- 0.1690      # 58 -- 0.1693      # 59 -- 0.1697      # 60 -- 0.1700      
# 61 -- 0.1703      # 62 -- 0.1707      # 63 -- 0.1710      # 64 -- 0.1713      # 65 -- 0.1717      # 66 -- 0.1720      
# 67 -- 0.1723      # 68 -- 0.1727      # 69 -- 0.1730      # 70 -- 0.1733      # 71 -- 0.1737      # 72 -- 0.1740      
# 73 -- 0.1743      # 74 -- 0.1747      # 75 -- 0.1750      # 76 -- 0.1753      # 77 -- 0.1757      # 78 -- 0.1760      
# 79 -- 0.1763      # 80 -- 0.1767      # 81 -- 0.1770      # 82 -- 0.1773      # 83 -- 0.1777      # 84 -- 0.1780      
# 85 -- 0.1783      # 86 -- 0.1787      # 87 -- 0.1790      # 88 -- 0.1793      # 89 -- 0.1797      # 90 -- 0.1800      
# 91 -- 0.1803      # 92 -- 0.1807      # 93 -- 0.1810      # 94 -- 0.1813      # 95 -- 0.1817      # 96 -- 0.1820      
# 97 -- 0.1823      # 98 -- 0.1827      # 99 -- 0.1830      #100 -- 0.1833      #101 -- 0.1837      #102 -- 0.1840      
#103 -- 0.1843      #104 -- 0.1847      #105 -- 0.1850      #106 -- 0.1853      #107 -- 0.1857      #108 -- 0.1860      
#109 -- 0.1863      #110 -- 0.1867      #111 -- 0.1870      #112 -- 0.1873      #113 -- 0.1877      #114 -- 0.1880      
#115 -- 0.1883      #116 -- 0.1887      #117 -- 0.1890      #118 -- 0.1893      #119 -- 0.1897      #120 -- 0.1900      
#121 -- 0.1903      #122 -- 0.1907      #123 -- 0.1910      #124 -- 0.1913      #125 -- 0.1917      #126 -- 0.1920      
#127 -- 0.1923      #128 -- 0.1927      #129 -- 0.1930      #130 -- 0.1933      #131 -- 0.1937      #132 -- 0.1940      
#133 -- 0.1943      #134 -- 0.1947      #135 -- 0.1950      #136 -- 0.1953      #137 -- 0.1957      #138 -- 0.1960      
#139 -- 0.1963      #140 -- 0.1967      #141 -- 0.1970      #142 -- 0.1973      #143 -- 0.1977      #144 -- 0.1980      
#145 -- 0.1983      #146 -- 0.1987      #147 -- 0.1990      #148 -- 0.1993      #149 -- 0.1997      #150 -- 0.2000 

• 150种股票的投资回报率的${\sigma }_i$

-----------------------------------------以下是每支股票的收益标准差----------------------------
#  1 -- 0.0236      #  2 -- 0.0334      #  3 -- 0.0410      #  4 -- 0.0473      #  5 -- 0.0529      #  6 -- 0.0579      
#  7 -- 0.0626      #  8 -- 0.0669      #  9 -- 0.0709      # 10 -- 0.0748      # 11 -- 0.0784      # 12 -- 0.0819      
# 13 -- 0.0853      # 14 -- 0.0885      # 15 -- 0.0916      # 16 -- 0.0946      # 17 -- 0.0975      # 18 -- 0.1003      
# 19 -- 0.1031      # 20 -- 0.1058      # 21 -- 0.1084      # 22 -- 0.1109      # 23 -- 0.1134      # 24 -- 0.1159      
# 25 -- 0.1182      # 26 -- 0.1206      # 27 -- 0.1229      # 28 -- 0.1251      # 29 -- 0.1274      # 30 -- 0.1295      
# 31 -- 0.1317      # 32 -- 0.1338      # 33 -- 0.1359      # 34 -- 0.1379      # 35 -- 0.1399      # 36 -- 0.1419      
# 37 -- 0.1438      # 38 -- 0.1458      # 39 -- 0.1477      # 40 -- 0.1496      # 41 -- 0.1514      # 42 -- 0.1533      
# 43 -- 0.1551      # 44 -- 0.1569      # 45 -- 0.1586      # 46 -- 0.1604      # 47 -- 0.1621      # 48 -- 0.1638      
# 49 -- 0.1655      # 50 -- 0.1672      # 51 -- 0.1689      # 52 -- 0.1705      # 53 -- 0.1722      # 54 -- 0.1738      
# 55 -- 0.1754      # 56 -- 0.1770      # 57 -- 0.1785      # 58 -- 0.1801      # 59 -- 0.1816      # 60 -- 0.1832      
# 61 -- 0.1847      # 62 -- 0.1862      # 63 -- 0.1877      # 64 -- 0.1892      # 65 -- 0.1907      # 66 -- 0.1921      
# 67 -- 0.1936      # 68 -- 0.1950      # 69 -- 0.1964      # 70 -- 0.1979      # 71 -- 0.1993      # 72 -- 0.2007      
# 73 -- 0.2021      # 74 -- 0.2034      # 75 -- 0.2048      # 76 -- 0.2062      # 77 -- 0.2075      # 78 -- 0.2089      
# 79 -- 0.2102      # 80 -- 0.2115      # 81 -- 0.2128      # 82 -- 0.2141      # 83 -- 0.2154      # 84 -- 0.2167      
# 85 -- 0.2180      # 86 -- 0.2193      # 87 -- 0.2206      # 88 -- 0.2218      # 89 -- 0.2231      # 90 -- 0.2244      
# 91 -- 0.2256      # 92 -- 0.2268      # 93 -- 0.2281      # 94 -- 0.2293      # 95 -- 0.2305      # 96 -- 0.2317      
# 97 -- 0.2329      # 98 -- 0.2341      # 99 -- 0.2353      #100 -- 0.2365      #101 -- 0.2377      #102 -- 0.2388      
#103 -- 0.2400      #104 -- 0.2412      #105 -- 0.2423      #106 -- 0.2435      #107 -- 0.2446      #108 -- 0.2458      
#109 -- 0.2469      #110 -- 0.2480      #111 -- 0.2492      #112 -- 0.2503      #113 -- 0.2514      #114 -- 0.2525      
#115 -- 0.2536      #116 -- 0.2547      #117 -- 0.2558      #118 -- 0.2569      #119 -- 0.2580      #120 -- 0.2591      
#121 -- 0.2601      #122 -- 0.2612      #123 -- 0.2623      #124 -- 0.2633      #125 -- 0.2644      #126 -- 0.2655      
#127 -- 0.2665      #128 -- 0.2676      #129 -- 0.2686      #130 -- 0.2696      #131 -- 0.2707      #132 -- 0.2717      
#133 -- 0.2727      #134 -- 0.2738      #135 -- 0.2748      #136 -- 0.2758      #137 -- 0.2768      #138 -- 0.2778      
#139 -- 0.2788      #140 -- 0.2798      #141 -- 0.2808      #142 -- 0.2818      #143 -- 0.2828      #144 -- 0.2838      
#145 -- 0.2848      #146 -- 0.2857      #147 -- 0.2867      #148 -- 0.2877      #149 -- 0.2887      #150 -- 0.2896      

• 150种股票的投资回报率的波动范围如下：
  $$\begin{array}{r}{\tilde{r}}_i\in \left[0.15+i\frac{0.05}{150}-\frac{0.05}{450}\sqrt{2ni\left(n+1\right)},0.15+i\frac{0.05}{150}+\frac{0.05}{450}\sqrt{2ni\left(n+1\right)}\right]\end{array}$$

-----------------------------------------以下是r的变动范围----------------------------
#  1 --[ 0.1267,  0.1740]     #  2 --[ 0.1172,  0.1841]     #  3 --[ 0.1100,  0.1920]     #  4 --[ 0.1040,  0.1986]     
#  5 --[ 0.0988,  0.2045]     #  6 --[ 0.0941,  0.2099]     #  7 --[ 0.0898,  0.2149]     #  8 --[ 0.0858,  0.2196]     
#  9 --[ 0.0821,  0.2239]     # 10 --[ 0.0785,  0.2281]     # 11 --[ 0.0752,  0.2321]     # 12 --[ 0.0721,  0.2359]     
# 13 --[ 0.0691,  0.2396]     # 14 --[ 0.0662,  0.2432]     # 15 --[ 0.0634,  0.2466]     # 16 --[ 0.0607,  0.2499]     
# 17 --[ 0.0582,  0.2532]     # 18 --[ 0.0557,  0.2563]     # 19 --[ 0.0533,  0.2594]     # 20 --[ 0.0509,  0.2624]     
# 21 --[ 0.0486,  0.2654]     # 22 --[ 0.0464,  0.2683]     # 23 --[ 0.0443,  0.2711]     # 24 --[ 0.0421,  0.2739]     
# 25 --[ 0.0401,  0.2766]     # 26 --[ 0.0381,  0.2793]     # 27 --[ 0.0361,  0.2819]     # 28 --[ 0.0342,  0.2845]     
# 29 --[ 0.0323,  0.2870]     # 30 --[ 0.0305,  0.2895]     # 31 --[ 0.0287,  0.2920]     # 32 --[ 0.0269,  0.2944]     
# 33 --[ 0.0251,  0.2969]     # 34 --[ 0.0234,  0.2992]     # 35 --[ 0.0218,  0.3016]     # 36 --[ 0.0201,  0.3039]     
# 37 --[ 0.0185,  0.3062]     # 38 --[ 0.0169,  0.3084]     # 39 --[ 0.0153,  0.3107]     # 40 --[ 0.0138,  0.3129]     
# 41 --[ 0.0122,  0.3151]     # 42 --[ 0.0107,  0.3173]     # 43 --[ 0.0093,  0.3194]     # 44 --[ 0.0078,  0.3215]     
# 45 --[ 0.0064,  0.3236]     # 46 --[ 0.0049,  0.3257]     # 47 --[ 0.0035,  0.3278]     # 48 --[ 0.0022,  0.3298]     
# 49 --[ 0.0008,  0.3319]     # 50 --[-0.0006,  0.3339]     # 51 --[-0.0019,  0.3359]     # 52 --[-0.0032,  0.3379]     
# 53 --[-0.0045,  0.3398]     # 54 --[-0.0058,  0.3418]     # 55 --[-0.0070,  0.3437]     # 56 --[-0.0083,  0.3456]     
# 57 --[-0.0095,  0.3475]     # 58 --[-0.0108,  0.3494]     # 59 --[-0.0120,  0.3513]     # 60 --[-0.0132,  0.3532]     
# 61 --[-0.0144,  0.3550]     # 62 --[-0.0155,  0.3569]     # 63 --[-0.0167,  0.3587]     # 64 --[-0.0179,  0.3605]     
# 65 --[-0.0190,  0.3623]     # 66 --[-0.0201,  0.3641]     # 67 --[-0.0212,  0.3659]     # 68 --[-0.0223,  0.3677]     
# 69 --[-0.0234,  0.3694]     # 70 --[-0.0245,  0.3712]     # 71 --[-0.0256,  0.3729]     # 72 --[-0.0267,  0.3747]     
# 73 --[-0.0277,  0.3764]     # 74 --[-0.0288,  0.3781]     # 75 --[-0.0298,  0.3798]     # 76 --[-0.0308,  0.3815]     
# 77 --[-0.0318,  0.3832]     # 78 --[-0.0329,  0.3849]     # 79 --[-0.0339,  0.3865]     # 80 --[-0.0349,  0.3882]     
# 81 --[-0.0358,  0.3898]     # 82 --[-0.0368,  0.3915]     # 83 --[-0.0378,  0.3931]     # 84 --[-0.0387,  0.3947]     
# 85 --[-0.0397,  0.3964]     # 86 --[-0.0406,  0.3980]     # 87 --[-0.0416,  0.3996]     # 88 --[-0.0425,  0.4012]     
# 89 --[-0.0434,  0.4028]     # 90 --[-0.0444,  0.4044]     # 91 --[-0.0453,  0.4059]     # 92 --[-0.0462,  0.4075]     
# 93 --[-0.0471,  0.4091]     # 94 --[-0.0479,  0.4106]     # 95 --[-0.0488,  0.4122]     # 96 --[-0.0497,  0.4137]     
# 97 --[-0.0506,  0.4152]     # 98 --[-0.0514,  0.4168]     # 99 --[-0.0523,  0.4183]     #100 --[-0.0532,  0.4198]     
#101 --[-0.0540,  0.4213]     #102 --[-0.0548,  0.4228]     #103 --[-0.0557,  0.4243]     #104 --[-0.0565,  0.4258]     
#105 --[-0.0573,  0.4273]     #106 --[-0.0581,  0.4288]     #107 --[-0.0590,  0.4303]     #108 --[-0.0598,  0.4318]     
#109 --[-0.0606,  0.4332]     #110 --[-0.0614,  0.4347]     #111 --[-0.0622,  0.4362]     #112 --[-0.0629,  0.4376]     
#113 --[-0.0637,  0.4391]     #114 --[-0.0645,  0.4405]     #115 --[-0.0653,  0.4419]     #116 --[-0.0660,  0.4434]     
#117 --[-0.0668,  0.4448]     #118 --[-0.0676,  0.4462]     #119 --[-0.0683,  0.4476]     #120 --[-0.0691,  0.4491]     
#121 --[-0.0698,  0.4505]     #122 --[-0.0705,  0.4519]     #123 --[-0.0713,  0.4533]     #124 --[-0.0720,  0.4547]     
#125 --[-0.0727,  0.4561]     #126 --[-0.0735,  0.4575]     #127 --[-0.0742,  0.4588]     #128 --[-0.0749,  0.4602]     
#129 --[-0.0756,  0.4616]     #130 --[-0.0763,  0.4630]     #131 --[-0.0770,  0.4643]     #132 --[-0.0777,  0.4657]     
#133 --[-0.0784,  0.4671]     #134 --[-0.0791,  0.4684]     #135 --[-0.0798,  0.4698]     #136 --[-0.0805,  0.4711]     
#137 --[-0.0811,  0.4725]     #138 --[-0.0818,  0.4738]     #139 --[-0.0825,  0.4751]     #140 --[-0.0831,  0.4765]     
#141 --[-0.0838,  0.4778]     #142 --[-0.0845,  0.4791]     #143 --[-0.0851,  0.4805]     #144 --[-0.0858,  0.4818]     
#145 --[-0.0864,  0.4831]     #146 --[-0.0871,  0.4844]     #147 --[-0.0877,  0.4857]     #148 --[-0.0884,  0.4870]     
#149 --[-0.0890,  0.4883]     #150 --[-0.0896,  0.4896]     

鲁棒优化模型

鲁棒优化模型，是要考虑在一定的风险范围之内，要使得最坏的情况(worst case)的解最好。这个在一定的风险范围之内如何来衡量，就是鲁棒优化研究中的一个重点，也就是大家熟知的不确定集（uncertainty set）（在分布鲁棒优化中，这个相应的变成了模糊集ambiguity set）。

描绘不确定集的方法有多种，在基本的鲁棒优化中，常见的有：

1. 预算不确定集(Budgeted uncertainty set)
2. 凸包不确定集(Convex hull uncertainty set)
3. 盒子不确定集(Box uncertainty set)
4. CVaR不确定集(CVaR uncertainty Set)
5. 椭球不确定集(Ellipsoidal uncertainty Set)

当然啦，除了用不确定集来处理，还可以用机会约束(Chance Constraints)来处理不确定性，这些就等日后再来写相关的科普文章。这里我们仅用简单的例子来入门。

另外，我们也这里也放出鲁棒优化的经典文献：

• Bertsimas D, Sim M. The price of robustness[J]. Operations research, 2004, 52(1): 35-53.

本文中，我们以预算不确定集(Budgeted uncertainty set)为例，来对上述问题进行建模。我们设置：至多有$\Gamma$个股票的投资回报率会偏离其均值$\mu$。为了方便叙述，取$\Gamma =4$。这里的$\Gamma$就是预算，通俗的理解，就是我们考虑明天股票走势最坏的情况，150支股票的收益相较平均值的波动的比例的总和，不会高于400%。当然了，如果设置$\Gamma =150$，表示你考虑了最点背的情况，你假想你是全世界最点背的人，明天有可能150支都跌倒最低点，哈哈。反之，如果设置$\Gamma$较小，说明你是一个乐观主义者，你比较看好明天股票的走势。

基于此，我们给出鲁棒优化模型如下：

$$\begin{array}{rlrl}\underset{x\in {\mathbb{R}}^n}{\max }\underset{z}{\min }& \sum _{i=1}^n\left({\mu }_i+{\sigma }_i{z}_i\right)x_i& & \\ & \sum _{i=1}^n{x}_i=100\%& & \\ & \sum _{i=1}^n|z_i|⩽\Gamma & & \\ & x_i⩾0,& & \forall i=1,\cdots ,n\\ & z_i\in \left[-1,1\right],& & \forall i=1,\cdots ,n\end{array}$$

其中，目标函数即为最大化收益，${\mu }_i+{\sigma }_i{z}_i$即为第$i$支股票的真实收益。

上述模型还有一个学术名称，叫Robust Counterpart。具体我也没看到官方翻译，我姑且翻译它为鲁棒对应模型。

注意，上面的写法是一种完全展开的形式，其中，不确定集即为

$$\begin{array}{rlrl}& \sum _{i=1}^n|z_i|⩽\Gamma & & \\ & z_i\in \left[-1,1\right],& & \forall i=1,\cdots ,n\end{array}$$
围成的可行域，我们一般将不确定集写作
Z = { z ∣ ∑ i = 1 n ∣ z i ∣ ⩽ Γ , z i ∈ [ − 1 , 1 ] , ∀ i = 1 , ⋯ , n } \mathcal{Z} = \{z|\sum_{i=1}^n{|z_i|}\leqslant \Gamma, z_i\in \left[ -1,1 \right] ,\forall i=1,\cdots ,n\} Z={z∣i=1∑n​∣zi​∣⩽Γ,zi​∈[−1,1],∀i=1,⋯,n}

这里吐槽一句：阻碍初学者学习鲁棒优化的一个重要原因就是鲁棒优化中各种奇怪的符号，其实搞定了符号表示什么，已经解决了很大一部分障碍了。

上述模型也可以等价的写成另外一种形式

$$\begin{array}{rlrl}\underset{x\in {\mathbb{R}}^n}{\max }\underset{z:|z_i|⩽1,{\sum }_i|z_i|⩽\Gamma }{\min }& \sum _{i=1}^n\left({\mu }_i+{\sigma }_i{z}_i\right)x_i& & \\ & \sum _{i=1}^n{x}_i=100\%& & \\ & x_i⩾0,& & \forall i=1,2,\cdots ,n\end{array}$$
或者进一步等价为
$$\begin{array}{rlrl}\underset{x\in {\mathbb{R}}^n}{\max }\underset{z\in \mathcal{Z}}{\min }& \sum _{i=1}^n\left({\mu }_i+{\sigma }_i{z}_i\right)x_i& & \\ & \sum _{i=1}^n{x}_i=100\%& & \\ & x_i⩾0,& & \forall i=1,2,\cdots ,n\end{array}$$
其中， Z = { z ∣ ∑ i = 1 n ∣ z i ∣ ⩽ Γ , z i ∈ [ − 1 , 1 ] , ∀ i = 1 , ⋯ , n } \mathcal{Z} = \{z|\sum_{i=1}^n{|z_i|}\leqslant \Gamma, z_i\in \left[ -1,1 \right] ,\forall i=1,\cdots ,n\} Z={z∣∑i=1n​∣zi​∣⩽Γ,zi​∈[−1,1],∀i=1,⋯,n}.
上面是几种鲁棒优化模型的等价写法，希望初学者们特别注意，看明白如上的各种形式，在今后读论文的时候会排除很多障碍。

该模型可以分为两个阶段。

• 第一阶段，将$x_i$看做常量，变量只有$z_i$，此时最小化目标函数，由于目标函数是收益，因此最小化目标函数就相当于得到worst-case，其实，第一阶段就是决策不确定性变量
  $z_i$ 的取值。
• 第二阶段，根据第一阶段的结果，决策$x_i$，使得在worst-case下的收益最大化。

Matlab调用ROME包求解Robust Counterpart

求解鲁棒优化一般有两大类方法：

1. 将鲁棒优化模型重构（reformulate）为一阶段线性规划，然后直接调用求解器求解；
2. 直接使用鲁棒优化求解器(ROME, RSOME等)建模求解。

至于两种方法的适用情况，这里我们先不做详细介绍，我们直接傻瓜式使用matlab调用ROME求解上述鲁棒优化模型。

求解之前，我们先介绍一下将要用到的几个函数：

• newvar: 用于创建决策变量，例如newvar z(n) uncertain;这里，uncertain为变量类型，共有下面几种选项：

  • nonneg: 连续非负决策变量
  • nonneg integer: 非负整数变量
  • integer: 整数变量
  • uncertain: 不确定变量，这个是鲁棒优化特有的，也就是反映不确定参数取值的辅助决策变量$z$。

• h = rome_begin;: 创建模型对象

• rome_maximize(6200 * Drug1 + 6900*Drug2 - (100 * Raw1 + 199 * Raw2 + 700 * Drug1 + 800 * Drug2));：创建目标函数

• rome_constraint(Raw1 + Raw2 <= 1000);：添加约束

• h.solve;：求解模型

• h.objective： 获得目标函数

• h.eval(Drug1);:获取最优解中，变量Drug1的值

关于不确定变量的创建，以下语句非常重要：

-newvar z(n) uncertain;：创建不确定变量

• rome_box(z, -1, 1); rome_constraint(norm1(z) <= Gamma);：创建不确定集

下面我们给出完整代码

% begin rome
clc;
clear;
h = rome_begin;n = 150;
mu = zeros(n, 1);
sigma = zeros(n, 1);
r = zeros(n, 2);
for i = 1:nmu(i, 1) = 0.15 + i * 0.05/150;sigma(i, 1) = (0.05/450) * sqrt(2 * i * n *(n + 1));r(i, 1) = mu(i, 1) - sigma(i, 1);r(i, 2) = mu(i, 1) + sigma(i, 1);
end
Gamma = 4; % 最大可以出现的偏离数
% declare uncertain parameters
newvar z(n) uncertain;
rome_box(z, -1, 1);
rome_constraint(norm1(z) <= Gamma);% Portfolio weights
newvar x(n) nonneg;
newvar y nonneg;% Objective to maximize the return
rome_maximize(y); %注意这里需要用到点乘% constraint to invest all the wealth avaiable
rome_constraint(sum(x) == 1);
rome_constraint(y <= (mu + sigma.*z)' * x);% solve the model
h.solve;
optobj_rob = h.objective;
xx_rob = h.eval(x);
y_rob = h.eval(y);
zz = h.eval(norm1(z))A5 = xx_rob;
save A5;Sum = 0;
for i = 1:length(r(:,1))if(i > length(r(:, 1)) - Gamma)Sum = Sum + r(i, 1) * xx_rob(i, 1);elseSum = Sum + mu(i, 1) * xx_rob(i, 1);end
endfor i = 1:nif(rem(i, 6) == 0)fprintf('#%3d -- %5.4f      \n', i, mu(i, 1));elsefprintf('#%3d -- %5.4f      ', i, mu(i, 1));end
end
fprintf('\n\n-----------------------------------------以下是r的变动范围----------------------------\n')
for i = 1:nif(rem(i, 6) == 0)fprintf('#%3d --[%7.4f, %7.4f]     \n', i, r(i, 1),r(i, 2));elsefprintf('#%3d --[%7.4f, %7.4f]     ', i, r(i, 1),r(i, 2));end
endfprintf('\n\n-----------------------------------------以下 是最优解----------------------------\n')
for i = 1:nif(rem(i, 6) == 0)fprintf('#%3d -- %5.4f      \n', i, xx_rob(i, 1));elsefprintf('#%3d -- %5.4f      ', i, xx_rob(i, 1));end
end
fprintf('Objective:\t  %10.3f\n', optobj_rob);
rome_end;

最后的最优解为：

-----------------------------------------以下 是最优解----------------------------
#  1 -- 0.0000      #  2 -- 0.0000      #  3 -- 0.0000      #  4 -- 0.0000      #  5 -- 0.0000      #  6 -- 0.0000      
#  7 -- 0.0000      #  8 -- 0.0000      #  9 -- 0.0000      # 10 -- 0.0000      # 11 -- 0.0000      # 12 -- 0.0000      
# 13 -- 0.0000      # 14 -- 0.0000      # 15 -- 0.0000      # 16 -- 0.0000      # 17 -- 0.0000      # 18 -- 0.0000      
# 19 -- 0.0000      # 20 -- 0.0000      # 21 -- 0.0000      # 22 -- 0.0000      # 23 -- 0.0000      # 24 -- 0.0000      
# 25 -- 0.0000      # 26 -- 0.0000      # 27 -- 0.0000      # 28 -- 0.0000      # 29 -- 0.0000      # 30 -- 0.0000      
# 31 -- 0.0000      # 32 -- 0.0000      # 33 -- 0.0000      # 34 -- 0.0000      # 35 -- 0.0000      # 36 -- 0.0000      
# 37 -- 0.0000      # 38 -- 0.0000      # 39 -- 0.0000      # 40 -- 0.0000      # 41 -- 0.0000      # 42 -- 0.0000      
# 43 -- 0.0000      # 44 -- 0.0000      # 45 -- 0.0000      # 46 -- 0.0000      # 47 -- 0.0000      # 48 -- 0.0000      
# 49 -- 0.0000      # 50 -- 0.0000      # 51 -- 0.0000      # 52 -- 0.0000      # 53 -- 0.0000      # 54 -- 0.0000      
# 55 -- 0.0000      # 56 -- 0.0000      # 57 -- 0.0000      # 58 -- 0.0000      # 59 -- 0.0000      # 60 -- 0.0000      
# 61 -- 0.0000      # 62 -- 0.0000      # 63 -- 0.0000      # 64 -- 0.0000      # 65 -- 0.0000      # 66 -- 0.0000      
# 67 -- 0.0000      # 68 -- 0.0000      # 69 -- 0.0000      # 70 -- 0.0000      # 71 -- 0.0000      # 72 -- 0.0155      
# 73 -- 0.0154      # 74 -- 0.0152      # 75 -- 0.0151      # 76 -- 0.0150      # 77 -- 0.0149      # 78 -- 0.0149      
# 79 -- 0.0148      # 80 -- 0.0147      # 81 -- 0.0146      # 82 -- 0.0145      # 83 -- 0.0144      # 84 -- 0.0143      
# 85 -- 0.0142      # 86 -- 0.0141      # 87 -- 0.0141      # 88 -- 0.0140      # 89 -- 0.0139      # 90 -- 0.0138      
# 91 -- 0.0137      # 92 -- 0.0137      # 93 -- 0.0136      # 94 -- 0.0135      # 95 -- 0.0135      # 96 -- 0.0134      
# 97 -- 0.0133      # 98 -- 0.0132      # 99 -- 0.0132      #100 -- 0.0131      #101 -- 0.0131      #102 -- 0.0130      
#103 -- 0.0129      #104 -- 0.0129      #105 -- 0.0128      #106 -- 0.0127      #107 -- 0.0127      #108 -- 0.0126      
#109 -- 0.0126      #110 -- 0.0125      #111 -- 0.0124      #112 -- 0.0124      #113 -- 0.0123      #114 -- 0.0123      
#115 -- 0.0122      #116 -- 0.0122      #117 -- 0.0121      #118 -- 0.0121      #119 -- 0.0120      #120 -- 0.0120      
#121 -- 0.0119      #122 -- 0.0119      #123 -- 0.0118      #124 -- 0.0118      #125 -- 0.0117      #126 -- 0.0117      
#127 -- 0.0116      #128 -- 0.0116      #129 -- 0.0115      #130 -- 0.0115      #131 -- 0.0115      #132 -- 0.0114      
#133 -- 0.0114      #134 -- 0.0113      #135 -- 0.0113      #136 -- 0.0112      #137 -- 0.0112      #138 -- 0.0112      
#139 -- 0.0111      #140 -- 0.0111      #141 -- 0.0110      #142 -- 0.0110      #143 -- 0.0110      #144 -- 0.0109      
#145 -- 0.0109      #146 -- 0.0109      #147 -- 0.0108      #148 -- 0.0108      #149 -- 0.0107      #150 -- 0.0107      
Objective:	       0.174

可以得知，最优解为17.4%，即考虑$\Gamma =4$的情况下，worst-case下的最优回报率为17.4%。

如果我们调整$Gamma=10$，则最优解变化为

-----------------------------------------以下 是最优解----------------------------
#  1 -- 0.0000      #  2 -- 0.0000      #  3 -- 0.0000      #  4 -- 0.0000      #  5 -- 0.0000      #  6 -- 0.0000      
#  7 -- 0.0000      #  8 -- 0.0000      #  9 -- 0.0000      # 10 -- 0.0000      # 11 -- 0.0000      # 12 -- 0.0000      
# 13 -- 0.0000      # 14 -- 0.0000      # 15 -- 0.0000      # 16 -- 0.0000      # 17 -- 0.0000      # 18 -- 0.0000      
# 19 -- 0.0000      # 20 -- 0.0000      # 21 -- 0.0000      # 22 -- 0.0000      # 23 -- 0.0000      # 24 -- 0.0000      
# 25 -- 0.0000      # 26 -- 0.0000      # 27 -- 0.0000      # 28 -- 0.0000      # 29 -- 0.0000      # 30 -- 0.0000      
# 31 -- 0.0133      # 32 -- 0.0131      # 33 -- 0.0129      # 34 -- 0.0127      # 35 -- 0.0125      # 36 -- 0.0124      
# 37 -- 0.0122      # 38 -- 0.0120      # 39 -- 0.0119      # 40 -- 0.0117      # 41 -- 0.0116      # 42 -- 0.0114      
# 43 -- 0.0113      # 44 -- 0.0112      # 45 -- 0.0111      # 46 -- 0.0109      # 47 -- 0.0108      # 48 -- 0.0107      
# 49 -- 0.0106      # 50 -- 0.0105      # 51 -- 0.0104      # 52 -- 0.0103      # 53 -- 0.0102      # 54 -- 0.0101      
# 55 -- 0.0100      # 56 -- 0.0099      # 57 -- 0.0098      # 58 -- 0.0097      # 59 -- 0.0097      # 60 -- 0.0096      
# 61 -- 0.0095      # 62 -- 0.0094      # 63 -- 0.0093      # 64 -- 0.0093      # 65 -- 0.0092      # 66 -- 0.0091      
# 67 -- 0.0091      # 68 -- 0.0090      # 69 -- 0.0089      # 70 -- 0.0089      # 71 -- 0.0088      # 72 -- 0.0087      
# 73 -- 0.0087      # 74 -- 0.0086      # 75 -- 0.0086      # 76 -- 0.0085      # 77 -- 0.0084      # 78 -- 0.0084      
# 79 -- 0.0083      # 80 -- 0.0083      # 81 -- 0.0082      # 82 -- 0.0082      # 83 -- 0.0081      # 84 -- 0.0081      
# 85 -- 0.0080      # 86 -- 0.0080      # 87 -- 0.0079      # 88 -- 0.0079      # 89 -- 0.0079      # 90 -- 0.0078      
# 91 -- 0.0078      # 92 -- 0.0077      # 93 -- 0.0077      # 94 -- 0.0076      # 95 -- 0.0076      # 96 -- 0.0076      
# 97 -- 0.0075      # 98 -- 0.0075      # 99 -- 0.0075      #100 -- 0.0074      #101 -- 0.0074      #102 -- 0.0073      
#103 -- 0.0073      #104 -- 0.0073      #105 -- 0.0072      #106 -- 0.0072      #107 -- 0.0072      #108 -- 0.0071      
#109 -- 0.0071      #110 -- 0.0071      #111 -- 0.0070      #112 -- 0.0070      #113 -- 0.0070      #114 -- 0.0069      
#115 -- 0.0069      #116 -- 0.0069      #117 -- 0.0069      #118 -- 0.0068      #119 -- 0.0068      #120 -- 0.0068      
#121 -- 0.0067      #122 -- 0.0067      #123 -- 0.0067      #124 -- 0.0067      #125 -- 0.0066      #126 -- 0.0066      
#127 -- 0.0066      #128 -- 0.0066      #129 -- 0.0065      #130 -- 0.0065      #131 -- 0.0065      #132 -- 0.0065      
#133 -- 0.0064      #134 -- 0.0064      #135 -- 0.0064      #136 -- 0.0064      #137 -- 0.0063      #138 -- 0.0063      
#139 -- 0.0063      #140 -- 0.0063      #141 -- 0.0062      #142 -- 0.0062      #143 -- 0.0062      #144 -- 0.0062      
#145 -- 0.0062      #146 -- 0.0061      #147 -- 0.0061      #148 -- 0.0061      #149 -- 0.0061      #150 -- 0.0061      
Objective:	       0.160

可以看到，最优解降低到16.0%。因为$\Gamma$越大，表示越保守。

好啦，今天的例子就讲解到这里，下次我们通过reformulation的方法，来求解该鲁棒优化模型，敬请期待。

参考文献

1. E. Guslitzer A. Ben-Tal, A. Goryashko and A. Nemirovski. Robust optimization.2009.
2. https://robustopt.com/
3. ROME: Joel Goh and Melvyn Sim. Robust Optimization Made Easy with ROME. Operations Research, 2011, 59(4), pp.973-985.
4. Joel Goh and Melvyn Sim. Distributionally Robust Optimization and its Tractable Approximations. Operations Research, 2010, 58(4), pp. 902-917.
5. Chen, Zhi, and Peng Xiong. 2021. RSOME in Python: an open-source package for robust stochastic optimization made easy. Optimization Online.
6. Chen, Zhi, Melvyn Sim, Peng Xiong. 2020. Robust stochastic optimization made easy with RSOME. Management Science 66(8) 3329–3339.