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

已于 2022-02-04 22:28:16 修改
阅读量2.5k 收藏 12
点赞数 3
分类专栏： 最优化求解器(Gurobi+CPLEX+SCIP) 运筹优化理论与实践 文章标签： 鲁棒优化 运筹优化 运筹学 算法 Gurobi+python
于 2022-02-03 01:31:04 首次发布
版权声明：本文为博主原创文章，遵循 CC 4.0 BY-SA 版权协议，转载请附上原文出处链接和本声明。
本文链接：https://blog.csdn.net/HsinglukLiu/article/details/122773390
版权

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

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

投资组合的例子

上一篇我们介绍了一个鲁棒优化求解投资组合的例子，文章链接为：
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∈Rnmax​z∈Zmin​​i=1∑n​(μi​+σi​zi​)xi​i=1∑n​xi​=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)。

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

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

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

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

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

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

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

利用对偶进行reformulation

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

  • 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​+σi​zi​)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​+σi​zi​)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​+σi​zi​)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.​∑ce​xe​ 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.​∑de​ye​ follower的约束​
      比如，共享出行中，平台想最大化整个系统的利润，但是对于每个个体而言，他们只想最大化自己的收入。这两个目标函数是不同的，但是所做的决策之间却有一部分耦合。
      关于这一部分，本文只简单提及，不再做详细介绍。

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

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

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

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

接下来，我们来进行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∈Rnmax​zi+​,zi−​min​​i=1∑n​[μi​+σi​(zi+​−zi−​)]xi​i=1∑n​xi​=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−​min​​C+i=1∑n​σi​xi​zi+​−i=1∑n​σi​xi​zi−​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​μi​xi​，且约束(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​⩽σi​xi​,π+θi​⩽−σi​xi​,π,λ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,π,λ,θmax​​i=1∑n​μi​xi​+πΓ+i=1∑n​λi​+i=1∑n​θi​i=1∑n​xi​=1π+λi​⩽σi​xi​,π+θi​⩽−σi​xi​,π,λ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 plt

fig = 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 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是正确的。开心~

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

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

利用上镜图reformulation

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

  • 引入辅助决策变量 t t t

则上述模型可以转化为

max ⁡ x ∈ R n t t ⩽ ∑ i = 1 n ( μ i + σ i z i ) x i ∑ i = 1 n x i = 100 % ∑ i = 1 n ∣ z i ∣ ⩽ Γ x i ⩾ 0 , ∀ i = 1 , ⋯   , n z i ∈ [ − 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∈Rnmax​​tt⩽i=1∑n​(μi​+σi​zi​)xi​i=1∑n​xi​=100%i=1∑n​∣zi​∣⩽Γxi​⩾0,zi​∈[−1,1],​​∀i=1,⋯,n∀i=1,⋯,n​

其中，约束 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​+σi​zi​)xi​就等价于 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​+σi​zi​)xi​。

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

Python调用gurobi求解reformulation后的模型

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

t ⩽ ∑ i = 1 n ( μ i + σ i z i ) x i ∑ i = 1 n ∣ z i ∣ ⩽ Γ \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​+σi​zi​)xi​i=1∑n​∣zi​∣⩽Γ​​​

其中 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​+σi​zi​)xi​包含bi-linear项 z i x i z_ix_i zi​xi​，这两都是连续变量，无法等价线性化(两个连续变量相乘，不存在等价的线性化方法，只有近似方法)，但是gurobi可以求解此类非线性问题。具体方法如下：

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

gurobi还可以求解 ∑ i = 1 n ∣ z i ∣ ⩽ Γ \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 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] = 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 seconds

    Nodes    |    Current Node    |     Objective Bounds      |     Work
 Expl Unexpl |  Obj  Depth IntInf | Incumbent    BestBd   Gap | It/Node Time

*    0     0               0       0.4896000    0.48960  0.00%     -    0s

Explored 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.4896

Optimal 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号股票上达到了最大，相应的 z 150 = 1 z_{150}=1 z150​=1，而150号股票的收益变动范围为 [ − 0.0896 , 0.4896 ] [-0.0896, 0.4896] [−0.0896,0.4896]。

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

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

小结

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

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

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

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

  1. 对偶reformulation，但是需要注意，这种方法应用的前提是内层模型一定需要是线性规划，以及reformulate之后一定要验证强对偶成立。
  2. 上镜图reformulation，但是该种方法需要注意，reformulation之后的非线性项的问题。以及，需要关注reformulation之后模型是否等价。(个人认为本文介绍的reformulation不等价)。
  3. 其他方法：这里暂不介绍，等到后面有涉及到再做详细介绍。

