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

投资组合的例子
鲁棒优化模型的reformulation: 利用对偶进行reformulation
- 利用对偶进行reformulation
- Python调用gurobi求解对偶reformulation后的模型
鲁棒优化模型的reformulation: 利用上镜图reformulation （个人认为略有问题）
- 利用上镜图reformulation
- Python调用gurobi求解reformulation后的模型
小结
参考文献

** 作者：刘兴禄，清华大学，清华-伯克利深圳学院，博士在读**

投资组合的例子

上一篇我们介绍了一个鲁棒优化求解投资组合的例子，文章链接为：
https://blog.csdn.net/HsinglukLiu/article/details/122769519

其中，该例子的Robust Counterpart模型如下：

max⁡x∈Rnmin⁡z∈Z∑i=1n(μi+σizi)xi∑i=1nxi=100%xi⩾0,∀i=1,2,⋯,n\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} x∈Rnmaxz∈Zmini=1∑n(μi+σizi)xii=1∑nxi=100%xi⩾0,∀i=1,2,⋯,n
其中，Z={z∣∑i=1n∣zi∣⩽Γ,zi∈[−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}.

该例中，我们采用了预算不确定集(Budgeted uncertainty set)。

我们也在上一篇文章中提到，求解鲁棒优化模型的两大类方法：

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

将鲁棒优化模型重构（reformulate）为一阶段线性规划，然后直接调用求解器求解；

直接使用鲁棒优化求解器(ROME, RSOME等)建模求解。

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

reformulation又分为很多种，本文仅仅介绍两种：

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

鲁棒优化模型的reformulation: 利用对偶进行reformulation

利用对偶进行reformulation

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

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

在这种设定下，庄家和玩家的目标函数表达式是一致的，都是∑i=1n(μi+σizi)xi\sum_{i=1}^n{\left( \mu _i+\sigma _iz_i \right) x_i}∑i=1n(μi+σizi)xi，但是二者的目的(或者说，优化方向)是相反的，庄家想要min⁡z∑i=1n(μi+σizi)xi\underset{z}{\min} \sum_{i=1}^n{\left( \mu _i+\sigma _iz_i \right) x_i}zmin∑i=1n(μi+σizi)xi，但是玩家想要max⁡x∑i=1n(μi+σizi)xi\underset{x}{\max} \sum_{i=1}^n{\left( \mu _i+\sigma _iz_i \right) x_i}xmax∑i=1n(μi+σizi)xi。

此外，其实这里有一个隐含的假设：庄家的决策和玩家的决策没有先后顺序。都是基于预判对方的决策，直接做出了自己的决策，后续无论庄家怎么变，玩家都不再改变自己的决策。而不是：庄家真正的先做决策，玩家再去根据对方的决策调整自己的策略。

小拓展：
还有一种bi-level的形式，是两个主体的目标函数或优化方向不一致(也可以目标函数和优化方向都不一致)。这种形式是更为一般的bi-level决策，是典型的的博弈问题。这种情况一般有一个leader和一个follower组成。比如

leader的决策： max⁡∑cexes.t.leader的约束\begin{aligned}\max &\sum c_e x_e\\ s.t. &\text{ leader的约束}\end{aligned}maxs.t.∑cexe leader的约束

follower的决策： max⁡∑deyes.t.follower的约束\begin{aligned}\max &\sum d_e y_e\\ s.t. &\text{ follower的约束}\end{aligned}maxs.t.∑deye follower的约束
比如，共享出行中，平台想最大化整个系统的利润，但是对于每个个体而言，他们只想最大化自己的收入。这两个目标函数是不同的，但是所做的决策之间却有一部分耦合。
关于这一部分，本文只简单提及，不再做详细介绍。

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

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

注意：本文介绍的这种基于对偶的reformulation方法有两个非常强的前提条件：

inner level的模型需要是线性模型，不能有binary和integer的变量以及非线性项。原因是:
-a) MIP是不能以这种方法直接对偶的（但是可以采用其他方法处理，不过也并不是完全等价转化）。
-b) 有非线性项的情况下，如果决策变量都是连续的，可以用KKT条件等价转化。

强对偶是成立的。只有强对偶是成立的，这种转化才是等价的。验证强对偶成立与否的方法很简单，对偶完成以后，看看对偶问题是否存在可行解，如果存在至少一个有解的可行解，则强对偶必然成立。这也是基于对偶定理得来的，即：如果一个问题有可行解，并且目标函数是有界的(因此也有最优解)，那么另一个问题也有可行解和有界的目标函数以及最优解，此时强对偶理论和弱对偶理论都是成立的。

