鲁棒优化 | 以Coding入门鲁棒优化：以一个例子引入（二）

IT圈老男孩1 · 发表于 2022-3-31 19:18

作者：刘兴禄，清华大学，清华-伯克利深圳学院，博士在读 邮箱：hsinglul@163.com
1.投资组合的例子

上一篇我们介绍了一个鲁棒优化求解投资组合的例子，文章链接为：
其中，该例子的Robust Counterpart模型如下：

$\begin{aligned} \underset{x\in \mathbb{R}^n}{\max}\quad\underset{z\in \mathcal{Z}}{\min} \qquad &\sum_{i=1}^n{\left( \mu _i+\sigma _iz_i \right) x_i}&& \\ &\sum_{i=1}^n{x_i}=100\%&& \\ &x_i\geqslant 0, &&\qquad \forall i=1,2,\cdots ,n \end{aligned}$
其中， $\mathcal{Z} = {z|\sum_{i=1}^n{|z_i|}\leqslant \Gamma, z_i\in \left[ -1,1 \right] ,\forall i=1,\cdots ,n}。$

该例中，我们采用了预算不确定集(Budgeted uncertainty set)。
我们也在上一篇文章中提到，求解鲁棒优化模型的两大类方法：

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

其中，第2类方法，我们已经在器上一篇文章中做了详细解释，本文我们来介绍第一类求解方法：reformulation。
reformulation又分为很多种，本文仅仅介绍两种：

利用对偶进行reformulation
利用上镜图(epigraph)进行reformulation

两种方法在操作上都需要对偶，有很多类似之处，但是思想却是不同的。
2.鲁棒优化模型的reformulation: 利用对偶进行reformulation

2.1利用对偶进行reformulation

上述模型其实是一个bi-level的模型，也就是 $\underset{x}{\max} \,\,\underset{z}{\min}$ 。上述鲁棒优化模型可以在一定程度上理解为是一个赌场博弈问题(只是为了方便理解，并不一定严谨)：

inner level或lower level的决策：不确定变量，可以看做是赌场的庄家，庄家不想让玩家赚钱，所以他控制收益 $\tilde{r}$ 波动(也就是控制风险)，尽可能让你赚的少(可以理解为给你制造worst case)
outer level或upper level的决策：投资决策。作为玩家，我们需要跟庄家博弈，我们要最大化收益，使得在worst case下的收益最大化。

在这种设定下，庄家和玩家的目标函数表达式是一致的，都是 $\sum_{i=1}^n{\left( \mu i+\sigma _iz_i \right) x_i}$ ，但是二者的目的(或者说，优化方向)是相反的，庄家想要 $\underset{z}{\min} \sum{i=1}^n{\left( \mu i+\sigma _iz_i \right) x_i}$ ，但是玩家想要 $\underset{x}{\max} \sum{i=1}^n{\left( \mu _i+\sigma _iz_i \right) x_i}$ 。
此外，其实这里有一个隐含的假设：庄家的决策和玩家的决策没有先后顺序。都是基于预判对方的决策，直接做出了自己的决策，后续无论庄家怎么变，玩家都不再改变自己的决策。而不是：庄家真正的先做决策，玩家再去根据对方的决策调整自己的策略。

小拓展：还有一种bi-level的形式，是两个主体的目标函数或优化方向不一致(也可以目标函数和优化方向都不一致)。这种形式是更为一般的bi-level决策，是典型的的博弈问题。这种情况一般有一个leader和一个follower组成。比如
leader的决策： $\begin{aligned}\max &\sum c_e x_e\ s.t. &\text{ leader的约束}\end{aligned}$
follower的决策： $\begin{aligned}\max &\sum d_e y_e\ s.t. &\text{ follower的约束}\end{aligned}$
比如，共享出行中，平台想最大化整个系统的利润，但是对于每个个体而言，他们只想最大化自己的收入。这两个目标函数是不同的，但是所做的决策之间却有一部分耦合。关于这一部分，本文只简单提及，不再做详细介绍。

接着介绍利用对偶进行reformulation。这种方法的主要思想是

基于bi-level的模型 $\max \, \min$ ，对inner level模型取对偶，将bi-level模型转化成为 $\max\,\max$ ，进而等价为一阶段问题进行求解

注意：本文介绍的这种基于对偶的reformulation方法有两个非常强的前提条件：
1. inner level的模型需要是线性模型，不能有binary和integer的变量以及非线性项。原因是:
  -a) MIP是不能以这种方法直接对偶的（但是可以采用其他方法处理，不过也并不是完全等价转化）。
  -b) 有非线性项的情况下，如果决策变量都是连续的，可以用KKT条件等价转化。
