鲁棒优化|基于ROME编程入门鲁棒优化：以一个例子引入（上）
鲁棒优化是处理数学规划中处理不确定性的重要方法论。这里我们以几个简单的例子来介绍鲁棒优化。同时，我们结合鲁棒优化工具包Robust Optimization Made Easy (ROME)来编程实现，方便读者理解和复现。
线性规划中的不确定性
在生产实际中，一些生产计划、运输问题可以建模为线性规划进行求解，从而达到利润的最大化或者成本的最小化。但是这些建模是基于一个重要假设：
这些线性规划模型，一般可以表示为下面的形式：
目标函数 max c T x \max c^T x maxcTx 表示要最大化收益，约束 A x ⩽ b Ax \leqslant b Ax⩽b 表示生产所用的原材料要小于等于可用的生产资源。
但是在实际生产生活中，价值系数
c
c
c，资源向量
b
b
b，生产要素矩阵
A
A
A都有可能是不确定的，也就是具有不确定性。例如，股票投资中，股票的预期收益是不确定的，无法提前预知的。再比如实际生产中，生产实际消耗的资源(对应
A
A
A)也是不确定的，会由于施工精度，施工条件和工人的操作产生偏差。
如果我们在建模中考虑这些不确定的因素，则价值系数 c c c，资源向量 b b b和生产要素矩阵 A A A则可以分别表示为 c ~ , b ~ , A ~ \tilde{c}, \tilde{b}, \tilde{A} c~,b~,A~. 因此，考虑不确定因素的线性规划模型可以写为
当然啦，实际中，上述不确定也许不会同时出现。或者，决策者不想同时考虑。因此有时决策者只会考虑其中一种(例如只考虑目标函数系数 c c c不确定)，或者多种。
处理这种带有不确定性的线性规划模型，一般有两大类方法论：
我们这期推文主要来讲解鲁棒优化。
一个生产的例子
本例子来自文献[1]。
我们将其建模为一个线性规划模型，如下
Matlab+ROME求解
网址：https://robustopt.com/
以及RSOME对应的论文
在这里向前辈大佬们致以崇高的敬意！！！
我们用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 solution
fprintf('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。本例说明了，除了约束以外，目标函数中有不确定性因素，也会影响到问题的解。
有 n n n种股票，在确定性问题下，每种股票的投资回报率为 r i r_i ri ，我们需要决策每种股票的投资占比，从而使得收益最大。确定型模型如下。
但是投资回报率是有波动的，确定型问题中的回报率 r i r_i ri只是一个期望值。
不确定回报率
为了建立robust optimization模型，我们考虑投资回报率 r i r_i ri是不确定的变量，即
我们考虑以下算例：有150种股票，
n
=
150
n=150
n=150，每种股票的均值和方差都是随着股票的ID依次递增。并且
μ
i
=
0.15
+
i
0.05
150
σ
i
=
0.05
450
2
n
i
(
n
+
1
)
\begin{aligned} \mu _i&=0.15+i\frac{0.05}{150} \\ \sigma _i&=\frac{0.05}{450}\sqrt{2ni\left( n+1 \right)} \end{aligned}
μiσi=0.15+i1500.05=4500.052ni(n+1)
因此，可以得到，150种股票的投资回报率的
μ
i
,
σ
i
\mu_i, \sigma_i
μi,σi以及投资回报率
r
~
i
\tilde{r}_i
r~i的波动范围。
- 150种股票的投资回报率的 μ i \mu_i μ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种股票的投资回报率的 σ i \sigma_i σ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种股票的投资回报率的波动范围如下：
r ~ i ∈ [ 0.15 + i 0.05 150 − 0.05 450 2 n i ( n + 1 ) , 0.15 + i 0.05 150 + 0.05 450 2 n i ( n + 1 ) ] \begin{aligned} \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{aligned} r~i∈[0.15+i1500.05−4500.052ni(n+1),0.15+i1500.05+4500.052ni(n+1)]
-----------------------------------------以下是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）。
描绘不确定集的方法有多种，在基本的鲁棒优化中，常见的有：
预算不确定集(Budgeted uncertainty set)凸包不确定集(Convex hull uncertainty set)盒子不确定集(Box uncertainty set)CVaR不确定集(CVaR uncertainty Set)椭球不确定集(Ellipsoidal uncertainty Set)
另外，我们也这里也放出鲁棒优化的经典文献：
本文中，我们以预算不确定集(Budgeted uncertainty set)为例，来对上述问题进行建模。我们设置：至多有
Γ
\Gamma
Γ个股票的投资回报率会偏离其均值
μ
\mu
μ。为了方便叙述，取
Γ
=
4
\Gamma = 4
Γ=4。这里的
Γ
\Gamma
Γ就是预算，通俗的理解，就是我们考虑明天股票走势最坏的情况，150支股票的收益相较平均值的波动的比例的总和，不会高于400%。当然了，如果设置
Γ
=
150
\Gamma=150
Γ=150，表示你考虑了最点背的情况，你假想你是全世界最点背的人，明天有可能150支都跌倒最低点，哈哈。反之，如果设置
Γ
\Gamma
Γ较小，说明你是一个乐观主义者，你比较看好明天股票的走势。
基于此，我们给出鲁棒优化模型如下：
其中，目标函数即为最大化收益， μ i + σ i z i \mu _i+\sigma _iz_i μi+σizi即为第 i i i支股票的真实收益。
上述模型还有一个学术名称，叫Robust Counterpart。具体我也没看到官方翻译，我姑且翻译它为鲁棒对应模型。
注意，上面的写法是一种完全展开的形式，其中，不确定集即为
上述模型也可以等价的写成另外一种形式
该模型可以分为两个阶段。
- 第一阶段，将
x
i
x_i
xi看做常量，变量只有
z
i
z_i
zi，此时最小化目标函数，由于目标函数是收益，因此最小化目标函数就相当于得到
worst-case，其实，第一阶段就是决策不确定性变量
z i z_i zi 的取值。 - 第二阶段，根据第一阶段的结果，决策
x
i
x_i
xi，使得在
worst-case下的收益最大化。
Matlab调用ROME包求解Robust Counterpart
至于两种方法的适用情况，这里我们先不做详细介绍，我们直接傻瓜式使用matlab调用ROME求解上述鲁棒优化模型。
求解之前，我们先介绍一下将要用到的几个函数：
-
newvar: 用于创建决策变量，例如newvar z(n) uncertain;这里，uncertain为变量类型，共有下面几种选项：nonneg: 连续非负决策变量nonneg integer: 非负整数变量integer: 整数变量uncertain: 不确定变量，这个是鲁棒优化特有的，也就是反映不确定参数取值的辅助决策变量z z 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:n
mu(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);
else
Sum = Sum + mu(i, 1) * xx_rob(i, 1);
end
end
for i = 1:n
if(rem(i, 6) == 0)
fprintf('#%3d -- %5.4f \n', i, mu(i, 1));
else
fprintf('#%3d -- %5.4f ', i, mu(i, 1));
end
end
fprintf('\n\n-----------------------------------------以下是r的变动范围----------------------------\n')
for i = 1:n
if(rem(i, 6) == 0)
fprintf('#%3d --[%7.4f, %7.4f] \n', i, r(i, 1),r(i, 2));
else
fprintf('#%3d --[%7.4f, %7.4f] ', i, r(i, 1),r(i, 2));
end
end
fprintf('\n\n-----------------------------------------以下 是最优解----------------------------\n')
for i = 1:n
if(rem(i, 6) == 0)
fprintf('#%3d -- %5.4f \n', i, xx_rob(i, 1));
else
fprintf('#%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%，即考虑
Γ
=
4
\Gamma=4
Γ=4的情况下，worst-case下的最优回报率为17.4%。
如果我们调整 G a m m a = 10 Gamma =10 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的方法，来求解该鲁棒优化模型，敬请期待。
参考文献
- E. Guslitzer A. Ben-Tal, A. Goryashko and A. Nemirovski. Robust optimization.2009.
- https://robustopt.com/
- 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.
- 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.