接下来，我们来进行reformulation。

首先我们将模型中的绝对值进行线性化(即∑i=1n∣zi∣⩽Γ\sum_{i=1}^n{|z_i|}\leqslant \Gamma∑i=1n∣zi∣⩽Γ)。

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

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

拓展
之前的推文中，有小伙伴纠结这个转化是否等价。一个肯定的回答是：完全等价。
这里我给出证明：

命题：
z+=max⁡{0,z}(1)z−=max⁡{0,−z}(2)\begin{aligned} &z^+=\max \left\{ 0,z \right\} \,\, \left( 1 \right) \\ &z^-=\max \left\{ 0,-z \right\} \,\, \left( 2 \right) \end{aligned} z+=max{0,z}(1)z−=max{0,−z}(2)
和
z=z+−z−(3)∣z∣=z++z−(4)\begin{aligned} &z=z^+-z^-\,\, \left( 3 \right) \\ &|z|=z^++z^-\,\,\left( 4 \right) \end{aligned} z=z+−z−(3)∣z∣=z++z−(4)
互为充要条件。

证明：
充分性：根据(1)(2)，我们有
z=z+−z−(3)∣z∣=z++z−(4)\begin{aligned} &z=z^+-z^-\,\, \left( 3 \right) \\ &|z|=z^++z^-\,\,\left( 4 \right) \end{aligned} z=z+−z−(3)∣z∣=z++z−(4)
充分性得证
必要性：根据(3)(4)我们有
(3)+(4)=z+∣z∣=2z+→z+=z+∣z∣2(3)−(4)=z−∣z∣=−2z−→z−=−z−∣z∣2情况一：ifz⩾0,z+=z+∣z∣2=z+z2=z=max⁡{0,z}z−=−z−∣z∣2=−z−z2=0=max⁡{0,−z}情况二：ifz⩽0,z+=z+∣z∣2=z−z2=0=max⁡{0,z}z−=−z−∣z∣2=−z−(−z)2=−2z2=−z=max⁡{0,−z}\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} (3)+(4)=z+∣z∣=2z+→z+=2z+∣z∣(3)−(4)=z−∣z∣=−2z−→z−=−2z−∣z∣情况一：ifz⩾0,z+=2z+∣z∣=2z+z=z=max{0,z}z−=−2z−∣z∣=−2z−z=0=max{0,−z}情况二：ifz⩽0,z+=2z+∣z∣=2z−z=0=max{0,z}z−=−2z−∣z∣=−2z−(−z)=−22z=−z=max{0,−z}
因此，必要性得证。
综上，(1)(2)和(3)(4)互为充要条件，以上转化完全等价。

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

max⁡x∈Rnmin⁡zi+,zi−∑i=1n[μi+σi(zi+−zi−)]xi∑i=1nxi=100%∑i=1n(zi++zi−)⩽Γxi⩾0,∀i=1,⋯,nzi+,zi−∈[0,1],∀i=1,⋯,n\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} x∈Rnmaxzi+,zi−mini=1∑n[μi+σi(zi+−zi−)]xii=1∑nxi=100%i=1∑n(zi++zi−)⩽Γxi⩾0,zi+,zi−∈[0,1],∀i=1,⋯,n∀i=1,⋯,n

注意：由于在鲁棒优化中，决策者往往并不关心worst case究竟长什么样子，只在意自己的决策应该怎么做，因此决策者也就不关心zzz的取值。

下面我们来进行对偶。我们将内层模型由min⁡\minmin对偶成max⁡\maxmax问题，这样就可以把双层(bi-level)模型转化为单层(single-level)模型。这里注意，对于内层模型(决策worst-case)而言，决策变量就只有zi+,zi−z_i^{+}, z_i^{-}zi+,zi−，而第二层的决策变量xxx就要被视为固定值(固定值)，也就是视为参数。

因此，内层模型就变化为(因为xxx视为固定值，故相关约束全部可以暂时删除)

min⁡zi+,zi−C+∑i=1nσixizi+−∑i=1nσixizi−∑i=1n(zi++zi−)⩽Γ→πzi+⩽1,∀i=1,⋯,n→λizi−⩽1,∀i=1,⋯,n→θizi+,zi−⩾0∀i=1,⋯,n\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} zi+,zi−minC+i=1∑nσixizi+−i=1∑nσixizi−i=1∑n(zi++zi−)⩽Γzi+⩽1,∀i=1,⋯,nzi−⩽1,∀i=1,⋯,nzi+,zi−⩾0∀i=1,⋯,n→π→λi→θi

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