2. 强对偶是成立的。只有强对偶是成立的，这种转化才是等价的。验证强对偶成立与否的方法很简单，对偶完成以后，看看对偶问题是否存在可行解，如果存在至少一个有解的可行解，则强对偶必然成立。这也是基于对偶定理得来的，即：如果一个问题有可行解，并且目标函数是有界的(因此也有最优解)，那么另一个问题也有可行解和有界的目标函数以及最优解，此时强对偶理论和弱对偶理论都是成立的。

接下来，我们来进行reformulation。
首先我们将模型中的绝对值进行线性化(即 $\sum_{i=1}^n{|z_i|}\leqslant \Gamma$ )。

我们引入两组辅助变量 $z_i^{+}, z_i^{-}, \forall i =1,2,\cdots, n$ 。其中
$\begin{aligned} &z_i^{+} = \max{0, z_i} \ \\ &z_i^{-} = \max{0, -z_i} \end{aligned}$
因此，就有
$\begin{aligned} &z_i = z_i^{+} - z_i^{-} \\ &|z_i| = z_i^{+} + z_i^{-} \end{aligned}$

这个比较容易理解，例如 $z_i = -0.5$ ，则 $z_i^{+} = 0, z_i^{-} = 0.5$ ，因此有 $z_i = z_i^{+} - z_i^{-} = 0 - 0.5 = -0.5$ ； $|z_i| = z_i^{+} + z_i^{-} = 0 + 0.5 = 0.5$ 。

拓展之前的推文中，有小伙伴纠结这个转化是否等价。一个肯定的回答是：完全等价。这里我给出证明：
命题：
$\begin{aligned} &z^+=\max \left\{ 0,z \right\} \,\, \left( 1 \right) \\ &z^-=\max \left\{ 0,-z \right\} \,\, \left( 2 \right) \end{aligned}$
和
$\begin{aligned} &z=z^+-z^-\,\, \left( 3 \right) \\ &|z|=z^++z^-\,\,\left( 4 \right) \end{aligned}$
互为充要条件。

证明：充分性：根据(1)(2)，我们有
$\begin{aligned} &z=z^+-z^-\,\, \left( 3 \right) \\ &|z|=z^++z^-\,\,\left( 4 \right) \end{aligned}$
充分性得证必要性：根据(3)(4)我们有 $\begin{aligned} &\left( 3 \right) +\left( 4 \right) =z+|z|=2z^+\,\, \rightarrow \,\,z^+=\frac{z+|z|}{2} \\ &\left( 3 \right) -\left( 4 \right) =z-|z|=-2z^-\,\, \rightarrow \,\,z^-=-\frac{z-|z|}{2} \\ &\text{情况一：}if\,\,z\geqslant 0, z^+=\frac{z+|z|}{2}=\frac{z+z}{2}=z=\max \left\{ 0,z \right\} \\ &z^-=-\frac{z-|z|}{2}=-\frac{z-z}{2}=0=\max \left\{ 0,-z \right\} \\ &\text{情况二：}if\,\,z\leqslant 0, z^+=\frac{z+|z|}{2}=\frac{z-z}{2}=0=\max \left\{ 0,z \right\} \\ &z^-=-\frac{z-|z|}{2}=-\frac{z-\left( -z \right)}{2}=-\frac{2z}{2}=-z=\max \left\{ 0,-z \right\} \end{aligned}$
因此，必要性得证。综上，(1)(2)和(3)(4)互为充要条件，以上转化完全等价。

我们就借助这一步转化，将模型进行改写为

$\begin{aligned} \underset{x\in \mathbb{R}^n}{\max}\quad \underset{z_i^{+}, z_i^{-} }{\min}\quad &\sum_{i=1}^n{\left[ \mu _i+\sigma _i (z_i^{+} - z_i^{-}) \right] x_i} && \\ &\sum_{i=1}^n{x_i}=100\% && \\ & \sum_{i=1}^n{ (z_i^{+} + z_i^{-}) }\leqslant \Gamma && \\ &x_i\geqslant 0, \qquad &&\forall i=1,\cdots ,n \\ &z_i^{+}, z_i^{-}\in \left[ 0,1 \right] , \qquad &&\forall i=1,\cdots ,n \end{aligned}$
注意：由于在鲁棒优化中，决策者往往并不关心worst case究竟长什么样子，只在意自己的决策应该怎么做，因此决策者也就不关心的取值。

