精华内容
下载资源
问答
  • 将区块链引入虚拟电厂(VPP)的调度运行机制中,针对新能源参与的电力系统模型,提出适用于VPP的实用拜占庭容错算法共识机制以实现区块链下半中心化的两阶段鲁棒优化调度模型,保留了VPP控制中心的导向作用。...
  • 包含储能,微型燃气机,功率平衡约束,配电网交互等约束,具体约束可参考《微电网两阶段鲁棒优化经济调度方法_刘一欣》,有一定相似性,但不包含全部。 程序保证稳定收敛,并且注释清晰,模型,推导过程均有文件和...
  • 两阶段鲁棒优化及列和约束生成算法

    千次阅读 多人点赞 2020-10-18 21:06:47
    两阶段鲁棒优化及列和约束生成算法1. 前言2. 阶段RO3. Benders对偶割平面法 1. 前言 有同学私信我两阶段鲁棒优化的问题,自己之前主要研究单阶段的鲁棒优化,对于阶段优化不太懂。查了点资料,通过翻译和自己的...

    1. 前言

    有同学私信我两阶段鲁棒优化的问题,自己之前主要研究单阶段的鲁棒优化,对于两阶段优化不太懂。查了点资料,通过翻译和自己的理解,写下这篇博文,抛砖引玉,以供大家共同学习交流。

    2. 两阶段RO

    本文主要翻译自2013年Bo Zeng 博士的高被引论文《Solving two-stage robust optimization problems using a column-and-constraint generation method》[1]。

    由于传统单阶段的鲁棒优化对于所有不确定性是完全免疫的,求出的结果一般过于保守、比较悲观。两阶段RO也称鲁棒可调优化或者自适应优化,它的引入就是为了应对传统RO存在的上述问题。两阶段RO是指,在作出第一阶段的决策和不确定性部分展现出来之后,再进行第二阶段的决策。 由于改进的建模能力,两阶段RO,广泛应用于网络/运输问题、投资组合优化和电力系统调度等问题。

    文献[1]假设一阶段和二阶段决策问题都是线性规划,并且不确定性集合 U \bf U U是离散有限的点集或者多面体。使用 y \bf y y表示第一阶段决策变量, x \bf x x表示第二阶段决策变量, u ∈ U u \in \bf U uU表示不确定矢量。在此假设下的两阶段RO的一般形式为:
    min ⁡ y c T y + max ⁡ u ∈ U min ⁡ x ∈ F ( y , u ) b T x (1) \min_{\bf y } \bf c^{\rm T}y+\rm \max_{u \in \bf U} \min_{\bf x \in F(y,\rm u)} \bf b^{\rm T}x \tag{1} ymincTy+uUmaxxF(y,u)minbTx(1) s . t . A T y ≥ d , y ∈ S y \rm s.t. \quad \bf A^{\rm T}y \geq d,y \in S_{y} s.t.ATyd,ySy

    其中, F ( y , u ) = { x : G x ≥ h − E y − M u , x ∈ S x } \bf F(y,\rm u)=\{\bf x: Gx \geq h-Ey-M \rm u,\bf x \in S_{x}\} F(y,u)={x:GxhEyMu,xSx} S y ⊆ R + n \bf S_{y}\subseteq \Bbb R_{+}^{n} SyR+n S x ⊆ R + m \bf S_{x}\subseteq \Bbb R_{+}^{m} SxR+m,向量 c , b , d , h \bf c,b,d,h c,b,d,h和矩阵 A , G , E , M \bf A,G,E,M A,G,E,M都是确定性的数值,不确定性体现在向量 u u u上。注意到第二阶段优化的约束条件 F ( y , u ) \bf F(y,\rm u) F(y,u)是关于不确定性 u u u的线性函数。

    两阶段RO有优势,当然有缺点。即便是最简单的两阶段RO,也可以是NP-hard问题,计算复杂度很高。为了减轻计算负担,一般有两种方法。第一种是使用近似算法,该种方法假设第二阶段的决策是关于不确定性的简单函数,例如仿射函数。第二种方法试图根据Benders分解得出精确解,即它们使用第二阶段决策问题的对偶解,逐步构造第一阶段决策的值函数(value function)。一般称为Benders对偶割平面算法。

    3. Benders对偶割平面法

    由于第二阶段的决策是关于 x \bf x x的线性规划问题,可以作以下假设(relatively complete recourse assumption):该线性规划对于任意给定的 y \bf y y u u u是可行的,也即是有解 (该假设不是很理解)。假设第二阶段的线性规划的对偶变量为 π \pi π,则将其转化为对偶问题为:
    S P 1 : O ( y ) = max ⁡ u , π { ( h − E y − M u ) T π } (2) \rm {SP_{1}}: \mathcal O(\bf y) =\rm \max_{u,\pi} \{ \bf (h-Ey-M \rm u)^{T}\pi\} \tag{2} SP1:O(y)=u,πmax{(hEyMu)Tπ}(2) s . t .    G T π ≤ b , u ∈ U , π ≥ 0 s.t. \; \bf G^{\rm T}\pi \leq b, \rm u \in \bf U, \pi \geq \bf 0 s.t.GTπb,uU,π0

    在对偶问题中,目标函数从原始的最小化转换为关于对偶变量 π \pi π的最大化,同时与(1)式中的最大化 u u u合并,得到上述子问题(2)。此时不确定性向量转化为对偶问题(2)的决策变量,需要注意到子问题(2)是存在 u T π u^{T}\pi uTπ的双线性项。此时问题(2)被称为关于 u , π u,\pi u,π的双线性规划。对于双线性规划的求解,放在后面再说。

    假设对于给定的 y k ∗ \bf y_{\mathit k}^{*} yk,子问题(2) 的最优解为 ( u k ∗ , π k ∗ ) (u_{ k}^{*},\pi_{ k}^{*}) (uk,πk),根据对偶定理,则可以构建以下割平面:
    η ≥ ( h − E y − M u k ∗ ) T π k ∗ (3) \eta \geq \bf (h-Ey-M \rm u_{\mathit k}^{*})^{T}\pi_{ \mathit k}^{*} \tag{3} η(hEyMuk)Tπk(3)
    其中, η = max ⁡ u ∈ U min ⁡ x ∈ F ( y , u ) b T x \eta = \max_{u \in \bf U} \min_{\bf x \in F(y,\rm u)} \bf b^{\rm T}x η=maxuUminxF(y,u)bTx为一维标量。注意该割平面是关于第一阶段决策变量 y \bf y y的约束。

    将割平面约束添加第一阶段优化中,可以得到:
    M P 1 : min ⁡ y c T y + η (4) \rm MP_{1}: \min_{\bf y } \bf c^{\rm T}y+ \eta \tag{4} MP1:ymincTy+η(4) s . t . A T y ≥ d , y ∈ S y \rm s.t. \quad \bf A^{\rm T}y \geq d,y \in S_{y} s.t.ATyd,ySy η ≥ ( h − E y − M u l ∗ ) T π l ∗ , ∀    l ≤ k \eta \geq \bf (h-Ey-M \rm u_{\mathit l}^*)^{T}\pi_{ \mathit l}^* ,\forall \; \mathit {l \leq k} η(hEyMul)Tπl,lk y ∈ S y , η ∈ R \bf y \in S_{y}, \eta \in \Bbb R ySy,ηR

    其中 ( u l ∗ , π l ∗ ) (u_{ l}^{*},\pi_{ l}^{*}) (ul,πl)已知。求解主问题(4)可以获得其最优解 ( y k + 1 ∗ , η k + 1 ∗ ) (\bf y_{\mathit k+1}^{*},\eta_{\mathit k+1}^{*}) (yk+1,ηk+1)

    此时问题(4)的目标函数值 c T y k + 1 ∗ + η k + 1 ∗ \bf c^{\rm T}y_{\mathit k+1}^{*}+\eta_{\mathit k+1}^{*} cTyk+1+ηk+1是原始两阶段问题(1)的下界(lower bound, LB)。如何理解?原始的两阶段RO是在worse-case情况下的目标值,也即是对于所有不确定性都具有包容性,可以理解成其最优决策在所有的不确定性下都成立。通过把不确定性转化为割平面,主问题(4)中添加的割平面只是在部分不确定性 u l u_{l} ul下的优化,也即是问题(4)是关于问题(1)的松弛问题,由于是求最小化,问题(4)的最优目标值一定是原问题(1)的一个下界。

    注意到 c T y k ∗ + O ( y k ∗ ) \bf c^{\rm T}y_{\mathit k}^{*}+\mathcal O(\bf y_{\mathit k}^{*}) cTyk+O(yk)是原始问题(1)的上界。如何理解? y k ∗ \bf y_{\mathit k}^{*} yk是第一阶段的一个可行解, O ( y k ∗ ) \mathcal O(\bf y_{\mathit k}^{*}) O(yk)是第二阶段优化在 y k ∗ \bf y_{\mathit k}^{*} yk下的最优目标值,也即是说存在二阶段的最优解 x k ∗ \bf x_{\mathit k}^* xk。此时的 y k ∗ , x k ∗ \bf y_{\mathit k}^{*},\bf x_{\mathit k}^{*} yk,xk只是问题(1)的可行解,而不一定是最优解。由于是求最小化问题, c T y k ∗ + O ( y k ∗ ) \bf c^{\rm T}y_{\mathit k}^{*}+\mathcal O(\bf y_{\mathit k}^{*}) cTyk+O(yk)是原始问题(1)的上界。

    y k + 1 ∗ \bf y_{\mathit k+1}^{*} yk+1带入到子问题(2)中,求解问题(2)得到最优解 ( u k + 1 ∗ , π k + 1 ∗ ) (u_{ k+1}^{*},\pi_{ k+1}^{*}) (uk+1,πk+1) 和最优值 O ( y k + 1 ∗ ) \mathcal O(\bf y_{\mathit k+1}^{*}) O(yk+1),进而可以构造新的割平面。因此,迭代的引入割平面(3),计算主问题(1),上届和下界逐渐收敛到问题(1)的最优解。

    计算复杂度[1]
    命题1:如果不确定集合 U \bf U U是多面体,用 p p p表示极点的数量,如果 U \bf U U是离散点集,用 p p p表示集合的势(集合中元素的数量)。用 q q q表示多面体 { π : G T π ≤ b , π ≥ 0 } \{\pi: \bf G^{\rm T}\pi \leq b, \pi \geq 0 \} {π:GTπb,π0} 的极点数量。Benders对偶算法需要经过 o ( p q ) o(pq) o(pq)次迭代才能求出问题(1)的最优解。

    Benders 对偶割平面法主要的问题是求解问题(2)的双线性规划问题

    4. 列和约束生成(C&CG)算法

    列和约束生成(column-and-constraint generation) 算法是基于以下事实:
    假设 x l \bf x^{\mathit l} xl是不确定性 u = u l u=u_{l} u=ul(一个实例)下 η = max ⁡ u ∈ U min ⁡ x ∈ F ( y , u ) b T x \eta = \max_{u \in \bf U} \min_{\bf x \in F(y,\rm u)} \bf b^{\rm T}x η=maxuUminxF(y,u)bTx的最优解,则一定存在如下约束成立:
    η ≥ b T x l \eta \geq \bf b^{\rm T}x^{\mathit l} ηbTxl G x l ≥ h − E y − M u l \bf Gx^{\mathit l} \geq h-Ey-M \rm u_{\mathit l} GxlhEyMul
    可以根据以上约束构造割平面,形成C&CG算法。

    Step 1 :设置 L B = − ∞ \rm LB=-\infty LB= U B = + ∞ \rm UB=+\infty UB=+,索引 k = 0 k=0 k=0,集合 O = ∅ \bf O=\emptyset O=
    Step 2: 求解如下主问题:
    M P 2 : min ⁡ y c T y + η (5) \rm MP_{2}: \min_{\bf y } \bf c^{\rm T}y+ \eta \tag{5} MP2:ymincTy+η(5) s . t . A T y ≥ d , y ∈ S y \rm s.t. \quad \bf A^{\rm T}y \geq d,y \in S_{y} s.t.ATyd,ySy η ≥ b T x l , ∀ l ∈ O \eta \geq \bf b^{\rm T}x^{\mathit l}, \forall \mathit {l \in \bf O} ηbTxl,lO G x l ≥ h − E y − M u l ∗ ∀ l ≤ k \bf Gx^{\mathit l} \geq h-Ey-M \rm u_{\mathit l}^* \forall \mathit {l \leq k} GxlhEyMullk y ∈ S y , η ∈ R , x l ∈ S x    ∀ l ≤ k \bf y \in S_{y}, \eta \in \Bbb R, x^{\mathit l} \in S_{x} \; \forall \mathit {l \leq k} ySy,ηR,xlSxlk
    求出最优解 ( y k + 1 ∗ , η k + 1 ∗ , x 1 ∗ , ⋯   , x k ∗ ) (\bf y_{\mathit k+1}^*,\eta_{\mathit k+1}^*,x^{\mathit 1*},\cdots,x^{\mathit k*}) (yk+1,ηk+1,x1,,xk),更新下界 L B = c T y k + 1 ∗ + η k + 1 ∗ \rm LB=\bf c^{\rm T}y_{\mathit k+1}^{*}+\eta_{\mathit k+1}^{*} LB=cTyk+1+ηk+1

    Step 3 :代入 y = y k + 1 ∗ \bf y=y_{\mathit k+1}^* y=yk+1,求解如下子问题:
    S P 2 : O ( y ) = max ⁡ u ∈ U min ⁡ x ∈ F ( y , u ) { b T x : G x ≥ h − E y − M u , x ∈ S x } (6) \rm SP2: \mathcal O(\bf y)= \rm \max_{u \in \bf U} \min_{\bf x \in F(y,\rm u)} \{ \bf b^{\rm T}x :\bf Gx \geq h-Ey-M \rm u, \bf x \in S_{x} \} \tag{6} SP2:O(y)=uUmaxxF(y,u)min{bTx:GxhEyMu,xSx}(6)
    更新上界 U B = min ⁡ { U B , c T y k + 1 ∗ + O ( y k + 1 ∗ ) } \rm UB=\min\{UB,\bf c^{\rm T}y_{\mathit k+1}^{*}+\mathcal O(\bf y_{\mathit k+1}^{*})\} UB=min{UB,cTyk+1+O(yk+1)}

    Step 4 : 如果 U B − L B ≤ ϵ \rm UB-LB\leq \epsilon UBLBϵ, 返回 y k + 1 ∗ \bf y_{\mathit k+1}^* yk+1,程序终止。否则:
    (a) 如果 O ( y k + 1 ∗ ) < + ∞ \mathcal O(\bf y_{\mathit k+1}^{*}) < +\infty O(yk+1)<+,添加变量 x k + 1 \bf x^\mathit {k+1} xk+1,添加如下约束:
    η ≥ b T x k + 1 (7) \eta \geq \bf b^{\rm T}x^{\mathit {k+1}} \tag{7} ηbTxk+1(7) G x k + 1 ≥ h − E y − M u k + 1 ∗ (8) \bf Gx^{\mathit {k+1}} \geq h-Ey-M \rm u_{\mathit {k+1}}^* \tag{8} Gxk+1hEyMuk+1(8)
    到问题(5)中。其中 u k + 1 ∗ u_{\mathit {k+1}}^* uk+1是问题(6)在 y = y k + 1 ∗ \bf y=y_{\mathit k+1}^* y=yk+1下的最优场景(不确定性),对问题(6)可以利用数据库(Call the oracle to solve subproblem SP2)进行枚举求得。
    然后,更新 k = k + 1 k=k+1 k=k+1,集合 O = O ⋃ k + 1 \bf O= O \bigcup {\mathit {k+1}} O=Ok+1,跳转至步骤Step 2

    (b) 如果 O ( y k + 1 ∗ ) = + ∞ \mathcal O(\bf y_{\mathit k+1}^{*}) =+\infty O(yk+1)=+ (对于某些 u ∗ ∈ U u^*\in \bf U uU,如果第二阶段决策 O ( y k + 1 ∗ ) \mathcal O(\bf y_{\mathit k+1}^{*}) O(yk+1)不可行(infeasible),则把 O ( y k + 1 ∗ ) \mathcal O(\bf y_{\mathit k+1}^{*}) O(yk+1)记为 + ∞ +\infty +),添加变量 x k + 1 \bf x^\mathit {k+1} xk+1,添加如下约束:
    G x k + 1 ≥ h − E y − M u k + 1 ∗ (9) \bf Gx^{\mathit {k+1}} \geq h-Ey-M \rm u_{\mathit {k+1}}^* \tag{9} Gxk+1hEyMuk+1(9) 到问题(5)中。其中 u k + 1 ∗ u_{\mathit {k+1}}^* uk+1是问题(6)在 y = y k + 1 ∗ \bf y=y_{\mathit k+1}^* y=yk+1下不可行的不确定性 u u u的值。
    然后,更新 k = k + 1 k=k+1 k=k+1,跳转至步骤Step 2

    对于步骤Step 4(a)的约束(7)和(8)被称为最优割(optimality cuts),而Step 4(b)的约束(9)被称为可性割(feasibility cuts)。对于可行割不是很明白,虽然在某些 y = y k + 1 ∗ \bf y=y_{\mathit k+1}^* y=yk+1下不可行,但是仍然是其中一个不确定性 u k + 1 ∗ u_{\mathit {k+1}}^* uk+1对二阶段决策的一个反映,因为二阶段是包含所有的不确定性。 不知以上理解是否正确。 由于约束(7)对于不可行的场景仍然是成立的,此时的步骤4(a)和4(b)可以合并,因此可以用统一的方法和程序去处理最优性和可行性。

    由于生成的割平面是由第二阶段决策变量(wait-and-see decision variables/recourse decision variables)以二阶段决策约束的形式定义的,整个过程其实就是列和约束生成,其中列生产是指在主问题中添加第二阶段决策变量,约束生成是指添加割平面约束。

    C&CG算法的复杂性
    命题2:如果不确定集合 U \bf U U是多面体,用 p p p表示极点的数量,如果 U \bf U U是离散点集,用 p p p表示集合的势(集合中元素的数量)。C&CG算法需要经过 o ( p ) o(p) o(p)次迭代才能求出问题(1)的最优解。

    C&CG算法与Benders对偶的区别:
    – (1) 主问题中决策的决策变量。C&CG算法每次迭代都添加决策变量 x l \bf x^\mathit l xl,而Benders对偶每次迭代决策变量不变。

    – (2)计算复杂度。由以上两个命题可以看出,在( the relatively complete recourse assumption)假设下,C&CG算法具有更低的复杂度。

    – (3)求解问题的能力。Benders分解要求第二阶段的决策问题为行性规划问题,C&CG算法对于第二阶段优化问题的变量的类型不敏感,可以是混合整数规划问题。

    5 存在的难点

    以上两种算法存在的难点在于:Benders对偶算法求解子问题(2)的双线性规划;C&CG算法求解子问题(6)。

    针对相对简单的多面体不确定性集,一般使用外部近似算法(outer approximation algorithm)和混合整数线性重构(mixed integer linear reformulation)两种方法求解。 第一种是求解SP1的启发式算法。 第二种是使用不确定性集的特殊结构,将双线性规划SP1转换为等效的混合整数线性规划。 文献[1] 使用经典的Karush–Kuhn–Tucker(KKT)条件来处理一般的多面体不确定性集,只要( the relatively complete recourse assumption)假设成立即可。文献[1]对于问题(6)SP2,通过引入对偶变量、利用KKT、结合大M方法,将鲁棒优化问题(6)转化为确定性的整数线性规划问题,再利用现成的求解器进行求解。具体转化就不翻译了,自己没看懂。

    文献[2]中表述:根据文献[3]的结论,问题(2)的双线性规划,变量 u u u的最优解为不确定集合 U U U的极点。对于多面体的不确定集合,其极点是有限的,需要找出所有的极点,然后通过枚举就可以求出最优的 u u u。同样根据文献[4],在问题(2)达到最优时,变量 u u u 和对偶变量 π \pi π取值总是其各自可行集的极点,因此也需要求出 π \pi π的极点。如果变量 u u u π \pi π的维数比较大的话,枚举的复杂度仍然比较高。

    对于问题(6)本实质上是单阶段的RO问题,可以使用求解一般RO的求解器进行求解。但是它们与传统的RO仍有区别,不论是问题(2)还是(6)都需要求出具体的不确定性 u u u的值。可以认为,传统的线性RO其最大的不确定性,在不确定性集合的极点处取得,因此还是需要计算不确定集合 U U U的极点。先使用单阶段RO求解出决策变量 x l \bf x^{\mathit l} xl,然后在枚举所有的 U U U的极点,找到满足临界约束条件的点,应该就是具体的 u u u。不知理解是否正确,后面有时间再写代码吧。

    参考文献
    [1] Zeng B , Zhao L . Solving Two-stage Robust Optimization Problems by A Constraint-and-Column Generation Method[J]. Operations Research Letters, 2013, 41(5):457-461.
    [2]刘一欣, 郭力, 王成山. 微电网两阶段鲁棒优化经济调度方法[J]. 中国电机工程学报, 2018, 038(014):4013-4022.
    [3] D. Bertsimas, E. Litvinov, X.A. Sun, Jinye Zhao, Tongxin Zheng, Adaptive robust optimization for the security constrained unit commitment problem, IEEE Transactions on Power Systems 28 (1) (2013) 52–63.
    [4] Bo Zeng. Solving Two-stage Robust Optimization Problems by A Constraint-and-Column Generation Method. 2011. http://www.optimization-online.org/DB_FILE/2011/06/3065.pdf

    展开全文
  • 两阶段鲁棒优化CCG列于约束生成和Benders代码,可扩展改编,复现自原论文。文件中附源代码以及论文。使用matlab-yalmip编
  • 两阶段鲁棒优化程序 采用微网为模型,主要将安装成本、运营成本以及综合效益三个方面纳入考虑范围,建立两阶段鲁棒优化模型,采用的是CCG方法,本程序为matlab编制,另外本程序考虑发电设备、风光储的容量配置和...

    两阶段鲁棒优化程序

    采用微网为模型,主要将安装成本、运营成本以及综合效益三个方面纳入考虑范围,建立两阶段鲁棒优化模型,采用的是CCG方法,本程序为matlab编制,另外本程序考虑发电设备、风光储的容量配置和出力情况,考虑风光负荷的不确定性

    展开全文
  • #资源达人分享计划#
  • 行业分类-物理装置-一种基于两阶段鲁棒优化模型的手术调度优化方法及系统.zip
  • 行业分类-物理装置-一种微电网能量两阶段鲁棒优化方法及系统.zip
  • 采用yalmip编的两阶段鲁棒优化,目标函数主要考虑了投资成本(第一阶段)和运行成本(第二阶段)部分,其中,投资成本主要为储能的等年值投资成本,运行成本则包括配电网交互成本(购售电成本)、各单元运维成本...
  • 在这个案例中,我们用RSOME建模求解一个两阶段批量生产问题。这里我们演示如何通过定义不同形式的含糊集合来实现相应的分布鲁棒优化模型。


    RSOME: Robust Stochastic Optimization Made Easy


    Adaptive Distributionally Robust Lot-Sizing

    In this section, we are using RSOME to replicate numerical case studies presented in Bertsimas et al. (2021). A capacitated network with n n n locations is considered, where each location i i i has an unknown demand d i d_i di, and the demand can be satisfied by the existing local stock x i x_i xi or by transporting an amount y i j y_{ij} yij of units from another location j j j, which is determined after the demand is realized. Such a lot-sizing problem can be written as the following two-stage formulation,

    min ⁡   ∑ i = 1 N c i x i + E [ Q ( x , d ~ ) ] , s.t.  0 ≤ x i ≤ K , i = 1 , 2 , . . . , n \begin{aligned} \min~& \sum\limits_{i=1}^Nc_ix_i + \mathbb{E}\left[Q(\pmb{x}, \tilde{\pmb{d}})\right], \\ \text{s.t.}~ &0 \leq x_i \leq K, & i = 1, 2, ..., n \end{aligned} min s.t. i=1Ncixi+E[Q(xxx,ddd~)],0xiK,i=1,2,...,n

    where the recourse problem Q ( x , d ) Q(\pmb{x}, \pmb{d}) Q(xxx,ddd) is written as

    Q ( x , z ) = min ⁡   ∑ i = 1 n ∑ j = 1 n c i j y i j s.t.  x i − ∑ j = 1 N y j i + ∑ j = 1 N y i j ≥ d i i = 1 , 2 , . . . , n 0 ≤ y ≤ b . \begin{aligned} Q(\pmb{x}, \pmb{z}) = \min~& \sum\limits_{i=1}^n\sum\limits_{j=1}^nc_{ij}y_{ij} &\\ \text{s.t.}~&x_i - \sum\limits_{j=1}^Ny_{ji} + \sum\limits_{j=1}^Ny_{ij} \geq d_i & i = 1, 2, ..., n \\ &\pmb{0} \leq \pmb{y} \leq \pmb{b}. \end{aligned} Q(xxx,zzz)=min s.t. i=1nj=1ncijyijxij=1Nyji+j=1Nyijdi000yyybbb.i=1,2,...,n

    Parameters of the case studies are given below:

    • The number of locations: n = 10 n=10 n=10;
    • Cost of local stock: a i = 1 a_i=1 ai=1, i = 1 , 2 , . . . , n i=1, 2, ..., n i=1,2,...,n;
    • Cost of shipping units: c i j c_{ij} cij is the Euclidean distance between location i i i and j j j.
    • Stock capacity: K = 20 K=20 K=20;
    • Transportation capacity: b i j = K / ( n − 1 ) u i j b_{ij}=K/(n-1)u_{ij} bij=K/(n1)uij, where each u i j u_{ij} uij is a random variable generated from a standard uniform distribution;
    • Sample size of the training set: N = 20 N=20 N=20;
    • Robustness parameter: ϵ = 10 N − 1 / 10 \epsilon=10N^{-1/10} ϵ=10N1/10.
    n = 10
    a = 1
    K = 20
    N = 20
    epsilon = 10 * N**(-1/10)
    
    b = K / (n-1) / rd.rand(n, n)
    
    xy = 10*rd.rand(2, n)
    c = ((xy[[0]] - xy[[0]].T) ** 2
         + (xy[[1]] - xy[[1]].T) ** 2) ** 0.5
    

    Here we generate a training dataset ds assuming the demand at each location follows a uniform distribution. The average sample demand at each location is shown by the following figure. The lot-sizing problem is then solved by approaches described in Bertsimas et al. (2021): sample average approximation (SAA), single-policy affine approximation (SP affine), multi-policy affine approximation (MP affine), and the Wasserstein scenario-wise adaptation (Wass SW).

    在这里插入图片描述

    SAA

    A direct means of implementing the SAA model is to define the recourse decision y \pmb{y} yyy as a three-dimensional array, where the first dimension represents the sample size, and the remaining two dimensions indicate the from and to locations.

    from rsome import dro
    import rsome.grb_solver as grb
    import rsome as rso
    
    model = dro.Model()             # define a model
    x = model.dvar(n)               # here-and-now location decisions
    y = model.dvar((N, n, n))       # wait-and-see transportation decisions
    
    # Define model objective and constraints
    model.min((a*x).sum() + (1/N)*(c*y.sum(axis=0)).sum())
    model.st(ds <= y.sum(axis=1) - y.sum(axis=2) + x)
    model.st(y >= 0, y <= b)
    model.st(x >= 0, x <= K)
    
    model.solve(grb)
    
    Being solved by Gurobi...
    Solution status: 2
    Running time: 0.0119s
    

    Alternatively, the SAA method could be implemented by defining a scenario-wise ambiguity set F \mathbb{F} F with singleton supports for each scenario,

    F = { P ∈ P ( R n × [ N ] ) ∣   ( d ~ , s ) ∼ P P [ d ~ = d ^ s ∣ s ~ = 1 ] = 1 , ∀ s ∈ [ N ] P [ s ~ = s ] = 1 , ∀ s ∈ [ N ] for some  p ∈ P } \begin{aligned} \mathbb{F} = \left\{ \mathbb{P} \in \mathcal{P}(\mathbb{R}^n\times [N])\left| \begin{array}{ll} ~(\tilde{\pmb{d}}, s) \sim \mathbb{P} & \\ \mathbb{P}[\tilde{\pmb{d}} = \hat{\pmb{d}}_s|\tilde{s}=1]=1, &\forall s\in [N] \\ \mathbb{P}[\tilde{s}=s] = 1, &\forall s \in [N] \\ \text{for some }\pmb{p} \in \mathcal{P} & \end{array} \right. \right\} \end{aligned} F=PP(Rn×[N]) (ddd~,s)PP[ddd~=ddd^ss~=1]=1,P[s~=s]=1,for some pppPs[N]s[N]

    with d ^ s \hat{\pmb{d}}_s ddd^s being each sample record of demands. Model with such an ambiguity set can be implemented by the code below.

    model = dro.Model(N)    # define a model with N scenarios
    x = model.dvar(n)       # here-and-now location decisions
    y = model.dvar((n, n))  # wait-and-see transportation decisions
    d = model.rvar(n)       # define random variables as the demand
    
    for s in range(N):
        y.adapt(s)          # the decision rule y adapts to all scenarios
    
    # define the ambiguity set
    fset = model.ambiguity()
    for s in range(N):
        fset[s].suppset(d == ds[s])
    pr = model.p
    fset.probset(pr == 1/N)
    
    # define model objective and constraints
    model.minsup((a*x).sum() + E((c*y).sum()), fset)
    model.st(d <= y.sum(axis=0) - y.sum(axis=1) + x)
    model.st(y >= 0, y <= b)
    model.st(x >= 0, x <= K)
    
    model.solve(grb)
    
    Being solved by Gurobi...
    Solution status: 2
    Running time: 0.0193s
    

    The optimal stock decision provided by the SAA model is shown by the following figure.

    在这里插入图片描述

    SP affine

    The second data-driven approach is referred to as sample robust optimization, where an uncertainty set around each data sample is considered. Here the recourse decisions are approximated by a single-policy linear decision rule, thus the name SP affine. The ambiguity set of the sample robust model is presented below.

    F = { P ∈ P ( R n × [ N ] ) ∣   ( d ~ , s ) ∼ P P [ ∥ d ~ − d ^ s ∥ ≤ ϵ ∣ s ~ = 1 ] = 1 , ∀ s ∈ [ N ] P [ s ~ = s ] = 1 , ∀ s ∈ [ N ] for some  p ∈ P } \begin{aligned} \mathbb{F} = \left\{ \mathbb{P} \in \mathcal{P}(\mathbb{R}^n\times [N])\left| \begin{array}{ll} ~(\tilde{\pmb{d}}, s) \sim \mathbb{P} & \\ \mathbb{P}[\|\tilde{\pmb{d}} - \hat{\pmb{d}}_s\|\leq \epsilon|\tilde{s}=1]=1, &\forall s\in [N] \\ \mathbb{P}[\tilde{s}=s] = 1, &\forall s \in [N] \\ \text{for some }\pmb{p} \in \mathcal{P} & \end{array} \right. \right\} \end{aligned} F=PP(Rn×[N]) (ddd~,s)PP[ddd~ddd^sϵs~=1]=1,P[s~=s]=1,for some pppPs[N]s[N]

    The code for implementing such an ambiguity set is given below.

    model = dro.Model(N)    # define a model with N scenarios
    x = model.dvar(n)       # here-and-now location decisions
    y = model.dvar((n, n))  # wait-and-see transportation decisions
    d = model.rvar(n)       # define random variables as the demand
    
    y.adapt(d)              # the decision rule y affinely depends on d
    
    # define the ambiguity set
    fset = model.ambiguity()
    for s in range(N):
        fset[s].suppset(d >= 0, d <= K, rso.norm(d - ds[s]) <= epsilon)
    pr = model.p
    fset.probset(pr == 1/N)
    
    # define model objective and constraints
    model.minsup((a*x).sum() + E((c*y).sum()), fset)
    model.st(d <= y.sum(axis=0) - y.sum(axis=1) + x)
    model.st(y >= 0, y <= b)
    model.st(x >= 0, x <= K)
    
    model.solve(grb)
    
    Being solved by Gurobi...
    Solution status: 2
    Running time: 38.8390s
    

    The optimal lot-sizing decisions of the SP affine sample robust optimization method would be more conservative compared with the SAA method, as shown by the following figure.

    在这里插入图片描述

    MP affine

    The sample robust optimization model above could be extended to the MP affine version, where the recourse decision y \pmb{y} yyy not only affinely depends on d \pmb{d} ddd, but also adapts to various scenarios. The code together with the optimal allocation of stocks are presented below. It can be seen that the solution is less conservative compared with the SP affine case, because the recourse decision has a higher level of flexibility in approximating the actual recourse decisions.

    model = dro.Model(N)    # define a model with N scenarios
    x = model.dvar(n)       # here-and-now location decisions
    y = model.dvar((n, n))  # wait-and-see transportation decisions
    d = model.rvar(n)       # define random variables as the demand
    
    y.adapt(d)              # the decision rule y affinely depends on d
    for s in range(N):
        y.adapt(s)          # the decision rule y also adapts to each scenario s
    
    # define the ambiguity set
    fset = model.ambiguity()
    for s in range(N):
        fset[s].suppset(d >= 0, d <= K, rso.norm(d - ds[s]) <= epsilon)
    pr = model.p
    fset.probset(pr == 1/N)
    
    # define model objective and constraints
    model.minsup((a*x).sum() + E((c*y).sum()), fset)
    model.st(d <= y.sum(axis=0) - y.sum(axis=1) + x)
    model.st(y >= 0, y <= b)
    model.st(x >= 0, x <= K)
    
    model.solve(grb)
    
    Being solved by Gurobi...
    Solution status: 2
    Running time: 19.4076s
    

    在这里插入图片描述

    Wass SW

    As pointed out by the paper Bertsimas et al. (2021), the distributionally robust optimization model considering the Wasserstein ambiguity set proposed in Chen et al. (2020) only yield one feasible first-stage decision x = ( K , K , . . . , K ) \pmb{x}=(K, K, ..., K) xxx=(K,K,...,K) because this is the only solution that guarantees the feasibility of the second-stage problem over the entire support of d \pmb{d} ddd. In order to make the numerical implementation more meaningful, we introduced an emergency order w i w_i wi to fill up the shortage at each location i i i. The cost of the emergency order is assumed to be five times of stock cost a i a_i ai. The code for implementing this distributionally robust optimization model is given below.

    model = dro.Model(N)    # define a model
    x = model.dvar(n)       # define location decisions
    y = model.dvar((n, n))  # define decision rule as the shifted quantities
    w = model.dvar(n)       # define decision rule as the emergency cost
    d = model.rvar(n)       # define random variables as the demand
    u = model.rvar()        # define an auxiliary random variable
    
    y.adapt(d)              
    y.adapt(u)
    w.adapt(d)
    w.adapt(u)
    for s in range(N):
        y.adapt(s)
        w.adapt(s)
    
    # define the uncertainty set
    fset = model.ambiguity()
    for s in range(N):
        fset[s].suppset(d >= 0, d <= K, rso.norm(d - ds[s], 1) <= u)
    fset.exptset(E(u) <= epsilon)
    pr = model.p
    fset.probset(pr == 1/N)
    
    # define model objective and constraints
    model.minsup((a*x).sum() + E((c*y).sum() + (5*a*w).sum()), fset)
    model.st(d <= y.sum(axis=0) - y.sum(axis=1) + x + w)
    model.st(y >= 0, y <= b)
    model.st(x >= 0, x <= K)
    model.st(w >= 0)
    
    model.solve(grb)
    
    Being solved by Gurobi...
    Solution status: 2
    Running time: 19.1652s
    

    The optimal lot-sizing decision is illustrated by the following figure.

    在这里插入图片描述

    Reference

    Bertsimas, Dimitris, Shimrit Shtern, and Bradley Sturt. 2021. Two-stage sample robust optimization. Operations Research.

    Chen, Zhi, Melvyn Sim, Peng Xiong. 2020. Robust stochastic optimization made easy with RSOME. Management Science 66(8) 3329–3339.

    展开全文
  • 为了全面和准确地考虑风电出力的不确定性和消纳能力,并兼顾系统运行的经济性和可靠性,通过在风电不确定区间可优化鲁棒区间经济调度模型中引入常规机组和储能系统运行状态的离散决策变量,建立风储联合运行的双层...
  • 鲁棒优化研究综述

    2016-10-17 22:11:24
    鲁棒优化研究综述
  • 考虑下游手术排程的两阶段 鲁棒优化
  • distributionally_robust_optimization 论文中实现的方法: 约束随机系统的分布鲁棒控制 使用Wasserstein指标的数据驱动的分布式鲁棒优化:性能保证和易于重构
  • 参考论文:Rodrigues, Filipe, and Agostinho Agra. "An exact robust approach for the integrated berth allocation and quay crane scheduling problem under uncertain arrival times."European Journal of ...

    参考论文:Rodrigues, Filipe, and Agostinho Agra. "An exact robust approach for the integrated berth allocation and quay crane scheduling problem under uncertain arrival times." European Journal of Operational Research (2021).

    展开全文
  • #资源达人分享计划#
  • 使用公共数据集,结果表明:(i)基于场景的鲁棒优化方法优于基于确定性,随机和鲁棒优化的基准; (ii)数据驱动的优化优于忽略协变量信息的基准; (iii)我们的剖切面算法通过在较短的计算时间内返回更好的解决...
  • #资源达人分享计划#
  • 鲁棒优化(一)

    2021-09-11 21:08:38
    我们在做优化工作过程中,往往面临着众多的不确定性,各种因素均会造成不确定,这使得传统运筹学中的确定性优化技术无法适应不确定的环境,那么比如说存在着下面几种造成不确定性的因素(实际中,会根据研究优化的...
  • 将逆向物流回收数量、再生产成本和市场需求作为不确定参数,以成本最小为目标,建立了包含生产成本、设备运作成本、库存成本在内的多周期多产品两阶段逆向物流网络鲁棒优化模型。通过算例验证了模型的有效性。
  • 综合考虑SNOP的功能特性与运行边界,建立了面向配电网弹性提升的SNOP配置三层防御-攻击-防御优化模型,并提出了求解该模型的两阶段鲁棒优化方法——列约束生成(CCG)算法。以IEEE 33节点为例,对所述模型和求解算法...
  • 数据驱动的鲁棒优化

    千次阅读 2020-01-16 14:52:33
    These data have motivated a shift in thinking—away from a priori reasoning and assumptions and towards a new data-centered paradigm. 这些数据促使人们从先验推理和假设转向以数据为中心的新范式。...
  • 通过仿真算例分析了多类型储能配置对综合能源系统规划问题的影响和不确定性参数的场景数量在阶段随机优化问题中的影响,并对比分析了阶段随机优化与两阶段鲁棒优化的优缺点。分析结果为综合能源系统灵活配置电/...
  • 对于这些误差目前有种手段来处理:随机优化和鲁棒优化。随机优化有一个很重要的假设:随机数据的分布是已知的或者是估计的。如果这个条件满足那么随机优化问题就是可行的。 鲁棒优化,并不需要随机数据的概率分布...
  • 然后,分析了鲁棒优化模型,随机优化模型和确定性优化模型之间的关系,并提出了种算法:枚举法和遗传算法。 这种算法的代码由Visual Studio 6.0上的Visual C ++实现。 代码中使用了优化软件Lingo 9.0来解决确定...
  • 针对工期不确定的资源受限项目调度问题,将鲁棒性资源分配和时间缓冲插入种方法进行有效地结合,通过设计两阶段集成优化算法构建抗干扰能力较强的鲁棒性项目调度计划....
  • 因此,提出了一种鲁棒两阶段动态多目标车辆路由方法。 该新方法的三个亮点是:(i)在第一阶段通过采用多目标粒子群优化为所有客户找到最佳的健壮虚拟路线后,通过从健壮的虚拟路线中删除所有动态客户来形成静态...
  • 已有的聚类算法大多仅考虑单一的目标,导致对某些形状的数据集性能较弱,为此提出一种基于改进粒子群优化的无标记数据鲁棒聚类算法。优化阶段:首先采用多目标粒子群优化的经典形式生成聚类解集合;然后使用K-means...
  • 说明列车到发信息、车流信息、调机作业时间信息等阶段计划编制依据的信息存在一定的不确定性,针对不同信息的不确定性对阶段计划编制和实施的影响,提出基于调度信息技术特征研究、阶段计划鲁棒优化方法和阶段计划...
  • 最新的matlab+yalmip+cplex综合能源及微电网优化运行研究相关论文,供研究学习参考。

空空如也

空空如也

1 2 3 4 5 ... 20
收藏数 13,818
精华内容 5,527
关键字:

两阶段鲁棒优化