max⁡π,λ,θπΓ+∑i=1nλi+∑i=1nθiπ+λi⩽σixi,∀i=1,⋯,nπ+θi⩽−σixi,∀i=1,⋯,nπ,λi,θi⩽0,∀i=1,⋯,n\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} π,λ,θmaxπΓ+i=1∑nλi+i=1∑nθiπ+λi⩽σixi,π+θi⩽−σixi,π,λi,θi⩽0,∀i=1,⋯,n∀i=1,⋯,n∀i=1,⋯,n

这里，我们来验证模型的强对偶是否成立。我们发现当所有π,λ,θ\pi, \lambda, \thetaπ,λ,θ均为0，是对偶模型的一个可行解，且目标函数值为0。也就是说，所有变量均为0，是一个有界的可行解。因此，强对偶成立。

结合内层的对偶问题，原问题可以转化为单层模型，如下

max⁡x,π,λ,θ∑i=1nμixi+πΓ+∑i=1nλi+∑i=1nθi∑i=1nxi=1π+λi⩽σixi,∀i=1,⋯,nπ+θi⩽−σixi,∀i=1,⋯,nπ,λi,θi⩽0,∀i=1,⋯,nxi⩾0,∀i=1,⋯,n\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} x,π,λ,θmaxi=1∑nμixi+πΓ+i=1∑nλi+i=1∑nθii=1∑nxi=1π+λi⩽σixi,π+θi⩽−σixi,π,λi,θi⩽0,xi⩾0,∀i=1,⋯,n∀i=1,⋯,n∀i=1,⋯,n∀i=1,⋯,n

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[i] = 0.15 + (i + 1) * (0.05/150)sigma[i] = (0.05/450) * math.sqrt(2 * stock_num * (i + 1) * (stock_num + 1)) reward[i][0] = mu[i] - sigma[i] reward[i][1] = mu[i] + sigma[i]
mu = np.array(mu)
sigma = np.array(sigma)
reward = np.array(reward)  model = Model('Robust portfolio Dual')
x = {}
lam = {}
theta = {}
pi = {}
pi[0] = model.addVar(lb = -GRB.INFINITY, ub = 0, obj = 0, vtype = GRB.CONTINUOUS, name = 'pi')
for i in range(stock_num):x[i] = model.addVar(lb = 0, ub = 1, obj = 0, vtype = GRB.CONTINUOUS, name = 'x_' + str(i+1))lam[i] = model.addVar(lb = -GRB.INFINITY, ub = 0, obj = 0, vtype = GRB.CONTINUOUS, name = 'lam_' + str(i+1))  theta[i] = model.addVar(lb = -GRB.INFINITY, ub = 0, obj = 0, vtype = GRB.CONTINUOUS, name = 'theta_' + str(i+1))    # set the objective function
obj = LinExpr(0)
obj.addTerms(Gamma, pi[0])
for i in range(stock_num):obj.addTerms(mu[i], x[i]) obj.addTerms(1, lam[i])  obj.addTerms(1, theta[i])
model.setObjective(obj, GRB.MAXIMIZE)  # constraints 1 : invest budget
lhs = LinExpr(0)
for i in range(stock_num):lhs.addTerms(1, x[i])
model.addConstr(lhs == 1, name = 'invest budget') # 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[i])  rhs.addTerms(sigma[i], x[i])  model.addConstr(lhs <= rhs, name = 'dual_cons1_' + 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[i])  rhs.addTerms(-sigma[i], x[i])   model.addConstr(lhs <= rhs, name = 'dual_cons2_' + str(i+1))model.write('model.lp') model.optimize() # print out the solution
print('-----------------------------------------------')
print('-----------   optimal solution   -------------')
print('-----------------------------------------------')
print('objective : ', round(model.ObjVal, 10))
for key in x.keys():  if(x[key].x > 0):print(x[key].varName, '\t = ', round(x[key].x, 10)) x_sol = []
for key in x.keys():  x_sol.append(round(x[key].x, 4))import matplotlib.pyplot as pltfig = plt.figure(figsize=(12,8))plt.plot(range(1, stock_num+1), x_sol,linewidth=2, color='b')
plt.xlabel('Stocks')
plt.ylabel('Fraction of investment')
plt.show()

求解结果如下：

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 nonzerosIteration    Objective       Primal Inf.    Dual Inf.      Time0    2.0000000e-01   2.896358e-01   0.000000e+00      0s79    1.7378554e-01   0.000000e+00   0.000000e+00      0sSolved 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是正确的。开心~