下面我们来进行对偶。我们将内层模型由$\min$对偶成$\max$问题，这样就可以把双层(bi-level)模型转化为单层(single-level)模型。这里注意，对于内层模型(决策worst-case)而言，决策变量就只有 $z_i^{+}, z_i^{-}$ ，而第二层的决策变量就要被视为固定值(固定值)，也就是视为参数。
因此，内层模型就变化为(因为视为固定值，故相关约束全部可以暂时删除)

$\begin{aligned} \underset{z_{i}^{+},z_{i}^{-}}{\min}\quad &C+\sum_{i=1}^n{\sigma _ix_iz_{i}^{+}-}\sum_{i=1}^n{\sigma _ix_iz_{i}^{-}} && \\ &\sum_{i=1}^n{\left( z_{i}^{+}+z_{i}^{-} \right)}\leqslant \Gamma && \qquad \rightarrow \quad \pi \\ &z_{i}^{+}\leqslant 1,\qquad \forall i=1,\cdots ,n && \qquad \rightarrow \quad \lambda_i \\ &z_{i}^{-}\leqslant 1,\qquad \forall i=1,\cdots ,n && \qquad \rightarrow \quad \theta_i \\ &z_{i}^{+},z_{i}^{-}\geqslant 0 \qquad\forall i=1,\cdots ,n && \end{aligned}$

其中$C$为常数，且$C = \sum_{i=1}^n{\mu _ix_i}$ $C = \sum_{i=1}^n{\mu _ix_i}$ ，且约束(1)、约束(2)和约束(3)的对偶变量分别是 $\pi$ 和 $\lambda$ 和 $\theta$ 。因此，内层模型的对偶问题为:

$\begin{aligned} \underset{\pi ,\lambda ,\theta}{\max}\qquad &\pi \Gamma +\sum_{i=1}^n{\lambda _i}+\sum_{i=1}^n{\theta _i}&& \\ &\pi +\lambda _i\leqslant \sigma _ix_i, &&\qquad \forall i=1,\cdots ,n \\ &\pi +\theta _i\leqslant -\sigma _ix_i, &&\qquad \forall i=1,\cdots ,n \\ &\pi ,\lambda _i,\theta _i\leqslant 0, &&\qquad \forall i=1,\cdots ,n \end{aligned}$
这个对偶并不难写出来，此处有困难的读者朋友可以参考《运小筹》公众号的往期推文：

这里，我们来验证模型的强对偶是否成立。我们发现当所有 $\pi, \lambda, \theta$ 均为0，是对偶模型的一个可行解，且目标函数值为0。也就是说，所有变量均为0，是一个有界的可行解。因此，强对偶成立。
结合内层的对偶问题，原问题可以转化为单层模型，如下:

$\begin{aligned} \underset{x,\pi ,\lambda ,\theta}{\max}\qquad &\sum_{i=1}^n{\mu _ix_i}+\pi \Gamma +\sum_{i=1}^n{\lambda _i}+\sum_{i=1}^n{\theta _i} && \\ &\sum_{i=1}^n{x_i}=1 && \\ &\pi +\lambda _i\leqslant \sigma _ix_i, &&\qquad \forall i=1,\cdots ,n \\ &\pi +\theta _i\leqslant -\sigma _ix_i, &&\qquad \forall i=1,\cdots ,n \\ &\pi ,\lambda _i,\theta _i\leqslant 0, &&\qquad \forall i=1,\cdots ,n \\ &x_i\geqslant 0,&& \qquad \forall i=1,\cdots ,n \end{aligned}$

2.2Python调用gurobi求解对偶reformulation后的模型

我们来验证该模型是否正确
# author: Xinglu Liu

from gurobipy import *
import numpy as np

# stock number
stock_num = 150
Gamma = 4
# input parameters
mu = [0] * stock_num
sigma = [0] * stock_num
reward = [ [0, 0] ] * stock_num

# initialize the parameters
for i in range(stock_num):
mu = 0.15 + (i + 1) * (0.05/150)
sigma = (0.05/450) * math.sqrt(2 * stock_num * (i + 1) * (stock_num + 1))
reward[0] = mu - sigma
reward[1] = mu + sigma
mu = np.array(mu)
sigma = np.array(sigma)
reward = np.array(reward)