最近的两篇推文：

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

是笔者学习鲁棒优化过程中精心整理的学习笔记，理论+模型+推导+代码都是一整套的，目前鲁棒优化的基础入门资料相对较少，希望这些资料可以对读者朋友们有所帮助。

由于笔者水平有限，推文中难免有错误，如果读者朋友们发现有错误，请及时联系我更正，非常感谢。

参考文献

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

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

刘兴禄
关注
  • 3
    点赞
  • 踩
  • 12
    收藏
    觉得还不错? 一键收藏
  • 1
    评论
专栏目录
鲁棒优化（一）
qq_25018077的博客
09-11 1万+
初步认识鲁棒优化
ruolubang_鲁棒约束_鲁棒_鲁棒优化_鲁棒优化模型_鲁棒优化cplex
09-11
书写代码，并在matlab中使用cplex求解鲁棒优化模型，考虑了各种模型的约束条件
【鲁棒优化笔记】Comment to: 以Coding入门鲁棒优化：以一个例子引入（二）
HsinglukLiu的博客
02-03 1451
Comment to: 以Coding入门鲁棒优化：以一个例子引入（二）证明绝对值线性化的等价性 ** 作者：刘兴禄， 清华大学，清华-伯克利深圳学院，博士在读** 证明绝对值线性化的等价性 下面我们来证明一下绝对值转化的等价性： 命题： z+=max⁡{0,z}  (1)z−=max⁡{0,−z}  (2) z^+=\max \left\{ 0,z \right\} \,\, \left( 1 \right) \\ z^-=\max \left\{ 0,-z \right\} \,\
优化｜鲁棒优化基础
weixin_53463894的博客
12-08 6815
作者：苏向阳 本文讲解了鲁棒优化的基础内容，内容主要分为3个部分：1.不确定性最优化。2.不确定集。3.对等式转换理论。在末尾附上了一个简单的算例及代码便于读者理解。该文也是【运筹OR帷幄】鲁棒优化电子书系列的文章内容之一。 1.不确定性最优化（Optimization under uncertainty） 1.1 不确定优化方法 实际生活中很多问题都具有不确定性（Uncertainty），随着最优化理论的不断发展和计算机能力的提高，不确定性优化受到了学界前所未有的重视。早在20世纪50年代，B..
【鲁棒优化笔记】基于ROME编程入门鲁棒优化：以一个例子引入(一)
HsinglukLiu的博客
02-02 4869
鲁棒优化|基于ROME编程入门鲁棒优化：以一个例子引入（上）线性规划中的不确定性一个生产的例子Matlab+ROME求解鲁棒优化简单案例：一个投资组合的例子不确定回报率鲁棒优化模型Matlab调用ROME包求解Robust Counterpart参考文献 鲁棒优化是处理数学规划中处理不确定性的重要方法论。这里我们以几个简单的例子来介绍鲁棒优化。同时，我们结合鲁棒优化工具包Robust Optimization Made Easy (ROME)来编程实现，方便读者理解和复现。 线性规划中的不确定性 在生产实际
鲁棒优化入门（二）——基于matlab+yalmip求解鲁棒优化问题
配电网和matlab的博客
07-10 9570
鲁棒优化的含义就是在最恶劣的情况下（不确定变量的取值使目标函数最大），求出满足约束条件，并且能使目标函数最优的决策变量。 yalmip工具箱可以用来求解鲁棒优化问题，但还是有一定局限性的，并不能处理任意形式的不确定集下的鲁棒优化问题，一般来说，当鲁棒优化问题的不确定集合为箱型不确定集、椭球不确定集以及多面体不确定集时，都可以用yalmip工具箱求解（具体细节可参考官方文档）。本文介绍了利用yalmip求鲁棒优化问题的一般方法。......
鲁棒优化入门(5)-Matlab+Yalmip求解鲁棒优化编程实战
08-20
包含以下内容： 系统介绍了如何通过Matlab+Yalmip工具箱求解单阶段鲁棒优化的方法，总共包括三种方法和两个算例： ...5.代码使用matlab+Yalmip求解，第二个鲁棒优化问题中还用到了MATPOWER工具箱。
MATLAB代码：计及条件风险价值的电-气综合能源系统能量备用-分布鲁棒优化
06-02
关键词： wasserstein距离 CVAR条件风险价值 分布鲁棒优化 电-气综合能源 能量-备用调度 完美复现：SCI《Energy and Reserve Dispatch with Distributionally Robust Joint Chance Constraints》。 主要内容：代码...
详细介绍电力系统鲁棒经济调度，可作为学习电力系统鲁棒优化的入门资料
最新发布
03-25
然而，可再生能源的引入给电力系统的稳定运行带来了前所未有的挑战，主要体现在以下几个方面： - **间歇性**：可再生能源发电量受天气条件的影响较大，具有很强的波动性和间歇性。 - **不确定性**：可再生能源发电的...
PPT周克敏-鲁棒控制最优化 共17章 bilibili鲁棒控制配套同步PPT
08-26
PPT周克敏-鲁棒控制最优化 bilibili鲁棒控制PPT
【Matlab代码+文档讲解】鲁棒优化由入门到精通一文讲透
05-09
【Matlab代码+文档讲解】鲁棒优化由入门到精通一文讲透 %% 下层用户侧目标函数:最低运行成本 %把购电与售电价格的差异给去掉，统一用购电价格。 Cost_EHbuy = EH_Price(1,1:24)*Ebuy+EH_Price(1,25:48)*Hbuy; %IES运行成本 Cost_YX = 0.7*sum(Eice) + 0.3*sum(Hgb); %IES维护成本 % Cost_WH = 0.015*sum(Epv+Ewt+abs(EessC)+abs(HessC)); Cost_WH = 0.015*sum(Epv+Ewt); Cost_IES = Cost_EHbuy + Cost_YX + Cost_WH ;
鲁棒优化线性规划（yalmip工具包）
03-14
%鲁棒线性优化 %利用yalmip求解鲁棒线性优化 % 1.简单线性规划示例 % % 让我们开始一个只有一个标量决策变量x和一个不确定的标量w的小问题。我们用一个不确定约束建立一个问题和一个简单的不确定模型。 % % sdpvar x w % F = [x+w <= 1]; % W = [-0.5 <= w <= 0.5, uncertain(w)]; % objective = -x; % 显然，最佳的x是0.5，因为如果x较大，就存在w使得不确定约束违例。 % 为了解决这个问题，我们称为 solvesdp。鲁棒对等式将自动生成和求解（一般来说，具有多面体不确定性的线性约束使用枚举方法求解，然而对于本例中简单方形不确定性，yalmip明确执行最大化，并导致一个更高效的最糟情况模型）。 %
鲁棒优化入门（一）——工具箱Xprog和RSOME的安装与使用
配电网和matlab的博客
07-07 6899
Xprog是由新加坡国立大学的Peng Xiong于2016年发布的一款matlab工具箱，可以用于求解确定性优化、随机优化、鲁棒优化和分布式鲁棒优化问题。 还有一个针对鲁棒随机优化设计的matlab工具箱——RSOME(Robust Stochastic Optimization Made Easy)，也可以用来解决一般的线性规划、随机规划、鲁棒优化和数据驱动的分布式鲁棒优化问题。...............
鲁棒优化（1）
Yudibrother的博客
04-26 1万+
前言 打算做一个鲁棒优化的系列，先从Bertsimas，2004的文章开始学习，最后的目标时做到分布鲁棒. 一、什么是鲁棒优化 简而言之，当你的模型存在不确定系数，而你需要免疫这些不确定系数，那你就需要考虑一下鲁棒优化. 二、基本的模型和假设 1.线性模型 我们假设，c矩阵不包含不确定变量，如果你的模型包含不确定变量，可以将目标函数化成约束的形式.扰动变量只在A阵出现.，约束的右端都是确定参数. 更具体的形式如上，公式（2）. 我们假设随机变量a，在一定范围内扰动，有上下界，且独立对称分布 a_hat
鲁棒优化入门(4)-两阶段鲁棒优化及行列生成算法(C&CG)超详细讲解(附matlab代码)
热门推荐
配电网和matlab的博客
05-20 1万+
鲁棒优化是应对数据不确定性的一种优化方法，但单阶段鲁棒优化过于保守。为了解决这一问题，引入了两阶段鲁棒优化(Two-stage Robust Optimization)以及更一般的多阶段鲁棒优化，其核心思想是将决策问题分为两个阶段。第一阶段是进行初步决策，第二阶段是根据第一阶段的决策结果制定更好的决策策略，应对数据不确定性的影响。这种方法可以降低保守性，提高鲁棒性。假设一阶段和二阶段决策问题都是线性规划，并且不确定性集合U是一个有限的离散集合或者多面体集。
鲁棒优化入门(6)—Matlab+Yalmip两阶段鲁棒优化通用编程指南(上)
配电网和matlab的博客
08-31 6152
Matlab+Yalmip两阶段鲁棒优化通用编程指南(上篇)，内容包括： ①.拿到一个复杂的两阶段鲁棒优化问题的分析步骤和方法。 ②.采用Yalmip工具箱中的uncertain函数和鲁棒优化模块求解两阶段鲁棒优化的子问题。 ③.Yalmip工具箱中的鲁棒优化模块和常规的求解思路有什么异同。 ④.使用KKT条件求解两阶段鲁棒优化的子问题，并使用C&CG算法进行迭代求解。
鲁棒优化（2）
Yudibrother的博客
04-27 4926
提示：文章写完后，目录可以自动生成，如何生成可参考右边的帮助文档 文章目录前言一、新的鲁棒优化模型二、使用步骤1.引入库2.读入数据总结 前言 上一篇文章以及讨论，最原始鲁棒优化，结果非常保守，但是它有一个非常重要的特点，它是线性模型，不需要很高深的知识就能看懂和求解.现在可以做鲁棒优化的工具箱很多，我们只需要理清思路，具体的实现可以交给工具箱。 言归正传，数学家针对，鲁棒优化的保守度，进行研究，分别提出了椭球不确定集合，和多面体不确定集。椭球不确定集合涉及到二阶锥，求解困难，而且不容易理解，目前在工程
分布式鲁棒优化基础知识学习 | Ref：《鲁棒优化入门》「运筹OR帷幄」
weixin_43905651的博客
10-15 9250
鲁棒：考虑最坏情况；分布：最坏情况的主体是环境参数的分布变量。从数学角度说，分布式鲁棒优化囊括随机规划和传统鲁棒优化两种形式。当分布式鲁棒优化下，环境变量的分布函数获知时，分布鲁棒优化退化为随机优化；仅知其不确定集时，退化为经典鲁棒优化。
利用gurobi写一个分布鲁棒的优化代码
03-30
以下是一个使用Gurobi实现分布鲁棒优化的示例代码： ``` from gurobipy import * import numpy as np # 定义数据 n = 1000 m = 10 np.random.seed(1) mu = np.random.rand(m) std = np.random.rand(m) A = np.random.rand(n, m) b = A.dot(mu) + 0.1 * np.random.randn(n) # 创建模型 model = Model("distributionally-robust-optimization") # 创建变量 x = model.addVars(m, lb=-GRB.INFINITY, name="x") # 定义目标函数 obj = quicksum(x[i] for i in range(m)) model.setObjective(obj, GRB.MINIMIZE) # 定义约束条件 for j in range(n): expr = quicksum(A[j, i] * x[i] for i in range(m)) model.addConstr(expr >= b[j] - quicksum(std[i] * x[i] for i in range(m))) # 解决模型 model.optimize() # 输出结果 if model.status == GRB.OPTIMAL: print("Optimal objective value:", model.objVal) print("Optimal solution:") for i in range(m): print("x[{}] = {}".format(i, x[i].x)) else: print("Optimization failed.") ``` 该代码解决了以下分布鲁棒优化问题： $$\begin{aligned} \min_{x} & \quad \sum_{i=1}^{m} x_i \\ \text{s.t.} & \quad \forall j \in \{1,\dots,n\}, \sum_{i=1}^{m} A_{ji} x_i \geq b_j - \sum_{i=1}^{m} \sigma_i x_i \end{aligned}$$ 其中，$x\in\mathbb{R}^m$为优化变量，$A\in\mathbb{R}^{n\times m}$和$b\in\mathbb{R}^n$为已知数据，$\sigma\in\mathbb{R}^m$为分布的不确定性参数，满足$\mathbb{E}[\sigma_i]=0$且$\|\sigma\|_2\leq 1$。该问题的目标是最小化$x$的$L_1$范数，使得所有约束条件都满足。 该代码使用了Gurobi中的快速求和函数`quicksum`和变量的下界设为负无穷`-GRB.INFINITY`。在定义约束条件时，使用了numpy中的矩阵乘法`dot`和元素乘法`*`，以及`quicksum`函数计算线性组合。最后，输出了最优解和目标函数的值。
评论 1
添加红包

请填写红包祝福语或标题

个

红包个数最小为10个

元

红包金额最低5元

当前余额3.43元 前往充值 >
需支付：10.00元
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付元
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明：

1.余额是钱包充值的虚拟货币，按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载，可以购买VIP、付费专栏及课程。

余额充值