鲁棒优化模型的reformulation: 利用上镜图reformulation （个人认为略有问题）

利用上镜图reformulation的方法，我个人认为有些小问题，这些也仅代表我自己的探索，希望跟大家探讨。求看出问题的读者朋友多多指教。

利用上镜图reformulation

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

引入辅助决策变量ttt

则上述模型可以转化为

max⁡x∈Rntt⩽∑i=1n(μi+σizi)xi∑i=1nxi=100%∑i=1n∣zi∣⩽Γxi⩾0,∀i=1,⋯,nzi∈[−1,1],∀i=1,⋯,n\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} && \\ &\sum_{i=1}^n{x_i}=100\% && \\ &\sum_{i=1}^n{|z_i|}\leqslant \Gamma && \\ &x_i\geqslant 0,& &\qquad \forall i=1,\cdots ,n \\ &z_i\in \left[ -1,1 \right] ,&&\qquad \forall i=1,\cdots ,n \end{aligned}x∈Rnmaxtt⩽i=1∑n(μi+σizi)xii=1∑nxi=100%i=1∑n∣zi∣⩽Γxi⩾0,zi∈[−1,1],∀i=1,⋯,n∀i=1,⋯,n

其中，约束t⩽∑i=1n(μi+σizi)xit\leqslant \sum_{i=1}^n{\left( \mu _i+\sigma _iz_i \right) x_i}t⩽∑i=1n(μi+σizi)xi就等价于min⁡z∈Z∑i=1n(μi+σizi)xi\underset{z\in \mathcal{Z}}{\min} \sum_{i=1}^n{\left( \mu _i+\sigma _iz_i \right) x_i}z∈Zmin∑i=1n(μi+σizi)xi。

讨论：利用上镜图reformulation是否是等价转化？
这里，我个人认为转化不等价。在bi-level的模型中，zzz和xxx相当于是两个主体所做的决策，并且分属不同level，他们之间的关系是通过max⁡min⁡\max \minmaxmin和其他约束进行耦合的。但是基于上镜图的reformulation，通过引入辅助决策变量ttt，将zzz和xxx放在了同一个level，虽然二者之间的约束层面的关系依然存在，但是从目标函数的对抗/博弈层面的耦合关系却被消除了。因此这种转化并不等价。下面我们也用gurobi验证了这一点，这种reformulation并不是等价转化。

Python调用gurobi求解reformulation后的模型

由于基于上镜图转化后的模型是非线性规划，模型两个非线性项，即

t⩽∑i=1n(μi+σizi)xi∑i=1n∣zi∣⩽Γ\begin{aligned} &t\leqslant \sum_{i=1}^n{\left( \mu _i+\sigma _iz_i \right) x_i} && \\ &\sum_{i=1}^n{|z_i|}\leqslant \Gamma && \end{aligned}t⩽i=1∑n(μi+σizi)xii=1∑n∣zi∣⩽Γ

其中t⩽∑i=1n(μi+σizi)xit\leqslant \sum_{i=1}^n{\left( \mu _i+\sigma _iz_i \right) x_i}t⩽∑i=1n(μi+σizi)xi包含bi-linear项zixiz_ix_izixi，这两都是连续变量，无法等价线性化(两个连续变量相乘，不存在等价的线性化方法，只有近似方法)，但是gurobi可以求解此类非线性问题。具体方法如下：

使用gurobi的QuadExpr对象以及addTerms(coef, var1, var2)来完成
求解的时候，设置模型的参数NonConvex=2，即model.setParam("NonConvex", 2)

gurobi还可以求解∑i=1n∣zi∣⩽Γ\sum_{i=1}^n{|z_i|}\leqslant \Gamma∑i=1n∣zi∣⩽Γ这种包含绝对值的约束。需要使用General Constraint Helper Functions或者addGenConstrAbs来完成，语法为

abs_函数，例如abs_z[i] == abs_(z[i])
addGenConstrAbs(b, a), 表示添加约束b=∣a∣b=|a|b=∣a∣.

完整代码如下：