model = Model(&#39;Robust portfolio Dual&#39;)
x = {}
lam = {}
theta = {}
pi = {}
pi[0] = model.addVar(lb = -GRB.INFINITY, ub = 0, obj = 0, vtype = GRB.CONTINUOUS, name = &#39;pi&#39;)
for i in range(stock_num):
x = model.addVar(lb = 0, ub = 1, obj = 0, vtype = GRB.CONTINUOUS, name = &#39;x_&#39; + str(i+1))
lam = model.addVar(lb = -GRB.INFINITY, ub = 0, obj = 0, vtype = GRB.CONTINUOUS, name = &#39;lam_&#39; + str(i+1))
theta = model.addVar(lb = -GRB.INFINITY, ub = 0, obj = 0, vtype = GRB.CONTINUOUS, name = &#39;theta_&#39; + str(i+1))

# set the objective function
obj = LinExpr(0)
obj.addTerms(Gamma, pi[0])
for i in range(stock_num):
obj.addTerms(mu, x)
obj.addTerms(1, lam)
obj.addTerms(1, theta)
model.setObjective(obj, GRB.MAXIMIZE)

# constraints 1 : invest budget
lhs = LinExpr(0)
for i in range(stock_num):
lhs.addTerms(1, x)
model.addConstr(lhs == 1, name = &#39;invest budget&#39;)

# constraints 2: dual constraint
for i in range(stock_num):
lhs = LinExpr(0)
rhs = LinExpr(0)
lhs.addTerms(1, pi[0])
lhs.addTerms(1, lam)
rhs.addTerms(sigma, x)
model.addConstr(lhs <= rhs, name = &#39;dual_cons1_&#39; + str(i+1))

# constraints 3: dual constraint
for i in range(stock_num):
lhs = LinExpr(0)
rhs = LinExpr(0)
lhs.addTerms(1, pi[0])
lhs.addTerms(1, lam)
rhs.addTerms(-sigma, x)
model.addConstr(lhs <= rhs, name = &#39;dual_cons2_&#39; + str(i+1))

model.write(&#39;model.lp&#39;)

model.optimize()

# print out the solution
print(&#39;-----------------------------------------------&#39;)
print(&#39;----------- optimal solution -------------&#39;)
print(&#39;-----------------------------------------------&#39;)
print(&#39;objective : &#39;, round(model.ObjVal, 10))
for key in x.keys():
if(x[key].x > 0):
      print(x[key].varName, &#39;\t = &#39;, round(x[key].x, 10))

x_sol = []
for key in x.keys():
x_sol.append(round(x[key].x, 4))

import matplotlib.pyplot as plt

fig = plt.figure(figsize=(12,8))

plt.plot(range(1, stock_num+1), x_sol,
      linewidth=2, color=&#39;b&#39;)
plt.xlabel(&#39;Stocks&#39;)
plt.ylabel(&#39;Fraction of investment&#39;)
plt.show()求解结果如下：

图2-1

Gurobi Optimizer version 9.0.1 build v9.0.1rc0 (win64)
Optimize a model with 301 rows, 451 columns and 1050 nonzeros
Model fingerprint: 0xb8185bc2
Coefficient statistics:
  Matrix range    [2e-02, 1e+00]
  Objective range  [2e-01, 4e+00]
  Bounds range    [1e+00, 1e+00]
  RHS range       [1e+00, 1e+00]
Presolve removed 150 rows and 150 columns
Presolve time: 0.01s
Presolved: 151 rows, 301 columns, 600 nonzeros

Iteration Objective    Primal Inf. Dual Inf.    Time
   0 2.0000000e-01 2.896358e-01 0.000000e+00    0s
   79 1.7378554e-01 0.000000e+00 0.000000e+00    0s

Solved in 79 iterations and 0.01 seconds
Optimal objective  1.737855434e-01
-----------------------------------------------
----------- optimal solution -------------
-----------------------------------------------
objective :  0.1737855434
x_72    =  0.0154576484
x_73    =  0.015351409
x_74    =  0.0152473305
x_75    =  0.0151453405
x_76    =  0.0150453702
x_77    =  0.0149473537
x_78    =  0.0148512282
x_79    =  0.0147569338
x_80    =  0.0146644129
x_81    =  0.0145736107
x_82    =  0.0144844746
x_83    =  0.0143969543
x_84    =  0.0143110016
x_85    =  0.0142265702
x_86    =  0.0141436157
x_87    =  0.0140620956
x_88    =  0.0139819691
x_89    =  0.0139031968
x_90    =  0.0138257411
x_91    =  0.0137495656
x_92    =  0.0136746355
x_93    =  0.0136009173
x_94    =  0.0135283785
x_95    =  0.0134569882
x_96    =  0.0133867162
x_97    =  0.0133175338
x_98    =  0.0132494129
x_99    =  0.0131823269
x_100 =  0.0131162496
x_101 =  0.0130511562
x_102 =  0.0129870223
x_103 =  0.0129238248
x_104 =  0.0128615409
x_105 =  0.012800149
x_106 =  0.0127396278
x_107 =  0.0126799571
x_108 =  0.0126211171
x_109 =  0.0125630887
x_110 =  0.0125058533
x_111 =  0.0124493932
x_112 =  0.0123936909
x_113 =  0.0123387297
x_114 =  0.0122844933
x_115 =  0.0122309658
x_116 =  0.012178132
x_117 =  0.0121259771
x_118 =  0.0120744865
x_119 =  0.0120236463
x_120 =  0.011973443
x_121 =  0.0119238633
x_122 =  0.0118748944
x_123 =  0.011826524
x_124 =  0.0117787399
x_125 =  0.0117315303
x_126 =  0.0116848839
x_127 =  0.0116387895
x_128 =  0.0115932363
x_129 =  0.0115482139
x_130 =  0.0115037119
x_131 =  0.0114597205
x_132 =  0.0114162299
x_133 =  0.0113732308
x_134 =  0.0113307139
x_135 =  0.0112886703
x_136 =  0.0112470913
x_137 =  0.0112059683
x_138 =  0.0111652931
x_139 =  0.0111250577
x_140 =  0.0110852542
x_141 =  0.0110458748
x_142 =  0.0110069122
x_143 =  0.0109683589
x_144 =  0.010930208
x_145 =  0.0108924524
x_146 =  0.0108550854
x_147 =  0.0108181004
x_148 =  0.0107814908
x_149 =  0.0107452504
x_150 =  0.010709373该结果与上篇推文中直接调用ROME求解的结果相同，说明我们的reformulation是正确的。开心~
3 鲁棒优化模型的reformulation: 利用上镜图reformulation

3.1 利用上镜图reformulation

由于上述模型是bi-level的模型，即 $\max \,\, \min$ 问题，我们可以通过上镜图的方法将其转化为single level的模型(这种做法也是参考一些文献的做法，参见参考文献7)，具体方法为：

引入辅助决策变量 $t$

则上述模型可以转化为

$\begin{aligned} \underset{x\in \mathbb{R}^n}{\max} \quad &t && \\ &t\leqslant \sum_{i=1}^n{\left( \mu _i+\sigma _iz_i \right) x_i}, && \forall \mathbf{z}=(z_1, \cdots, z_n)^T \in \mathcal{Z} \\ &\sum_{i=1}^n{x_i}=100\% && \\ &x_i\geqslant 0,& &\forall i=1,\cdots ,n \end{aligned}$

其中，约束 $t\leqslant \sum_{i=1}^n{\left( \mu i+\sigma _iz_i \right) x_i}, \,\, \forall \mathbf{z}=(z_1, \cdots, z_n)^T \in \mathcal{Z}$ 就等价于 $\underset{z\in \mathcal{Z}}{\min} \sum{i=1}^n{\left( \mu _i+\sigma _iz_i \right) x_i}$ 。

一定要注意约束 $t\leqslant \sum_{i=1}^n{\left( \mu _i+\sigma _iz_i \right) x_i}$ 后面的 $\forall \mathbf{z}=(z_1, \cdots, z_n)^T \in \mathcal{Z}$ ,这个一定不能少。

上述约束其实也可以写为

$\begin{aligned} t\leqslant \underset{z}{\min} \sum_{i=1}^n{\left( \mu _i+\sigma _iz_i \right) x_i}, \,\, \forall \mathbf{z}=(z_1, \cdots, z_n)^T \in \mathcal{Z} \qquad (1) \end{aligned}$
这种写法也上面的写法是等价的。因为，  右端项所有可能的取值，就等价于  右端项的最小值。

实际上，约束(1)自己就是一个小优化问题，并不是能够直接处理的约束。
接下来，我们的任务就是要将(1)转化成为可处理的形式。
我们需要将其作进一步处理。我们将(1)改写成

$\begin{aligned} t\leqslant \underset{z_i\in \left[ -1,1 \right] ,\forall i;\sum{|z_i|}\leqslant \Gamma}{\min}\sum_{i=1}^n{\left( \mu _i+\sigma _iz_i \right) x_i}, \qquad (2) \end{aligned}$

转化的方法，是将 $\underset{z_\in \mathcal{Z}}{\min}\sum_{i=1}^n{\left( \mu _i+\sigma _iz_i \right) x_i}$ 取对偶。方法与上面介绍的直接对偶的方法一致，我们可以直接把结果拿过来。

$\begin{aligned} \underset{\pi ,\lambda ,\theta}{\max}\qquad & \sum_{i=1}^n{\mu _ix_i} + \pi \Gamma +\sum_{i=1}^n{\lambda _i}+\sum_{i=1}^n{\theta _i}&& \\ &\pi +\lambda _i\leqslant \sigma _ix_i, &&\qquad \forall i=1,\cdots ,n \\ &\pi +\theta _i\leqslant -\sigma _ix_i, &&\qquad \forall i=1,\cdots ,n \\ &\pi ,\lambda _i,\theta _i\leqslant 0, &&\qquad \forall i=1,\cdots ,n \end{aligned}$
于是，(3)就可以写成

$\begin{aligned} &t\leqslant \underset{\pi ,\lambda ,\theta}{\max}\,\,\sum_{i=1}^n{\mu _ix_i}+\pi \Gamma +\sum_{i=1}^n{\lambda _i}+\sum_{i=1}^n{\theta _i}, \qquad (4) \\ &\pi +\lambda _i\leqslant \sigma _ix_i, \qquad \qquad \forall i=1,\cdots ,n \qquad (5) \\ &\pi +\theta _i\leqslant -\sigma _ix_i, \qquad \qquad \forall i=1,\cdots ,n \qquad (6) \\ &\pi ,\lambda _i,\theta _i\leqslant 0, \qquad \qquad \forall i=1,\cdots ,n \qquad (7) \end{aligned}$
当然，这里也需要验证强对偶成立(我们在上面已经验证过了)。对于(4)，我们可以直接将  拿掉，变成

$\begin{aligned} &t\leqslant \sum_{i=1}^n{\mu ix_i}+ \pi \Gamma +\sum{i=1}^n{\lambda i}+\sum{i=1}^n{\theta _i}, \qquad (8) \end{aligned}$

因为，  右端项的任意一个值，就能推出 $t\leqslant \underset{\pi ,\lambda ,\theta}{\max}\,\, \sum_{i=1}^n{\mu ix_i}+\pi \Gamma +\sum{i=1}^n{\lambda i}+\sum{i=1}^n{\theta _i}$ 。
基于上述讨论，我们可以整理出最终的等价形式。

利用epigraph转化后的最终形式：

$\begin{aligned} \underset{x\in \mathbb{R}^n}{\max} \quad &t && \\ > &t\leqslant \sum_{i=1}^n{\mu _ix_i}+ \pi \Gamma +\sum_{i=1}^n{\lambda _i}+\sum_{i=1}^n{\theta _i}, && > \\ &\pi +\lambda _i\leqslant \sigma _ix_i, &&\forall i=1,\cdots ,n \\ &\pi +\theta _i\leqslant -\sigma _ix_i, &&\forall i=1,\cdots ,n \\ &\pi ,\lambda _i,\theta _i\leqslant 0, &&\forall i=1,\cdots ,n \\ &\sum_{i=1}^n{x_i}=100\% && \\ &x_i\geqslant 0,& &\forall i=1,\cdots ,n && \\ & t \text{ free} \end{aligned}$
我们发现，实质上和直接使用对偶方法获得的结果是等价的。最后的差别只是多了一个辅助变量$t$而已。为了方便读者对照，我们将直接使用对偶得到的结果也放在这里。
结合内层的对偶问题，原问题可以转化为单层模型，如下

直接对偶得到的结果：

$\begin{aligned} \underset{x,\pi ,\lambda ,\theta}{\max}\qquad &\sum_{i=1}^n{\mu _ix_i}+\pi \Gamma +\sum_{i=1}^n{\lambda _i}+\sum_{i=1}^n{\theta _i} && \\ &\sum_{i=1}^n{x_i}=1 && \\ &\pi +\lambda _i\leqslant \sigma _ix_i, &&\qquad \forall i=1,\cdots ,n \\ &\pi +\theta _i\leqslant -\sigma _ix_i, &&\qquad \forall i=1,\cdots ,n \\ &\pi ,\lambda _i,\theta _i\leqslant 0, &&\qquad \forall i=1,\cdots ,n \\ &x_i\geqslant 0,&& \qquad \forall i=1,\cdots ,n \end{aligned}$
当然了，这个结果不用代码验证都能看出来是等价的。但是为了亲眼看一看，我们就多此一举地再来验证一下了，毕竟都到这个份儿上了。当然，也是顺便凑点字数。
3.2 Python调用gurobi求解reformulation后的模型

完整代码如下：
# author: Xinglu Liu

from gurobipy import *
import numpy as np

# stock number
stock_num = 150
Gamma = 4
# input parameters
mu = [0] * stock_num
sigma = [0] * stock_num
reward = [[0, 0]] * stock_num

# initialize the parameters
for i in range(stock_num):
mu = 0.15 + (i + 1) * (0.05 / 150)
sigma = (0.05 / 450) * math.sqrt(2 * stock_num * (i + 1) * (stock_num + 1))
reward[0] = mu - sigma
reward[1] = mu + sigma
mu = np.array(mu)
sigma = np.array(sigma)
reward = np.array(reward)

model = Model(&#39;Robust portfolio Dual&#39;)
x = {}
lam = {}
theta = {}
pi = {}
pi[0] = model.addVar(lb=-GRB.INFINITY, ub=0, obj=0, vtype=GRB.CONTINUOUS, name=&#39;pi&#39;)
t = model.addVar(lb=-GRB.INFINITY, ub=GRB.INFINITY, obj=1, vtype=GRB.CONTINUOUS, name=&#39;t&#39;)
for i in range(stock_num):
x = model.addVar(lb=0, ub=1, obj=0, vtype=GRB.CONTINUOUS, name=&#39;x_&#39; + str(i + 1))
lam = model.addVar(lb=-GRB.INFINITY, ub=0, obj=0, vtype=GRB.CONTINUOUS, name=&#39;lam_&#39; + str(i + 1))
theta = model.addVar(lb=-GRB.INFINITY, ub=0, obj=0, vtype=GRB.CONTINUOUS, name=&#39;theta_&#39; + str(i + 1))

# set the objective function
model.setObjective(t, GRB.MAXIMIZE)

# constraints 1 : dual constraint
rhs = LinExpr(0)
rhs.addTerms(Gamma, pi[0])
for i in range(stock_num):
rhs.addTerms(mu, x)
rhs.addTerms(1, lam)
rhs.addTerms(1, theta)
model.addConstr(t <= rhs, name=&#39;dual constraint&#39;)

# constraints 1 : invest budget
lhs = LinExpr(0)
for i in range(stock_num):
lhs.addTerms(1, x)
model.addConstr(lhs == 1, name=&#39;invest budget&#39;)

# constraints 2: dual constraint
for i in range(stock_num):
lhs = LinExpr(0)
rhs = LinExpr(0)
lhs.addTerms(1, pi[0])
lhs.addTerms(1, lam)
rhs.addTerms(sigma, x)
model.addConstr(lhs <= rhs, name=&#39;dual_cons1_&#39; + str(i + 1))

# constraints 3: dual constraint
for i in range(stock_num):
lhs = LinExpr(0)
rhs = LinExpr(0)
lhs.addTerms(1, pi[0])
lhs.addTerms(1, lam)
rhs.addTerms(-sigma, x)
model.addConstr(lhs <= rhs, name=&#39;dual_cons2_&#39; + str(i + 1))

model.write(&#39;model.lp&#39;)

model.optimize()

# print out the solution
print(&#39;-----------------------------------------------&#39;)
print(&#39;----------- optimal solution -------------&#39;)
print(&#39;-----------------------------------------------&#39;)
print(&#39;objective : &#39;, round(model.ObjVal, 10))
print(&#39;t =  : &#39;, t.x)
for key in x.keys():
if (x[key].x > 0):
      print(x[key].varName, &#39;\t = &#39;, round(x[key].x, 10))

x_sol = []
for key in x.keys():
x_sol.append(round(x[key].x, 4))

import matplotlib.pyplot as plt

fig = plt.figure(figsize=(12, 8))

plt.plot(range(1, stock_num + 1), x_sol,
      linewidth=2, color=&#39;b&#39;)
plt.xlabel(&#39;Stocks&#39;)
plt.ylabel(&#39;Fraction of investment&#39;)
plt.show()求解结果如下
Iteration Objective    Primal Inf. Dual Inf.    Time
   0 2.0000000e-01 2.896358e-01 0.000000e+00    0s
   79 1.7378554e-01 0.000000e+00 0.000000e+00    0s

Solved in 79 iterations and 0.00 seconds
Optimal objective  1.737855434e-01
-----------------------------------------------
----------- optimal solution -------------
-----------------------------------------------
objective :  0.1737855434
t =  :  0.17378554337248034
x_72    =  0.0154576484
x_73    =  0.015351409
x_74    =  0.0152473305
x_75    =  0.0151453405
x_76    =  0.0150453702
x_77    =  0.0149473537
x_78    =  0.0148512282
x_79    =  0.0147569338
x_80    =  0.0146644129
x_81    =  0.0145736107
x_82    =  0.0144844746
x_83    =  0.0143969543
x_84    =  0.0143110016
x_85    =  0.0142265702
x_86    =  0.0141436157
x_87    =  0.0140620956
x_88    =  0.0139819691
x_89    =  0.0139031968
x_90    =  0.0138257411
x_91    =  0.0137495656
x_92    =  0.0136746355
x_93    =  0.0136009173
x_94    =  0.0135283785
x_95    =  0.0134569882
x_96    =  0.0133867162
x_97    =  0.0133175338
x_98    =  0.0132494129
x_99    =  0.0131823269
x_100 =  0.0131162496
x_101 =  0.0130511562
x_102 =  0.0129870223
x_103 =  0.0129238248
x_104 =  0.0128615409
x_105 =  0.012800149
x_106 =  0.0127396278
x_107 =  0.0126799571
x_108 =  0.0126211171
x_109 =  0.0125630887
x_110 =  0.0125058533
x_111 =  0.0124493932
x_112 =  0.0123936909
x_113 =  0.0123387297
x_114 =  0.0122844933
x_115 =  0.0122309658
x_116 =  0.012178132
x_117 =  0.0121259771
x_118 =  0.0120744865
x_119 =  0.0120236463
x_120 =  0.011973443
x_121 =  0.0119238633
x_122 =  0.0118748944
x_123 =  0.011826524
x_124 =  0.0117787399
x_125 =  0.0117315303
x_126 =  0.0116848839
x_127 =  0.0116387895
x_128 =  0.0115932363
x_129 =  0.0115482139
x_130 =  0.0115037119
x_131 =  0.0114597205
x_132 =  0.0114162299
x_133 =  0.0113732308
x_134 =  0.0113307139
x_135 =  0.0112886703
x_136 =  0.0112470913
x_137 =  0.0112059683
x_138 =  0.0111652931
x_139 =  0.0111250577
x_140 =  0.0110852542
x_141 =  0.0110458748
x_142 =  0.0110069122
x_143 =  0.0109683589
x_144 =  0.010930208
x_145 =  0.0108924524
x_146 =  0.0108550854
x_147 =  0.0108181004
x_148 =  0.0107814908
x_149 =  0.0107452504
x_150 =  0.010709373我们发现，这个解与之前直接调用ROME求解的解，以及直接用对偶得到的结果都是相同的。
4 小结

求解鲁棒优化模型的两大类方法：

求解鲁棒优化一般有两大类方法：
1. 将鲁棒优化模型重构（reformulate）为单层(single-level)段线性规划模型，然后直接调用通用的求解器求解；
2. 直接使用鲁棒优化求解器(ROME, RSOME等)建模求解。

Reformulate鲁棒优化模型的几种方法：

对偶reformulation，但是需要注意，这种方法应用的前提是内层模型一定需要是线性规划，以及reformulate之后一定要验证强对偶成立。

上镜图reformulation，这种方法也涉及到对偶。最后得到的模型与第一种方法的结果是等价的。

其他方法：这里暂不介绍，等到后面有涉及到再做详细介绍。

最近的两篇推文：

【鲁棒优化笔记】基于ROME编程入门鲁棒优化：以一个例子引入（一）

【鲁棒优化笔记】以Coding入门鲁棒优化：以一个例子引入（二）

是笔者学习鲁棒优化过程中精心整理的学习笔记，理论+模型+推导+代码都是一整套的，目前鲁棒优化的基础入门资料相对较少，希望这些资料可以对读者朋友们有所帮助。
由于笔者水平有限，推文中难免有错误，如果读者朋友们发现有错误，请及时联系我更正，非常感谢。
5 参考文献

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.

Gorissen B L, Yankolu , den Hertog D. A practical guide to robust optimization[J]. Omega, 2015, 53: 124-137.

<hr/>运小筹公众号致力于分享运筹优化(LP、MIP、NLP、随机规划、鲁棒优化)、凸优化、强化学习等研究领域的内容以及涉及到的算法的代码实现。编程语言和工具包括Java、Python、Matlab、CPLEX、Gurobi、SCIP 等。欢迎广大同行投稿！欢迎大家关注我们的公众号“运小筹”以及粉丝群！
如果群满加不进去，可以加本文作者微信，然后备注姓名+机构+专业，然后拉您入群。
QQ群：
QQ群里有我们共享的一些书籍、论文、算例等资料，欢迎大家加入。
<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/> 作者：刘兴禄，清华大学，清华-伯克利深圳学院，博士在读
编辑：张瑞三，四川大学，硕士在读，邮箱：493810171@qq.com、sum3rebirth@gmail.com

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

鲁棒优化 | 以Coding入门鲁棒优化：以一个例子引入（二）

本帖子中包含更多资源