【鲁棒优化笔记】以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 ∈ R n min z ∈ Z ∑ i = 1 n ( μ i + σ i z i ) x i ∑ i = 1 n x i = 100 % x i ⩾ 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 = 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}.
该例中，我们采用了预算不确定集(Budgeted uncertainty set)。
我们也在上一篇文章中提到，求解鲁棒优化模型的两大类方法：
求解鲁棒优化一般有两大类方法：
- 将鲁棒优化模型重构（reformulate）为一阶段线性规划，然后直接调用求解器求解；
- 直接使用鲁棒优化求解器(ROME, RSOME等)建模求解。
其中，第2类方法，我们已经在器上一篇文章中做了详细解释，本文我们来介绍第一类求解方法：reformulation。
reformulation又分为很多种，本文仅仅介绍两种：
- 利用对偶进行reformulation
- 利用上镜图(
epigraph)进行reformulation
鲁棒优化模型的reformulation: 利用对偶进行reformulation
利用对偶进行reformulation
上述模型其实是一个bi-level的模型，也就是 max x min z \underset{x}{\max} \,\,\underset{z}{\min} xmaxzmin。上述鲁棒优化模型可以在一定程度上理解为是一个赌场博弈问题(只是为了方便理解，并不一定严谨)：
inner level或lower level的决策：不确定变量 z z z，可以看做是赌场的庄家，庄家不想让玩家赚钱，所以他控制收益 r ~ \tilde{r} r~波动(也就是控制风险)，尽可能让你赚的少(可以理解为给你制造worst case)outer level或upper level的决策：投资决策 x x x。作为玩家，我们需要跟庄家博弈，我们要最大化收益，使得在worst case下的收益最大化。
在这种设定下，庄家和玩家的目标函数表达式是一致的，都是 ∑ i = 1 n ( μ i + σ i z i ) x i \sum_{i=1}^n{\left( \mu _i+\sigma _iz_i \right) x_i} ∑i=1n(μi+σizi)xi，但是二者的目的(或者说，优化方向)是相反的，庄家想要 min z ∑ i = 1 n ( μ i + σ i z i ) x i \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 = 1 n ( μ i + σ i z i ) x i \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 ∑ c e x e s . t . leader的约束 \begin{aligned}\max &\sum c_e x_e\\ s.t. &\text{ leader的约束}\end{aligned} maxs.t.∑cexe leader的约束
- follower的决策： max ∑ d e y e s . 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 \, \min maxmin，对inner level模型取对偶，将bi-level模型转化成为 max max \max\,\max maxmax，进而等价为一阶段 max \max max问题进行求解
注意：本文介绍的这种基于对偶的reformulation方法有两个非常强的前提条件：
inner level的模型需要是线性模型，不能有binary和integer的变量以及非线性项。原因是:
-a)MIP是不能以这种方法直接对偶的（但是可以采用其他方法处理，不过也并不是完全等价转化）。
-b)有非线性项的情况下，如果决策变量都是连续的，可以用KKT条件等价转化。强对偶是成立的。只有强对偶是成立的，这种转化才是等价的。验证强对偶成立与否的方法很简单，对偶完成以后，看看对偶问题是否存在可行解，如果存在至少一个有解的可行解，则强对偶必然成立。这也是基于对偶定理得来的，即：如果一个问题有可行解，并且目标函数是有界的(因此也有最优解)，那么另一个问题也有可行解和有界的目标函数以及最优解，此时强对偶理论和弱对偶理论都是成立的。
接下来，我们来进行reformulation。
首先我们将模型中的绝对值进行线性化(即 ∑ i = 1 n ∣ z i ∣ ⩽ Γ \sum_{i=1}^n{|z_i|}\leqslant \Gamma ∑i=1n∣zi∣⩽Γ)。
我们引入两组辅助变量 z i + , z i − , ∀ i = 1 , 2 , ⋯ , n z_i^{+}, z_i^{-}, \forall i =1,2,\cdots, n zi+,zi−,∀i=1,2,⋯,n。其中
z i + = max { 0 , z i } z i − = max { 0 , − z i } \begin{aligned} &z_i^{+} = \max\{0, z_i\} \\ &z_i^{-} = \max\{0, -z_i\} \end{aligned} zi+=max{0,zi}zi−=max{0,−zi}
因此，就有
z i = z i + − z i − ∣ z i ∣ = z i + + z i − \begin{aligned} &z_i = z_i^{+} - z_i^{-} \\ &|z_i| = z_i^{+} + z_i^{-} \end{aligned} zi=zi+−zi−∣zi∣=zi++zi−
这个比较容易理解，例如 z i = − 0.5 z_i = -0.5 zi=−0.5， 则 z i + = 0 , z i − = 0.5 z_i^{+} = 0, z_i^{-} = 0.5 zi+=0,zi−=0.5， 因此有
z i = z i + − z i − = 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 i ∣ = z i + + z i − = 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 ∣ = 2 z + → z + = z + ∣ z ∣ 2 ( 3 ) − ( 4 ) = z − ∣ z ∣ = − 2 z − → z − = − z − ∣ z ∣ 2 情况一： i f z ⩾ 0 , z + = z + ∣ z ∣ 2 = z + z 2 = z = max { 0 , z } z − = − z − ∣ z ∣ 2 = − z − z 2 = 0 = max { 0 , − z } 情况二： i f z ⩽ 0 , z + = z + ∣ z ∣ 2 = z − z 2 = 0 = max { 0 , z } z − = − z − ∣ z ∣ 2 = − z − ( − z ) 2 = − 2 z 2 = − 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 ∈ R n min z i + , z i − ∑ i = 1 n [ μ i + σ i ( z i + − z i − ) ] x i ∑ i = 1 n x i = 100 % ∑ i = 1 n ( z i + + z i − ) ⩽ Γ x i ⩾ 0 , ∀ i = 1 , ⋯ , n z i + , z i − ∈ [ 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究竟长什么样子，只在意自己的决策应该怎么做，因此决策者也就不关心 z z z的取值。
下面我们来进行对偶。我们将内层模型由 min \min min对偶成 max \max max问题，这样就可以把双层(bi-level)模型转化为单层(single-level)模型。这里注意，对于内层模型(决策worst-case)而言，决策变量就只有 z i + , z i − z_i^{+}, z_i^{-} zi+,zi−，而第二层的决策变量 x x x就要被视为固定值(固定值)，也就是视为参数。
因此，内层模型就变化为(因为 x x x视为固定值，故相关约束全部可以暂时删除)
min z i + , z i − C + ∑ i = 1 n σ i x i z i + − ∑ i = 1 n σ i x i z i − ∑ i = 1 n ( z i + + z i − ) ⩽ Γ → π z i + ⩽ 1 , ∀ i = 1 , ⋯ , n → λ i z i − ⩽ 1 , ∀ i = 1 , ⋯ , n → θ i z i + , z i − ⩾ 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
其中 C C C为常数，且 C = ∑ i = 1 n μ i x i C = \sum_{i=1}^n{\mu _ix_i} C=∑i=1nμixi，且约束(1)、约束(2)和约束(3)的对偶变量分别是 π \pi π
和 λ \lambda λ和 θ \theta θ。因此，内层模型的对偶问题为
max π , λ , θ π Γ + ∑ i = 1 n λ i + ∑ i = 1 n θ i π + λ i ⩽ σ i x i , ∀ i = 1 , ⋯ , n π + θ i ⩽ − σ i x i , ∀ 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 = 1 n μ i x i + π Γ + ∑ i = 1 n λ i + ∑ i = 1 n θ i ∑ i = 1 n x i = 1 π + λ i ⩽ σ i x i , ∀ i = 1 , ⋯ , n π + θ i ⩽ − σ i x i , ∀ i = 1 , ⋯ , n π , λ i , θ i ⩽ 0 , ∀ i = 1 , ⋯ , n x i ⩾ 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
由于上述模型是bi-level的模型，即 max min \max \,\, \min maxmin问题，我们可以通过上镜图的方法将其转化为single level的模型(这种做法也是参考一些文献的做法，参见参考文献7)，具体方法为：
- 引入辅助决策变量 t t t
则上述模型可以转化为
max x ∈ R n t t ⩽ ∑ i = 1 n ( μ i + σ i z i ) x i , ∀ z = ( z 1 , ⋯ , z n ) T ∈ Z ∑ i = 1 n x i = 100 % x i ⩾ 0 , ∀ 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}, && \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} x∈Rnmaxtt⩽i=1∑n(μi+σizi)xi,i=1∑nxi=100%xi⩾0,∀z=(z1,⋯,zn)T∈Z∀i=1,⋯,n
其中，约束 t ⩽ ∑ i = 1 n ( μ i + σ i z i ) x i , ∀ z = ( z 1 , ⋯ , z n ) T ∈ Z 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} t⩽∑i=1n(μi+σizi)xi,∀z=(z1,⋯,zn)T∈Z就等价于 min z ∈ Z ∑ i = 1 n ( μ i + σ i z i ) x i \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。
- 一定要注意约束 t ⩽ ∑ i = 1 n ( μ i + σ i z i ) x i t\leqslant \sum_{i=1}^n{\left( \mu _i+\sigma _iz_i \right) x_i} t⩽∑i=1n(μi+σizi)xi后面的 ∀ z = ( z 1 , ⋯ , z n ) T ∈ Z \forall \mathbf{z}=(z_1, \cdots, z_n)^T \in \mathcal{Z} ∀z=(z1,⋯,zn)T∈Z,这个一定不能少。
上述约束其实也可以写为
t ⩽ min z ∑ i = 1 n ( μ i + σ i z i ) x i , ∀ z = ( z 1 , ⋯ , z n ) T ∈ Z ( 1 ) \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} t⩽zmini=1∑n(μi+σizi)xi,∀z=(z1,⋯,zn)T∈Z(1)
这种写法也上面的写法是等价的。因为， t ⩽ t\leqslant t⩽右端项所有可能的取值，就等价于 t ⩽ t\leqslant t⩽右端项的最小值。
实际上，约束(1)自己就是一个小优化问题，并不是能够直接处理的约束。
接下来，我们的任务就是要将(1)转化成为可处理的形式。
我们需要将其作进一步处理。我们将(1)改写成
t ⩽ min z i ∈ [ − 1 , 1 ] , ∀ i ; ∑ ∣ z i ∣ ⩽ Γ ∑ i = 1 n ( μ i + σ i z i ) x i , ( 2 ) \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} t⩽zi∈[−1,1],∀i;∑∣zi∣⩽Γmini=1∑n(μi+σizi)xi,(2)
当然，也可以简写为
t ⩽ min z ∈ Z ∑ i = 1 n ( μ i + σ i z i ) x i , ( 3 ) \begin{aligned} t\leqslant \underset{z_\in \mathcal{Z}}{\min}\sum_{i=1}^n{\left( \mu _i+\sigma _iz_i \right) x_i}, \qquad (3) \end{aligned} t⩽z∈Zmini=1∑n(μi+σizi)xi,(3)
转化的方法，是将 min z ∈ Z ∑ i = 1 n ( μ i + σ i z i ) x i \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取对偶。方法与上面介绍的直接对偶的方法一致，我们可以直接把结果拿过来。
max π , λ , θ ∑ i = 1 n μ i x i + π Γ + ∑ i = 1 n λ i + ∑ i = 1 n θ i π + λ i ⩽ σ i x i , ∀ i = 1 , ⋯ , n π + θ i ⩽ − σ i x i , ∀ i = 1 , ⋯ , n π , λ i , θ i ⩽ 0 , ∀ i = 1 , ⋯ , n \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} π,λ,θmaxi=1∑nμixi+πΓ+i=1∑nλi+i=1∑nθiπ+λi⩽σixi,π+θi⩽−σixi,π,λi,θi⩽0,∀i=1,⋯,n∀i=1,⋯,n∀i=1,⋯,n
于是，(3)就可以写成
t ⩽ max π , λ , θ ∑ i = 1 n μ i x i + π Γ + ∑ i = 1 n λ i + ∑ i = 1 n θ i , ( 4 ) π + λ i ⩽ σ i x i , ∀ i = 1 , ⋯ , n ( 5 ) π + θ i ⩽ − σ i x i , ∀ i = 1 , ⋯ , n ( 6 ) π , λ i , θ i ⩽ 0 , ∀ i = 1 , ⋯ , n ( 7 ) \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} t⩽π,λ,θmaxi=1∑nμixi+πΓ+i=1∑nλi+i=1∑nθi,(4)π+λi⩽σixi,∀i=1,⋯,n(5)π+θi⩽−σixi,∀i=1,⋯,n(6)π,λi,θi⩽0,∀i=1,⋯,n(7)
当然，这里也需要验证强对偶成立(我们在上面已经验证过了)。对于(4)，我们可以直接将 max \max max拿掉，变成
t ⩽ ∑ i = 1 n μ i x i + π Γ + ∑ i = 1 n λ i + ∑ i = 1 n θ i , ( 8 ) \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⩽i=1∑nμixi+πΓ+i=1∑nλi+i=1∑nθi,(8)
因为， t ⩽ t\leqslant t⩽右端项的任意一个值，就能推出 t ⩽ max π , λ , θ ∑ i = 1 n μ i x i + π Γ + ∑ i = 1 n λ i + ∑ i = 1 n θ i 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} t⩽π,λ,θmax∑i=1nμixi+πΓ+∑i=1nλi+∑i=1nθi。
基于上述讨论，我们可以整理出最终的等价形式。
利用epigraph转化后的最终形式：
max x ∈ R n t t ⩽ ∑ i = 1 n μ i x i + π Γ + ∑ i = 1 n λ i + ∑ i = 1 n θ i , π + λ i ⩽ σ i x i , ∀ i = 1 , ⋯ , n π + θ i ⩽ − σ i x i , ∀ i = 1 , ⋯ , n π , λ i , θ i ⩽ 0 , ∀ i = 1 , ⋯ , n ∑ i = 1 n x i = 100 % x i ⩾ 0 , ∀ i = 1 , ⋯ , n t free \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} x∈Rnmaxtt⩽i=1∑nμixi+πΓ+i=1∑nλi+i=1∑nθi,π+λi⩽σixi,π+θi⩽−σixi,π,λi,θi⩽0,i=1∑nxi=100%xi⩾0,t free∀i=1,⋯,n∀i=1,⋯,n∀i=1,⋯,n∀i=1,⋯,n
我们发现，实质上和直接使用对偶方法获得的结果是等价的。最后的差别只是多了一个辅助变量 t t t而已。为了方便读者对照，我们将直接使用对偶得到的结果也放在这里。
结合内层的对偶问题，原问题可以转化为单层模型，如下
- 直接对偶得到的结果： max x , π , λ , θ ∑ i = 1 n μ i x i + π Γ + ∑ i = 1 n λ i + ∑ i = 1 n θ i ∑ i = 1 n x i = 1 π + λ i ⩽ σ i x i , ∀ i = 1 , ⋯ , n π + θ i ⩽ − σ i x i , ∀ i = 1 , ⋯ , n π , λ i , θ i ⩽ 0 , ∀ i = 1 , ⋯ , n x i ⩾ 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 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] = 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')
t = model.addVar(lb=-GRB.INFINITY, ub=GRB.INFINITY, obj=1, 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))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
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[i], x[i])rhs.addTerms(1, lam[i])rhs.addTerms(1, theta[i])
model.addConstr(t <= rhs, name='dual constraint')# 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))
print('t = : ', t.x)
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()
求解结果如下
Iteration 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.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求解的解，以及直接用对偶得到的结果都是相同的。
小结
求解鲁棒优化模型的两大类方法：
求解鲁棒优化一般有两大类方法：
- 将鲁棒优化
模型重构（reformulate）为单层(single-level)段线性规划模型，然后直接调用通用的求解器求解；直接使用鲁棒优化求解器(ROME, RSOME等)建模求解。
Reformulate鲁棒优化模型的几种方法：
对偶reformulation，但是需要注意，这种方法应用的前提是内层模型一定需要是线性规划，以及reformulate之后一定要验证强对偶成立。上镜图reformulation，这种方法也涉及到对偶。最后得到的模型与第一种方法的结果是等价的。- 其他方法：这里暂不介绍，等到后面有涉及到再做详细介绍。
最近的两篇推文：
- 【鲁棒优化笔记】基于ROME编程入门鲁棒优化：以一个例子引入（一）
- 【鲁棒优化笔记】以Coding入门鲁棒优化：以一个例子引入（二）
是笔者学习鲁棒优化过程中精心整理的学习笔记，理论+模型+推导+代码都是一整套的，目前鲁棒优化的基础入门资料相对较少，希望这些资料可以对读者朋友们有所帮助。
由于笔者水平有限，推文中难免有错误，如果读者朋友们发现有错误，请及时联系我更正，非常感谢。
参考文献
- 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, Yanıkoğlu İ, den Hertog D. A practical guide to robust optimization[J]. Omega, 2015, 53: 124-137.
** 作者：刘兴禄， 清华大学，清华-伯克利深圳学院，博士在读**