# author: Xinglu Liufrom 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[i] = round(0.15 + (i + 1) * (0.05/150), 4)sigma[i] = round( (0.05/450) * math.sqrt(2 * stock_num * (i + 1) * (stock_num + 1)), 4)reward[i][0] = round(mu[i] - sigma[i], 4)reward[i][1] = round(mu[i] + sigma[i], 4)
mu = np.array(mu)
sigma = np.array(sigma)
reward = np.array(reward)  model = Model('Robust portfolio Model')
x = {}
z = {}
abs_z = {}
t = model.addVar(lb = 0, ub = GRB.INFINITY, obj = 0, vtype = GRB.CONTINUOUS, name = 't')
for i in range(stock_num):x[i] = model.addVar(lb = 0, ub = 1, obj = 0, vtype = GRB.CONTINUOUS, name = 'x_' + str(i+1))z[i] = model.addVar(lb = -1, ub = 1, obj = 0, vtype = GRB.CONTINUOUS, name = 'z_' + str(i+1)) abs_z[i] = model.addVar(lb = 0, ub = 1, obj = 0, vtype = GRB.CONTINUOUS, name = 'abs_z_' + str(i+1))   # set the objective function
model.setObjective(t, GRB.MAXIMIZE)# constraints 1: epigraph constraint
lhs = LinExpr(0)
rhs = QuadExpr(0)
lhs.addTerms(1, t)
for i in range(stock_num):rhs.addTerms(mu[i], x[i])  rhs.addTerms(sigma[i], z[i], x[i])
model.addQConstr(lhs <= rhs, name = 'epigraph_constraint')  # constraints 2 : invest budget
lhs = LinExpr(0)
for i in range(stock_num):lhs.addTerms(1, x[i])
model.addConstr(lhs == 1, name = 'invest budget')   # constraints 3 : risk budget
lhs = LinExpr(0)
for i in range(stock_num): lhs.addTerms(1, abs_z[i])
model.addConstr(lhs <= Gamma, name = 'risk budget') # constraints 4 : abs constraints
for i in range(stock_num): # x5 = abs(x1)model.addGenConstrAbs(abs_z[i], z[i], "absconstr_" + str(i+1))  # overloaded form
#     model.addConstr(abs_z[i] == abs_(z[i]), name="absconstr_" + str(i+1)) model.write('model.lp')
model.setParam("NonConvex", 2)
model.optimize() # print out the solution
print('-----------------------------------------------')
print('-----------   optimal solution   -------------')
print('-----------------------------------------------')
print('Obj: {}'.format(model.ObjVal))
for key in x.keys():  print(x[key].varName, round(x[key].x, 2), '\t,', z[key].varName, round(z[key].x, 2), '\t,', abs_z[key].varName, round(abs_z[key].x, 2))

求解结果如下

Root relaxation: objective 4.896000e-01, 531 iterations, 0.00 secondsNodes    |    Current Node    |     Objective Bounds      |     WorkExpl Unexpl |  Obj  Depth IntInf | Incumbent    BestBd   Gap | It/Node Time*    0     0               0       0.4896000    0.48960  0.00%     -    0sExplored 0 nodes (531 simplex iterations) in 0.01 seconds
Thread count was 16 (of 16 available processors)Solution count 1: 0.4896
No other solutions better than 0.4896Optimal solution found (tolerance 1.00e-04)
Best objective 4.896000000000e-01, best bound 4.896000000000e-01, gap 0.0000%
-----------------------------------------------
-----------   optimal solution   -------------
-----------------------------------------------
Obj: 0.48960000000000004
x_101 = 0.0    , z_101 = -1.0     , abs_z_101 = 1.0
x_110 = 0.0    , z_110 = -1.0     , abs_z_110 = 1.0
x_149 = 0.0    , z_149 = -1.0     , abs_z_149 = 1.0
x_150 = 1.0    , z_150 = 1.0  , abs_z_150 = 1.0

我们发现，这个解与之前直接调用ROME求解的并不相同。

在这个解中：

我们需要把100%的资产全部投放到第150号股票上，并且最终的收益为0.4896。
收益的波动量，在150号股票上达到了最大，相应的z150=1z_{150}=1z150=1，而150号股票的收益变动范围为[−0.0896,0.4896][-0.0896, 0.4896][−0.0896,0.4896]。

非常明显，这并不是worst case，相反，这应该是是best case才对。哈哈哈~

这种上镜图的reformulation方法虽然简单，但是由于reformulation后模型有时候会不等价，因此，还是建议用基于对偶的方法。

小结

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

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

将鲁棒优化模型重构（reformulate）为单层(single-level)段线性规划模型，然后直接调用通用的求解器求解；

直接使用鲁棒优化求解器(ROME, RSOME等)建模求解。

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

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

上镜图reformulation，但是该种方法需要注意，reformulation之后的非线性项的问题。以及，需要关注reformulation之后模型是否等价。(个人认为本文介绍的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.