精华内容
下载资源
问答
  • 两阶段鲁棒优化及约束生成算法
    万次阅读 多人点赞
    2020-10-18 21:06:47

    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

    更多相关内容
  • 本文尽量避免数学公式,使用文字解释列生成算法的原理,争取让读者能形成直观上的理解。 为什么需要了解列生成算法的原理 列生成算法无法简单地调用第三方库来使用,必须根据具体问题,构造不同的算法模型。 只有...

    本文尽量避免数学公式,使用文字解释列生成算法的原理,争取让读者能形成直观上的理解。

    为什么需要了解列生成算法的原理

    1. 列生成算法无法简单地调用第三方库来使用,必须根据具体问题,构造不同的算法模型。
    2. 只有了解了原理,才能在踩到各种坑时,有所针对地去优化各种细节。不然只能抓瞎或者抓腮。

    列生成算法原理

    列生成算法可以从两个视角来理解:对偶角度和单纯形算法角度。

    对偶角度

    啥是对偶

    这里简单过一下对偶的概念。

    假设有个长得很标准的线性规划问题:

    那么,它的对偶问题为:

    下面我们都以这个问题来讨论,即说到原问题时,默认是一个最小化问题;说到对偶问题时,默认是一个最大化问题。

    怎么理解这个对偶关系呢?借用经济学方面的话来说,假设原问题的目标是让成本最小,那么对偶就是让收入最大。更确切地讲,是:

    • 原问题丶:保证收入不低于某个值的条件下,使成本最小化。
    • 对偶问题:保证成本不高于某个值的条件下,使收入最大化。

    那个丶纯粹是为了对齐,忽略之……

    可以看到,原问题和对偶问题其实就是一个问题:目标净收益最大。只是一个是约束收入优化成本,一个是约束成本优化收入。角度不同而已。体现在公式上,就是原问题的变量对应对偶问题的约束,目标系数对应约束边界,约束矩阵倒转过来

    另外,关于对偶,一个比较重要的特性是:原问题的最优值与对偶问题的最优值相等

    从对偶角度看列生成算法

    列生成算法主要用途在于求解变量多,但是大部分变量将会取值为0的线性规划问题。总体思路是先忽略大部分变量,构造一个只使用小部分变量的模型(其余变量相当于值为0),这样就能很快求出一个解。然后寻找模型外的变量,找到能够让目标值更优的变量,加入模型再次求解。重复这个过程直至找不到更好的变量。

    这个过程的关键问题在于,怎么评估模型外的变量是否能让目标值更优。

    我们从对偶的角度来研究这个问题。

    原问题的变量对应对偶问题的约束。所以原问题新增变量,相当于对偶问题新增约束。

    原问题新增变量 -> 对偶问题新增约束

    由于对偶问题是个最大化问题,所以对偶问题新增约束后,显然最优值不变或变差,也就是不变或变小。从常理上看,约束越多,最优值越差嘛。

    而前面提到的,原问题的最优值等于对偶问题的最优值。也就是说,如果对偶问题最优值不变,那么原问题最优值也不变;如果对偶问题最优值变小,那么原问题最优值也变小。而我们需要的正是让原问题的最优值变小。

    所以问题变为如何尽量避免新增的约束没有改变最优值。设想一下,当加入新约束时,如果当前对偶的最优解没有违反新的约束,那么这个解仍然会是新增约束后的对偶问题的最优解,最优值将不变。

    因此,我们要找的新增的约束,要和当前最优解冲突

    整条逻辑链为:

    新增变量后原问题最优解变小 -> 新增约束后对偶问题最优解变小 -> 新增约束前的最优解不在新增约束后的可行域 -> 新增约束前的最优解不满足新增的约束

    一行对偶问题的约束的公式为:

    假设最优解为w*,那么违反约束的条件为:

    变换一下,变成:

    左侧的式子,叫做的reduced cost,也叫做检验数

    通过分析,我们知道,只要加入reduced cost小于0的对偶约束(从而加入了原问题对应的变量)即可

    很自然的想法是,我们更倾向于找到reduced cost最小的一个或几个变量加入,也就是最好能找到最小化reduced cost的新约束:

    这里就出现了一个新的最优化问题。这个问题叫做列生成的子问题(sub problem)。其中w*是已知的,未知量是c和a。c和a是和问题的应用场景有关的,需要根据实际场景来构造c和a的约束条件。所以子问题无法通用地求解,只能根据具体问题选择不同的方法求解。

    当所有未加入模型的变量的reduced cost都大于等于0时,目标值无法再优化,说明我们已得到最优解。

    另外,熟悉对偶问题经济学含义的同学会知道,reduced cost是指产品的差额成本。那么显然要新增的是差额成本为负的产品了。这是另一种理解列生成的思路。

    单纯形算法角度

    对偶角度给出了一个偏感性的方式来理解列生成算法。换个视角,从单纯形算法角度上看,则是单纯形算法本身,为了更高效地求解包含大量变量的问题,自然地扩展为列生成算法。

    相信有不少人被单纯形算法虐得有心理阴影——公式复杂,手工计算量也巨大……

    其实,如果我们先不看细节,单纯形的核心原理并没有那么难以理解。下面讲解时不会很严谨,理解算法框架就够了。严谨的过程请参阅运筹学相关书籍。

    单纯形算法

    众所周知,单纯形算法有一个几何上的解释:

    • 线性规划是一个凸优化问题,局部最优解就是全局最优解。
    • 线性规划的解空间是一个n维的凸多面体。最优解在这个凸多面体的某个顶点上。
    • 单纯形算法从一个初始顶点开始,不断沿着邻边找更好的顶点。
    • 当一个顶点四周没有更好的顶点时,这个顶点就是最优解。

    整个过程就像水沿着一条蜿蜒的沟渠流下,最终汇聚到最低点一样。

    问题是,这里面的几何概念和代数公式怎么对应?

    这里用不严谨(但更容易理解)的语言说明一下:

    • 边界:解空间是由不等式约束(包括变量非负这些约束)围起来的一块空间区域。当点p使得若干个不等式取等号时,那么点p就在约束边界的超平面上。这个边界可能是一个面、边、顶点。
    • 顶点:顶点会让尽量多的约束取等号。也就是说,顶点是由若干个改为等号的约束组成的方程组的解。我们叫这个方程组为约束边界方程组
    • 沿着边:约束边界方程组去掉一个方程,其解集就变成与顶点邻接的一条边。再取一条原方程组外的约束条件加入,所得到的解就是相邻的顶点。简单说,就是约束边界方程组中替换掉一个方程,形成的新方程组解出来就是相邻的顶点。

    这里涉及到通过让约束取等号来求边界的操作,而不等式乱糟糟地混在方程型的约束和变量非负约束里,会使这里的分析比较困难。所以使用单纯形算法之前,都会通过引入松弛变量、剩余变量和人工变量等方法(这一步在这里不重要,不详细展开了),将线性规划转换成如下标准形式:

    标准形式中只有变量非负约束包含不等式,其他约束都是等式。这样我们就可以很容易地做边界相关的计算了。假设变量数量为n,等式约束数量为m。通过转换而来的标准形式都会有n > m。那么,我们知道,只要让n-m个变量等于0,剩下的m个变量就可以通过这m个等式联立方程组(约束边界方程组)求出一个解(简单起见,不考虑无解,无数解这些边缘条件)。这个解就是一个顶点。

    这里约束边界方程组中的m个变量叫作基变量,固定值为0的n-m个变量叫作非基变量

    沿着边找相邻顶点,就是取一个被固定为0的非负约束,也就是一个非基变量(这个操作叫入基),替换掉一个基变量(这个操作叫出基,这个变量出基后就固定值为0),然后重新求解一个顶点。

    入基操作需要选择入基变量,选择的依据是这个变量在目标函数中的下降速度,也就是这个变量增加1时,目标值减少多少。经过推导可知,下降速度的计算公式刚好是检验数(reduced cost)。这里就和对偶的视角联系起来了。

    出基操作这里就不细说了,大致的思路是在约束条件下,旧的基变量有一部分会随着入基变量的增长而下降,其中最先下降到0的旧的基变量就会被选为出基变量。

    整个单纯形算法的计算步骤是:

    1. 选取基变量和非基变量,简单能出初始解就好。
    2. 计算所有非基变量的reduced cost,找到最小且为负值的那个作为入基变量。如果reduced cost都大于等于0,迭代终止。
    3. 选出基变量
    4. 解约束边界方程组,回到步骤2

    从单纯形算法角度看列生成算法

    在单纯形的步骤2,需要计算所有非基变量的RC。找到最小的那个。当变量个数很多的时候,这一步就成为了算法运行时间瓶颈。

    在一些情况下,通过巧妙构造问题,可以让这一步不需要遍历所有变量。甚至我们都不需要知道有多少变量,只要能在每次迭代的时候生成一个或者多个变量,提升优化效果就可以了。

    由于不需要遍历所有变量,所以一开始就不需要使用所有变量,只需要使用一组能产生初始解的初始变量构成线性优化问题即可。这种只使用部分变量的模型被称为原问题的restricted master problem(RMP)

    每次迭代时,生成一个或多个让reduced cost最小的变量加入RMP。这个生成步骤就是求解子问题。不断加入新变量直到没有小于0的reduced cost的变量时就达到最优解。

    到这里就和对偶角度分析的结果一致了。

    下面是单纯形算法与列生成算法简要流程图的对比,可以看到,两者的结构是一样的。

    一般来说,我们不会手搓单纯形算法,所以正常都是直接调用单纯形算法库解RMP,然后做列生成,再跑RMP,直到达到最优。

    一个经典例子:Cutting Stock Problem

    这是一个列生成算法的经典例子。

    原纸卷每个长17m,顾客们分别需要25个3m长,20个5m长,18个7m长的纸卷。
    问:如何切割使消耗的原纸卷数量最少?

    令一个原纸卷的切割方案集合为:

    P = {(a, b, c) | 3a + 5b + 7c <= 17}

    其中,a是一个原纸卷切割出的3m纸卷数量,b是5m纸卷数量,c是7m纸卷数量。

    我们用变量x(abc)表示使用切割方案(a, b, c)的原纸卷数量。

    显然,一个变量与一个原纸卷切割方案一一对应。建模如下:

    这里故意不适用传统的下标序号标记,意在突出我们不需要对变量编号,只需要知道变量在对应在什么集合上,如何通过集合中的元素生成变量就行了。

    初始解很好找。比如说我们可以取25个原纸卷按照方案(1, 0, 0)切割,20个原纸卷按照方案(0, 1, 0)切割,18个原纸卷按照方案(0, 0, 1)切割。这当然会有很多浪费。但是初始解可行就可以了,浪费的部分会在下面的迭代中优化掉。

    接下来要生成变量。变量与切割方案一一对应的。所以是要找出一个切割方案(a, b, c),使得reduced cost最小。

    其中w1、w2和w3分别为约束R1、R2和R3的对偶值。

    约束条件除了a、b、c非负外,还需要满足切割后的纸卷长度综合小于或者等于原纸卷的长度。

    这样子问题就构造好了。求解子问题得到新增变量。然后迭代直到最优。具体计算这里不展开了。

    整数规划求解

    前面提到的单纯形算法和列生成算法求解的都是线性规划。在实际应用中,一般还会需要求解整数规划。也就是变量都约束为整数的线性规划。

    这里先提一个概念:整数规划的线性松弛,整数规划问题,不考虑整数约束,剩下的约束条件和目标组成的线性规划问题。

    其实我们并没有很好的方法直接求解整数规划,通常都是不断地调整并求解线性松弛,最后找到最优整数解。

    分支定界

    分支定界是一个用来求整数优化问题的框架。其实思路很简单,就是采用类似二分法的技巧,在线性解空间中暴力搜索整数解。

    首先,求解线性松弛得到线性解。取一个变量x进行分支。比如x线性解值为1.2,那么产生 x <= 1 和 x >= 2 两个分支。将这两个条件分别加入到线性松弛,得到两个线性规划。再求解这两个线性规划,两个又分别分支……直到求得最优解。

    有些情况下判断一个解是否为最优解是有方法的,所以不用搜索所有分支。但是,分支定界在最坏情况下的时间复杂度仍是指数级。为了防止运行时间过长,一般使用分支定界时还会额外加一些终止条件,比如回溯次数限制、运行时间限制、找到第一个整数解就结束等。

    分支定价

    分支定价就是在分支定界框架中,使用列生成算法来求解每个分支节点的算法。

    不过这里,除了根节点,其他节点不用从头开始生成新变量,继承父节点用到的变量即可,这样可以节省很多重复生成变量的过程。

    在01整数规划中,还有更简化的方法,每次列生成得到线性松弛的最优解后,找出值最接近1的变量,新加这些变量等于1的约束,继续跑列生成,直到找出所以值为1的变量。剩下的自然都是0了。这个方法可以加入回溯,也可以不回溯,出现无解就直接结束……据说不回溯也很少出现无解的情况。

    一些相关问题

    退化问题/类退化问题

    通过RC找的新变量不一定能让目标值变得更好,仍然存在不变的可能。极小概率的情况下,单纯形算法可能会有入基变量和出基变量循环出现的情况。由于我们肯定是调用线性规划库来跑单纯形的,所以不用考虑这个……

    列生成没有出基操作,不会出现循环。但是有一些改进会剔除冗余变量,这时就会有极小概率会出现循环了。这种情况不需要费心去处理,玄学调参降低出现概率,并设置最大迭代次数等强制终止条件,确保能终止就好。

    最恶心的情况是没有循环,但是长时间没有提升目标值的情况。这其实是算法卡在一个拐点上了,只要过了拐点就能开始提升。特别是在一些约束较强的问题(比如密集的排班问题)中,使用某些启发式算法或者手工做出来的初始解就很容易出现这种情况。

    而我们为了避免算法跑太久,通常会设置多次迭代没有提升就结束的条件。这可能使算法从拐点出发后,几次迭代无优化就直接结束了。

    这种情况无法完美解决。简单的就是调参加更巧妙地设置结束条件,通过多次试验尽量让算法能跨过这个拐点。还有另一个技巧是可以适当地给约束边界加一下噪音,比如说小于等于1的约束,可以放宽到小于等于1.0001。这样从初始解出发迭代时,由于边界宽松了一些,变量可以有些许变化,会让目标值有一些微小的提升,帮助判断是否需要结束迭代。

    CPLEX计算reduced cost的问题

    使用CPLEX时,我们可以很方便的设置变量的上下界。比如设置 0 <= x <= 1。这时,x <= 1这部分是会影响reduced cost的值的。而CPLEX接口计算的是没有考虑这个条件的……所以可能你自己手搓代码出来reduced cost和CPLEX接口出来的reduced cost不相等。

    更严重的是,可能你会忘了 x <= 1 这个条件,导致列生成的过程中算错reduced cost。

    这个问题其实影响不大,主要是会干扰一些计算过程正确性的验算。

    如何验证子问题有没有严重问题

    没有做分支操作时,线性松弛的目标值如果变差,说明子问题可能出现了一些很蠢的问题。

    线性规划求解算法选择

    每次迭代求解线性规划时,选用不同的算法会影响求解时间。根据经验:

    • 增加少量列时(列生成),使用单纯形算法(Primal)。
    • 增加少量约束时(分支),使用对偶单纯形算法(Dual)。
    • 其他情况,酌情使用对偶单纯形算法或者内点法。通过试验决定。
      • 对偶单纯形算法快的时候很快,慢的时候很慢。
      • 内点法速度比较稳定。

    以上也并不是所有场景通用的。应当针对具体问题,反复试验来确定使用什么算法。

    迭代中剔除冗余变量

    Reduced cost可以用来评估变量的“有用”程度。越小表示变量越有用,越大表示变量越无关紧要。

    列生成迭代次数较多后,变量数量会越来越多,从而每次迭代的运行速度越来越低。可以设定一个变量规模上限,当变量数量大于上限时,从模型中去掉reduced cost最大的那些变量。

    标签: 整数规划对偶单纯形线性规划运筹学分支定界分支定价列生成

    展开全文
  • 本期小编为大家介绍单纯形法的另一种拓展算法—列生成算法,主要从其产生背景、基本思想、应用实例、平台实现以及与单纯形法的区别和联系五个方面展开讲解,一起来看一看吧! 列生成算法(Column Generation ...

    本期小编为大家介绍单纯形法的另一种拓展算法—列生成算法,主要从其产生背景、基本思想、应用实例、平台实现以及与单纯形法的区别和联系五个方面展开讲解,一起来看一看吧!

    列生成算法(Column Generation Algorithm)是一种用于求解大规模线性优化问题的高效算法,其理论基础由Danzig等人于1960年提出。本质上来讲,列生成算法是单纯形法的一种形式,用来求解线性规划问题。

    一、产生背景

    在某些线性规划问题的模型中,约束条件的数目有限,但变量的数目会随着问题规模的增长爆炸式的增长,很难把所有的变量都显性的在模型中表达出来,这类问题就是大规模的线性规划问题。面对这类问题,单纯形法虽然能保证在数次迭代后找到最优解,但由于其需要对众多变量进行基变换,求解过程会异常繁琐。此外,在用单纯形法求解这类问题时,基变量只与约束条件的个数相关,每次迭代只会有一个新的非基变量进基,换言之,在整个求解过程中其实只有很少一部分变量会被涉及到。在这种背景下,研究人员基于单纯形法提出了列生成算法。

    该算法不是直接同时处理所有的候选方案,而是基于当前生成的列的子集,通过限制主问题进行优化求解;其余的候选方案可以改善限制主问题当前最优解时,才会进入该子集。与单纯形法相比,列生成的进基变量是通过求解子问题生成的,而单纯形法的进基变量是模型存在的变量。

    目前,列生成算法已被应用于求解很多著名的NP-hard优化问题,如机组人员调度问题(Crew assignment problem)、切割问题(Cutting stock problem)、车辆路径问题(Vehicle routing problem)、单资源工厂选址问题(The single facility location problem)等。

    二、基本思想

    列生成算法通过求解子问题(sub problem)来找到可以进基的非基变量,该非基变量在模型中并没有显性的写出来(可以看成是生成了一个变量,每个变量其实等价于一列,所以该方法被称为列生成算法)。如果找不到一个可以进基的非基变量,那么就意味着所有的非基变量的检验数(Reduced Cost,RC)都满足最优解的条件,也就是说,该线性规划的最优解已被找到。其思路如下:

    1、先把原问题(Master Problem,MP)转换到一个规模更小(即变量数比原问题少)的问题上,这个只使用部分变量的模型被称为原问题的RMP 问题(Restricted Master Problem)。在RMP上用单纯形法求最优解,注意此时求得的最优解只是RMP上的,并不是MP的最优解。

    2、然后需要通过一个子问题去检测在那些未被考虑的变量中是否有使得RC小于零的情况,如果有,那么就把这个变量的相关系数列加入到RMP的系数矩阵中,返回第1步。

    经过反复迭代,直到子问题中的RC都大于等于零,此时就找到了MP的最优解。流程图如下:

    三、应用实例

    1、问题描述

    有三种长度分别为91416的木材,对应的成本价为5910,需要切割成长度为4的成品30个;长度为5的成品20个;长度为7的成品40个。求解切割方案,使得总体成本价最低。

    首先枚举所有可能切法:

     

    (2)建立模型:定义变量x_{j}表示使用第j种切割方案的次数,c_{j}表示第j种方案的成本,a_{ij}表示第j种方案生成第i种木材的数量(i=1,2,3),所有切割方案的集合为J。a_{j}=(a_{1j},a_{2j},a_{3j})^{^{T}}表示j方案下分别切出来尺寸为4、5、7的木材的个数,b=(30,20,40)^{T}表示所需尺寸为4、5、7木材的最小个数。该问题表示如下:

     (1)求解目标:找到最好的方案组合,使价格最小。

    2、列生成方法:

    (1)初始化

    在枚举出的14种方案中挑选较好(价格最低)的方案组合,即前三种方案:

    (2)主问题

    此时的主问题为:

    基于初始化的方案求解,利用gurobi求出最佳方案组合:x^{*}=(5,20,40)^{T},z^{*}=325,w^{*}=(2.5,2.5,5),即前三种方案分别使用15、20、40次,总成本为325, w^{*}表示生成这三种长度的木材的合理价格。

    (3)子问题

    节省成本:

     即找到一个新的方案,计算该方案较原始方案的节省成本(切出来的总价值-成本)。当节省的成本都小于等于0时,求解问题达到最优。由于原木材有三种:长度尺寸为9、14、16,所以可分解为三个子问题。

     子问题2求解出来的节省成本最大,加入到主问题中。

    (4)迭代

    求解新生成的主问题,直到主问题对应的三个子问题求解出来的节省成本都是负数,表明当前主问题的解已是最优解,停止计算并输出结果。

    四、平台实现

    基于以上求解步骤即思想,依靠Python语言,借助Pycharm平台和gurobi求解器,求解上述板材切割问题。部分代码如下(关注“运筹学”公众号→后台回复“列生成算法”即可获取完整代码):

    最终结果:

    可以得到原问题的最优解为:x_{1}=5,x_{2}=20,x_{4}=20,z^{*}=305,即采用第一种、第二种和第四种切割方案,使用次数分别为5、20、20,最小成本为305。

    五、与单纯形法对比

    1、区别

    1在求解过程上:单纯形法需要计算所有非基变量的RC,而列生成算法通过巧妙构造子问题,让这一步不需要遍历所有变量,甚至都不需要知道一共有多少变量,只要能在每次迭代的时候生成一个或者多个变量,提升优化效果就可以了,通过不断加入新变量直到没有小于0RC的变量来求得最优解。

    2在适用范围上:与单纯形法相比,列生成算法可以求解变量很多的线性规划问题,也就是前面所说的大规模的线性规划问题。

    2、联系

    从本质上来讲,列生成算法就是单纯形法的一种形式,是用来求解线性规划问题的。此外,列生成算法在求解每个子问题时,需要用到单纯形法。

    ★ 在python中安装Gurobi的详细教程可以参考以下网址:

        https://blog.csdn.net/weixin_41596280/article/details/89112302

    ★ 列生成算法参考资料1:https://blog.csdn.net/Rivalsx/article/details/97448943

    ★ 列生成算法参考资料2:https://blog.csdn.net/u014090505/article/details/89019327


    作者|陈怡敏、曹贵玲
    责编|何洋洋
    审核|徐小峰 

    展开全文
  • 列生成算法原理(单纯形基础)

    千次阅读 2021-09-22 20:29:29
    转载自:【学界】单纯形算法的原理和示例实现 一 线性规划问题的标准形式 线性规划问题的标准型: 其中 一般总假设的元素都为整数(线性规划多应用在实际问题,所以要求为整数也可以理解), 可行域: 最...

    一 、单纯形法基础

     1.1 定理和推论

    定理1: X 是可行域 P 的一个顶点 \Rightarrow x 的正分量对应的 A 中的各列是线性独立的。

    补充:

    (1)

     

     (2)

     *******推论*******:

    注: 基可行解对应有解,即列向量之间线性无关(定理1成立)

     满秩

    由于可能有退化的情况,所以基可行解和顶点不是一一映射

    任给一个基可行解,存在唯一的一个顶点与之对应;对于 P 的一个顶点,可能有多个基可行解与之对应

    1.2 单纯形算法

     

    则上述线性规划问题就变成:

     

     

     注:1. 检验数为正,说明 reduce cost可以继续下降 ;2. 当所有检验数为正值,reduce cost对于所有为正的x只能继续上升,无法下降,因此达到了最优解

    1.3 求基可行解方法

                                                                 原始问题标准型

                                                           对约束构造人工变量

    此时约束已经变了,目标问题求解也就变了。两种通用的方法就出来了---大M法和两阶段法,一般使用大M法,M是一个充分大的数,目标函数求最大就取减号,求最小就取加号。在这个问题中要有最优解,人工变量就必须取0,(大M作为惩罚因子使得最优解对应的变量一定为“0”,从而消除加入变量的影响)否则,原问题无解。这样就很顺利的得到了一个基可行解,然后代入单纯形算法进行迭代

    1.4 入基变量的选取:

    最简单的选取规则就是:目标函数求max,检验数大的为入基变量,目标函数求min,检验数小的为入基变量

    1.5 出基变量的选取:

    最小比值法选取出基变量,当选取完入基变量  \lambda_{m+k}后,取 x_{m+k}=min\frac{a_i}{\beta_{i,m+k}} (将出基变量变为“0”,从而得到入基变量的值,),相应的 x_i 做为出基变量置为0,取最小的原则是可以保证 [公式] ,防止出现非负变量。(不取最小值,则会使其他变量变为负数,从而违背变量x全为正的约束条件

    借用原作者的图,在进行出基、入基后,使用新基变量B^{-1}使得基变量系数为单位阵(表示线性独立的形式,即一个基变量均由非基变量和常数项b来表示),非基变量系数则使用B^{-1}N进行更新,从而重新表达新的变量之间的关系

    1.6 几何解释

    • 线性规划是一个凸优化问题,局部最优解就是全局最优解。
    • 线性规划的解空间是一个n维的凸多面体。最优解在这个凸多面体的某个顶点上。
    • 单纯形算法从一个初始顶点开始,不断沿着邻边找更好的顶点。
    • 当一个顶点四周没有更好的顶点时,这个顶点就是最优解。

    二、列生成

      单纯形需要计算所有非基变量的RC。找到最小的那个。当变量个数很多的时候,这一步就成为了算法运行时间瓶颈。

    而列生成只要能在每次迭代的时候生成一个或者多个变量,提升优化效果就可以了。由于不需要遍历所有变量,所以一开始就不需要使用所有变量,只需要使用一组能产生初始解的初始变量构成线性优化问题即可。这种只使用部分变量的模型被称为原问题的restricted master problem(RMP)

    每次迭代时,生成一个或多个让reduced cost最小的变量加入RMP。这个生成步骤就是求解子问题。不断加入新变量直到没有小于0的reduced cost的变量时就达到最优解。

    下面是单纯形算法与列生成算法简要流程图的对比,可以看到,两者的结构是一样的。

    2.1  对偶角度

    假设有个很标准的线性规划问题:

    那么,它的对偶问题为:

    借用经济学方面的话来说,假设原问题的目标是让成本最小,那么对偶就是让收入最大。更确切地讲,是:

    • 原问题   :保证收入不低于某个值的条件下,使成本最小化。
    • 对偶问题:保证成本不高于某个值的条件下,使收入最大化。

    可以看到,原问题和对偶问题其实就是一个问题:目标净收益最大。只是一个是约束收入优化成本,一个是约束成本优化收入。角度不同而已。体现在公式上,就是原问题的变量对应对偶问题的约束目标系数对应约束边界,约束矩阵倒转过来

    另外,关于对偶,一个比较重要的特性是:原问题的最优值与对偶问题的最优值相等

     从对偶角度看列生成算法

    列生成算法主要用途在于求解变量多,但是大部分变量将会取值为0的线性规划问题。总体思路是先忽略大部分变量,构造一个只使用小部分变量的模型(其余变量相当于值为0),这样就能很快求出一个解。然后寻找模型外的变量,找到能够让目标值更优的变量,加入模型再次求解。重复这个过程直至找不到更好的变量。

    这个过程的关键问题在于,怎么评估模型外的变量是否能让目标值更优

    我们从对偶的角度来研究这个问题。原问题的变量对应对偶问题的约束。所以原问题新增变量,相当于对偶问题新增约束。

    原问题新增变量 -> 对偶问题新增约束

    由于对偶问题是个最大化问题,所以对偶问题新增约束后,显然最优值不变或变差,也就是不变或变小。(约束越多限制越多,结果会越来越差,最好的结果仅仅是保持不变)

    而原问题的最优值等于对偶问题的最优值。也就是说,如果对偶问题最优值不变,那么原问题最优值也不变;如果对偶问题最优值变小,那么原问题最优值也变小。而我们需要的正是让原问题的最优值变小。

    所以问题变为如何尽量避免新增的约束没有改变最优值。设想一下,当加入新约束时,如果当前对偶的最优解没有违反新的约束,那么这个解仍然会是新增约束后的对偶问题的最优解,最优值将不变

    因此,我们要找的新增的约束,要和当前最优解冲突

    ***************

    新增变量后原问题最优解变小 -> 新增约束后对偶问题最优解变小 -> 新增约束前的最优解不在新增约束后的可行域 -> 新增约束前的最优解不满足新增的约束

    *****************

    一行对偶问题的约束的公式为:

    假设最优解为w*,那么违反约束的条件为:

    变换一下,变成:

    左侧的式子,叫做的reduced cost,也叫做检验数

    通过分析,我们知道,只要加入reduced cost小于0的对偶约束(从而加入了原问题对应的变量)即可

    我们更倾向于找到reduced cost最小的一个或几个变量加入,也就是最好能找到最小化reduced cost的新约束:

    这里就出现了一个新的最优化问题。这个问题叫做列生成的子问题(sub problem)。其中w*是已知的,未知量是c和a。c和a是和问题的应用场景有关的,需要根据实际场景来构造c和a的约束条件。所以子问题无法通用地求解,只能根据具体问题选择不同的方法求解

    当所有未加入模型的变量的reduced cost都大于等于0时,目标值无法再优化,说明我们已得到最优解

     

    3. Cutting Stock Problem

    原纸卷每个长17m,顾客们分别需要25个3m长,20个5m长,18个7m长的纸卷。
    问:如何切割使消耗的原纸卷数量最少?

    令一个原纸卷的切割方案集合为:

    P = {(a, b, c) | 3a + 5b + 7c <= 17}

    其中,a是一个原纸卷切割出的3m纸卷数量,b5m纸卷数量,c7m纸卷数量。

    我们用变量x(abc)表示使用切割方案(a, b, c)的原纸卷数量。

    显然,一个变量与一个原纸卷切割方案一一对应。建模如下:

     

    这里故意不使用传统的下标序号标记,意在突出我们不需要对变量编号,只需要知道变量在对应在什么集合上,如何通过集合中的元素生成变量就行了。

    初始解很好找。比如说我们可以取25个原纸卷按照方案(1, 0, 0)切割,20个原纸卷按照方案(0, 1, 0)切割,18个原纸卷按照方案(0, 0, 1)切割。这当然会有很多浪费。但是初始解可行就可以了,浪费的部分会在下面的迭代中优化掉。

    接下来要生成变量。变量与切割方案一一对应的。所以是要找出一个切割方案(a, b, c),使得reduced cost最小。

     其中w1、w2和w3分别为约束R1、R2和R3的对偶值。(RC 为\sigma_j=c_j-c_BB^{-1}a_j其中c为1

    c_BB^{-1}有两重含义:

    • 通过求解RMP问题得到的影子价格(shadow price)。
    • 通过求解RMP对偶问题得到的对偶变量(dual variable)。

    约束条件除了a、b、c非负外,还需要满足切割后的纸卷长度综合小于或者等于原纸卷的长度。

    这样子问题就构造好了。求解子问题得到新增变量。然后迭代直到最优。

    4、应用实例

    • P是所有可行的裁剪方案的集合,里面方案的总数为n(我们并不需要确切的知道这个值是多少,只需要知道它很大)。
    • a_{ij}表示第j种方案里类别i的个数。
    • y_j表示第 j 种方案的选择个数。

    于是,我们得到如下模型:

     使用LpsolveIDE,输入为:

                    ​​​​​​​        ​​​​​​​        ​​​​​​​        ​​​​​​​        

    首先,一个卷筒有三种切割方案:
    方案1:切成5个3m
    方案2:切成2个6m
    方案3:切成2个7m

    很容易得出,5个方案1、10个方案2、8个方案3,是能满足所有客户需求的。即得到MP的一个RMP(restricted master problem)如下:

                                                     

    其中:

                                                            

    这三列分别对应着5个方案1、10个方案2、8个方案3。还有一点需要注意的,对于每一列,都需要满足:

                                            3a_{1j}+6a{2j}+7 a{3j}\leq16

    也就是每一卷纸只有16的长度,不能超出这个长度。这个叫列生成规则,不同问题有不同的规则约束。

    正式迭代

    iteration 1

    RMP:

                                             

    将该模型输入lpsolve,得到对偶变量如下:

                             

     得到    c_BB^{-1}=[0.2,0.5,0.5]]   在要找一列加入RMP,是哪一列呢?现在还不知道,我们暂记为   a_4=[a_{14},a_{24},a_{34}]^T  。

    非基变量检验数       \sigma_4=c_4-c_BB^{-1}a_4=1-0.2a_{14}-0.5a_{24}-0.5a_{34}

    子问题:

                                    

    求解结果得  \alpha_4 = [1,2,0]^T, \sigma_4= -0.2 < 0      reducedcost为负数,因此将  \alpha_4 加入RMP,开始第二轮迭代。

    iteration 2

     RMP:

                                             

    将该模型输入lpsolve,得到对偶变量如下:

            ​​​​​​​        

    得到c_BB^{-1}=[0.2,0.4,0.5]]   ,现在要找一列加入RMP,我们暂记为a_5=[a_{15},a_{25},a_{35}]^T

     非基变量检验数       \sigma_5=c_5-c_BB^{-1}a_5=1-0.2a_{15}-0.4a_{25}-0.5a_{35}

    子问题:

                                     

    求解结果得\alpha_5 = [1,1,1]^T, \sigma_5= -0.1 < 0,  reducedcost为负数,因此将 \alpha_5加入RMP,开始第三轮迭代。

    iteration 3

    RMP:

                                     

     将该模型输入lpsolve,得到对偶变量如下:

                             

    得到c_BB^{-1}=[0.2,0.4,0.4]],现在要找一列加入RMP,我们暂记为a_6=[a_{16},a_{26},a_{36}]^T

     非基变量检验数    \sigma_6=c_6-c_BB^{-1}a_6=1-0.2a_{16}-0.4a_{26}-0.4a_{36}

    subproblem:

            ​​​​​​​        ​​​​​​​        

    求解结果得 \alpha_6 = [5,0,0]^T, \sigma_6 = 0,reducedcost不为负数,因此不用将 \alpha_6 加入RMP,列生成算法结束。

    最终,我们求解最后一次迭代的RMP:

                                    

            ​​​​​​​        ​​​​​​​        ​​​​​​​        ​​​​​​​        

    得到RM的最优解   y=[1.2,0,0,1,18] , 按理说y应该是整数才对。回到原问题

    RM:

            ​​​​​​​        ​​​​​​​        ​​​​​​​        ​​​​​​​        

    我们并没有加上y_i\in Z 这个约束,这是因为我们在用列生成的时候把这个模型给松弛为了线性模型,毕竟列生成是用于求解linear program的。如果要求解大规模整数规划问题,列生成是无法办到的,后面我们会介绍结合column generation的branch and price方法。

    限制主问题和子问题的关系为:
    在这里插入图片描述​​​​​​​ 

     

    5、列生成和分支定价

     在这里插入图片描述

     

    Reference:

    【1】【学界】单纯形算法的原理和示例实现 

    【2】  深入浅出列生成算法  

    【3】 算法介绍之列生成算法   

    【4】彻底了解column generation(列生成)算法的原理

    展开全文
  • 这篇博文主要是想记录一下我这几天的工作心得。...标题中的约束生成或者行生成是我自己翻译的,可能存在误差(单纯形法中用一行表示一个约束)。前人提出这个算法的主要原因是因为有的问题约束特别多,例...
  • 包含十大经典算法: 顶点覆盖近似算法、哈密尔顿回路、画等温线、模拟退火应用、生成全排列矩阵、随机数的产生、最大流和最小截、最短路和次短路、最短路径、最小生成树Prim算法
  • MATLAB代码:基于列约束生成法CCG的两阶段问题求解 ...主要内容:代码构建了两阶段鲁棒优化模型,并用文档中的相对简单的算例,进行CCG算法的验证,此篇文献是CCG算法或者列约束生成算法的入门级文献,其经典程度不
  • 根据问题特点,使用Danzig-Wolfe策略将这个模型分解为一个带有集划分约束的主问题和一个具有背包特征约束的价格子问题,开发了分支价格算法进行求解。计算结果表明所开发的分支价格算法能够最优求解生产实际问题。
  • 列生成column generation算法,是求解大规模线性规划问题的一种常用算法,在运筹优化领域有着非常广泛的应用。 本篇介绍它的最初形式,原理及应用。 1.最初是为了解决下料问题 即cut-stock problem。这个问题在cplex...
  • 列生成(Column Generation)算法

    千次阅读 2021-05-11 20:57:49
    列生成算法的背景 多年来,寻找大规模的、复杂的优化问题的最优解一直是决策优化领域重要的研究方向之一。 列生成算法通常被应用于求解大规模整数规划问题的分支定价算法(branch-and-price algorithm)中,其理论...
  • 列生成算法是求解大规模线性规划的有效手段之一,对于难以直接求解的大规模线性规划问题,生成在某些场景能够缩小问题规模,降低问题求解难度,在可接受的时间内精确求解问题。生成的基本思路是将一个原问题从一...
  • 讲解列生成算法的标配问题。可以百度,一百度一大把。 模型建立 代码实现 ref:改代码引用自gurobi官方示例。 个人感觉这个实现结构还是非常优秀的。 from gurobipy import * TypesDemand = [3, 7, 9, 16] # ...
  • 利用Python+Gurobi编写代码,复现文章:Solving two-stage robust optimization problems using a column-and- constraint generation method。
  • 浅析constraint generation(约束生成,行生成)和column generation(生成) 然而当两个算法需要统一成一个算法(Column-and-row Generation,以下简称CRG)的时候,发现相关的中资料特别少,那我就来抛砖引玉吧...
  • 优化|列生成算法及Java调用cplex实现

    千次阅读 2020-12-20 11:47:31
    优化|列生成算法及Java调用cplex实现Cutting Stock ProblemColumn Generation AlgorithmJava调用cplex实现CG算法 Cutting Stock Problem 本文中的课件来自清华大学深圳国际研究生院,物流与交通学部张灿荣教授...
  • 列生成算法的学习

    千次阅读 2019-11-04 10:14:44
    因此,一般软件用的是barrier method,也就是内点法求解,这个是多项式时间算法,这个算法好像可以很快求解得到最优解的。注意,是好像,没有验证过。 2、继续问题。如果是整数规划的大规模线性规划问题,还是要...
  • 从一维cutting问题看列生成算法

    千次阅读 2019-08-13 23:21:08
    从一维cutting问题看列生成算法列生成算法一维cutting问题新的改变功能快捷键合理的创建标题,有助于目录的生成如何改变文本的样式插入链接与图片如何插入一段漂亮的代码片生成一个适合你的列表创建一个表格设定内容...
  • cutting stock 问题; 原始整数规划模型; 列生成算法
  • OUTLINE 前言 预备知识预警 什么是column generation 相关概念科普 Cutting Stock Problem CG求解Cutting Stock ...到列生成算法这一块,看了好几天总算把这块硬骨头给啃下来了。然后发现网上关于生成的教...
  • 论文研究-多星联合对地观测调度问题的列生成算法.pdf, 多星联合对地观测调度问题作为一类大规模组合优化问题, 其求解算法往往采用启发式或超启发式. 运用生成思想对该...
  • YALMIP+IBM CPLEX12.8 matlab toolbox 支持包亲测可用,解压后添加到matlab启动路径即可使用。yalmiptest命令查看支持情况,内有测试程序。
  • 目录 1 CSP问题与模型 2 列生成方法理论 3 Cplex OPL演示列生成迭代过程 4 多种长度木材的例子 5 java实现代码
  • 列生成和分支定价

    万次阅读 热门讨论 2019-04-04 12:14:03
    列生成算法适用于求解一类每个决策方案对应整体规划模型中约束矩阵的一的组合优化问题。该算法不是直接同时处理所有的候选方案,而是基于当前生成的的子集,通过限制主问题进行优化求解;其余的候选方案可以改善...
  • 为进一步改善列生成算法存在的尾效应, 将基于次梯度优化的拉格朗日松弛算法嵌入列生成算法框架中, 构建了采用双重迭代的改进型生成(MCG) 算法. 最后, 通过理论分析和仿真实验表明了MCG算法是有效、可行的.</p>
  • 基于C&CG列约束生成算法和强对偶理论,可将原问题分解为具有混合整数线性特征的主问题和子问题进行交替求解,从而得到原问题的最优解。程序基于MATLAB yalmip调用CPLEX实现求解,每一行代码均有详细注释,同时,附带...
  • 二阶随机优化算法小结

    千次阅读 2018-07-07 18:48:10
    二阶随机优化算法 标签(空格分隔): L-BFGS NewSample Lissa 最近看了几篇二阶优化算法,现在总结一下,以便日后查阅 二阶算法 二阶优化算法又称为牛顿法,牛顿法是微积分学中, 通过迭代以求解可微函数f...
  • 基于Cplex的列生成技术

    千次阅读 2021-01-15 17:15:36
    首先介绍一下什么是列生成列生成技术是一种求解大型线性规划问题的技术,广泛应用于机组人员调度问题、切割问题、车辆路径规划问题。以上问题有一个共同点,即小规模的问题可以直接使用线性规划进行求解,但是对于...
  • 该方法的特点是 在每一步列生成算法的迭代中,通过机器学习训练出一个模型来生成新的(选择一个变量的子集添加到主问题构成新的一)。该方法通过选择最合适的来降低主问题的计算时间。我们在带有时间窗约束的...
  • 列生成(Column Generation)是一种把线性规划问题分解为小规模子问题的技巧1 2. 它的原理基于单纯形算法. 从一个基本解(Basic Solution)出发, 主问题的系数矩阵只考虑基本解对应的, 然后求解子问题来生成主问题系数...
  • 为了改善航班计划两阶段完成的次优性,将机型指派、路线选择与机尾号指派综合考虑,构建了...通过设计了一种基于约束编程思想的列生成算法对该问题进行求解。最后,通过航空公司实例数据验证了模型算法的正确性和有效性。

空空如也

空空如也

1 2 3 4 5 ... 20
收藏数 43,308
精华内容 17,323
关键字:

列约束生成算法