运筹系列8：Set partitioning问题的列生成方法

IE06

已于 2023-09-24 20:43:56 修改

阅读量1.7w

点赞数 19

分类专栏：运筹学文章标签：算法机器学习动态规划

于 2018-08-15 22:32:47 首次发布

本文链接：https://blog.csdn.net/kittyzc/article/details/81712448

版权

运筹学专栏收录该内容

93 篇文章 350 订阅

订阅专栏

1. 应用场景

列生成算法是不断添加求解变量的算法，可参考论文。上一节我们介绍了约束条件存在分块结构的情形下可以使用单纯型法的判断规则进行列生成，下面我们介绍另一种场景，set-partitionning问题的列生成算法。
具体来说，场景如下：有 $n$ 个0-1变量 $z_1...z_n$ ，每个 $z_i$ 带着很多中间变量 $x_{i,j}$ 用来进行约束，这是个变量很多，约束也很多的模型。我们首先对问题模型进行转化，将所有 $z_i$ 的组合枚举出来，在枚举的过程中，就把约束条件都过了一遍，只有满足所有约束条件的组合才会保留下来。最后的枚举结果分别对应 $y_1...y_m$ 。
不难看出， $m$ 是 $n$ 的指数形式，使用set-partitioning模型之后，变量的数量大的可怕。好消息是，原问题的约束条件比较紧，只有少数组合成立，可以通过一定的规则产生，这样m的数量比较小，问题就比较好解决。

2. 例1：cutting stock problem

来看网上流传很多的cutting stock problem的例子，问题如下：
我们有一些纸筒，每个纸筒16m。顾客需要裁剪的长度为：25个3m，20个6m，15个7m。要求在满足顾客需求的情况下，裁剪的纸筒数最小。
我们可以枚举出每筒可能的裁剪方案：
方案1：5个3m，即3+3+3+3+3
方案2：2个6m，即6+6
方案3：2个7m，即7+7
我们可以继续枚举下去：
4：6+6+3
5：6+3+3+3
6：7+3+3+3
7：7+6+3
…

仅使用前三种方案，我们用23筒（5个方案1，10个方案2，8个方案3）是肯定可以满足要求的，这算问题的一个上界。下面我们探索一下其他的切筒方式，看能不能给出下界。

一般性的，我们假设 $P$ 是所有可行的裁剪方案的集合（其中前3个可以设置为我们前面的3个裁剪方案），里面方案的总数为 $n$ (我们并不需要确切的知道这个值是多少，只需要知道它很大)。令 $a_{ij}$ 表示第 $j$ 种方案里类别 $i$ 的个数， $y_j$ 表示第 $j$ 种方案的选择个数，原问题可以变为：
$min y_1+...+y_n$
$a_{1,1}y_1+a_{1,2}y_2+...+a_{1,n}y_n \geq 25$
$a_{2,1}y_1+a_{2,2}y_2+...+a_{2,n}y_n \geq 20$
$a_{3,1}y_1+a_{3,2}y_2+...+a_{3,n}y_n \geq 15$
其中
$a_{1,1}=5$ ， $a_{2,1}=0$ ， $a_{3,1}=0$
$a_{1,2}=0$ ， $a_{2,2}=2$ ， $a_{3,2}=0$
$a_{1,3}=0$ ， $a_{2,3}=0$ ， $a_{3,3}=2$
对应最初的3种方案。
$a_{i,j}$ 满足条件： $3*a_{1,j}+6*a_{2,j}+7*a_{3,j}\leq 16$ （前面所说的列生成的规则）
注意我们在使用列生成算法的时候，将整数变量松弛为了连续变量。

用矩阵形式重新表达如下：
$\min \Sigma_{j} y_j$
s.t. $Ay ≥ b^T$
$A$ 是一个矩阵， $b=[b_1,...,b_m]$ 和 $y=[y_1,...,y_n]$ 是向量；并且有重要的一点：set-covering的变量 $y_j$ 对应的列 $a_j=[a_{1j},...,a_{mj}]^T$ 生成的规则，可以通过线性约束条件 $B*a_j ≤ D$ 得到。这使得我们可以用线性规划的方法来求解子问题~
注意由于列生成算法求解的是松弛线性规划问题，因此对整数规划模型需要结合branch and bound算法进行。

2.1 限制主问题

首先用启发式算法找出一部分 $A$ ，比如说选出了 $k$ 列。然后我们的线性规划问题就变成了：
min $y_1+y_2+...+y_k$
s.t. $a_{i,1}*y_1+a_{i,2}*y_2+...+a_{i,k}*y_k≥ b_i$ ， $i = 1... n$
相比原来的模型，相当于把 $y_{k+1}$ ~ $y_m$ 强制限制为非基变量了，称为限制主问题（Restricted Master Problem，RMP）。上面的限制主问题求解完成后，我们想使用单纯型法进行基变量的转换，看看 $y_{k+1}$ ~ $y_m$ 中，是否有可以转入基变量的列。检验数 $\sigma = c_N - \pi * a_j$ ，并且 $π = c_BB^{-1}$ 可以由求解器给出。我们要找出非基变量中最小的负数 $\sigma$ ，将其转入基变量。正如前面所述，我们这里使用线性规划来找这个新的列。

2.2 子问题

注意 $c = [1, ..., 1]$ ，因此 min $c_N - \pi a_j$ 等价于
min $\pi *a$
s.t. $B * a \leq D$ (列生成的规则)
称为子问题（sub problem）。如果目标函数最优值＜0，就将新生成的列y_k+1转入基变量，生成新的限制主问题进行求解。如此往复，直至子问题的目标函数≥0。

2.3 求解质量分析

一般性的，我们假设原问题为
$Z=\max cx$
$Ax=b,x\ge0$

RMP为
$\max c\bar{x}$
$A\bar{x}=b,\bar{x}\ge0$ 。
令 $\bar{z}$ 为RMP的最优解， $z^*$ 为原问题的最优解，我们有 $z^*\le \bar{z}$

SP求得的最佳残差（负数）为 $\bar{c}^*$ ，并令 $k\ge \Sigma x$ 。 $\bar{z}$ 同样也是对偶问题的最优解，因此我们有 $\bar{z}+ \bar{c}^*k\le z^*$

在set-partitioning问题中，我们可以取 $\Sigma x=x^*$ ，因此有 $\bar{z}/(1- \bar{c}^*)\le z^*\le \bar{z}$ 。同样，我们可以参考simplex的启发算法，令
在这里插入图片描述
则

2.4 迭代过程

2.4.1 第一轮

问题对应的初始解为 $y_1=5$ ， $y_2=10$ ， $y_3=8$ ， $y_4=...y_n=0$
限制主问题（restricted master problem）如下：
$min y_1+y_2+y_3$
$y_1+0 y_2+0y_3 \geq 25$
$0y_1+2y_2+0y_3 \geq 20$
$0y_1+0y_2+2y_3 \geq 15$

下面是pymprog的sensitivity输出：
在这里插入图片描述
对应 $\pi = (0.2, 0.5, 0.5)$ ，我们把新转入的变量记作 $a_4 = (a_{1,4},a_{2,4},a_{3,4})$ ，则子问题为：
min $1 - 0.2*a_{1,4}-0.5*a_{2,4}-0.5*a_{3,4}$
s.t. $3*a_{1,4}+6*a_{2,4}+7*a_{3,4}\leq 16$ (列生成的规则)
$a_{i,j}\in Z$
求解结果为 $a_4=(1,2,0)$ ，检验数 $\sigma = 1 - \pi * a_4 =-0.2\le 0$ 。此时问题目标为22.5。
将 $y_4$ 进基，把这个结果添加到主问题当中去，开始第二轮迭代。

2.4.2 第二轮

添加了 $y_4$ 后，限制主问题为：
$min y_1+y_2+y_3+y_4$
$y_1+0 y_2+0y_3+1y_4 \geq 25$
$0y_1+2y_2+0y_3+2y_4 \geq 20$
$0y_1+0y_2+2y_3+0y_4 \geq 15$

直观上来看，是添加了一个新列，所以称为列生成算法。上面的主问题求解结果为：
在这里插入图片描述
对应 $\pi = (0.2, 0.4, 0.5)$ ，我们把新转入的变量记作 $a_5 = (a_{1,5},a_{2,5},a_{3,5})$ ，则子问题为：
min $1 - 0.2*a_{1,5}-0.4*a_{2,5}-0.5*a_{3,5}$
s.t. $3*a_{1,5}+6*a_{2,5}+7*a_{3,5}\leq 16$ (列生成的规则)
$a_{i,j}\in Z$
求解结果为 $a_5=(3,0,1)$ ，检验数 $\sigma = 1 - \pi * a_5=-0.1< 0$ 。此时问题目标为20.5。
将 $y_5$ 进基，把这个结果添加到主问题当中去，开始第三轮迭代。

2.4.3 第三轮

添加了 $y_5$ 后，限制主问题为：
$min y_1+y_2+y_3+y_4+y_5$
$y_1+0 y_2+0y_3+1y_4 +3y_5 \geq 25$
$0y_1+2y_2+0y_3+2y_4 +0y_5\geq 20$
$0y_1+0y_2+2y_3+0y_4 +1y_5\geq 15$

上面的主问题求解结果为：
在这里插入图片描述
对应 $\pi = (0.166667, 0.416667, 0.5)$ ，我们把新转入的变量记作 $a_6 = (a_{1,6},a_{2,6},a_{3,6})$ ，则子问题为：
min $1 - 0.166667*a_{1,6}-0.416667*a_{2,6}-0.5*a_{3,6}$
s.t. $3*a_{1,6}+6*a_{2,6}+7*a_{3,6}\leq 16$ (列生成的规则)
$a_{i,j}\in Z$
求解结果为 $a_6=(1,1,1)$ ，检验数 $\sigma = 1 - \pi * a_6=-0.083334< 0$ 。此时问题目标为20。
将 $y_6$ 进基，把这个结果添加到主问题当中去，开始第三轮迭代。

2.4.4 第四轮

添加了 $y_6$ 后，限制主问题为：
$min y_1+y_2+y_3+y_4+y_5+y_6$
$y_1+0 y_2+0y_3+1y_4 +3y_5+y_6 \geq 25$
$0y_1+2y_2+0y_3+2y_4 +0y_5+y_6\geq 20$
$0y_1+0y_2+2y_3+0y_4 +1y_5+y_6\geq 15$

上面的主问题求解结果为：
在这里插入图片描述
对应 $\pi = (0.2, 0.4, 0.4)$ ，我们把新转入的变量记作 $a_7 = (a_{1,7},a_{2,7},a_{3,7})$ ，则子问题为：
min $1 - 0.2*a_{1,7}-0.4*a_{2,7}-0.4*a_{3,7}$
s.t. $3*a_{1,7}+6*a_{2,7}+7*a_{3,7}\leq 16$ (列生成的规则)
$a_{i,j}\in Z$
求解结果为 $a_7=(3,1,0)$ ，检验数 $\sigma = 1 - \pi * a_6= 0$ ，结束迭代，此时问题目标为19。

2.4.5 LP和MIP的差别

我们把下面的模型称为Kantorovich模型：
在这里插入图片描述

LP问题可以给出一个下界，bound为在这里插入图片描述
这个bound很直观，即假设纸可以拼接。
我们换一个例子

使用LP得到bound为：
在这里插入图片描述
和最优值156差别还是比较大的。

2.5 Julia调用GLPK代码示例

代码如下：

A = [[5,0,0],[0,2,0],[0,0,2]]
b = [25 20 15];
while true
    # main model
    model = Model(GLPK.Optimizer)
    @variable(model, y[1:length(A[1])]>=0)
    @objective(model, Min, sum(y))
    r=@constraint(model,[i=1:3],sum(A[i].*y)>=b[i]);
    optimize!(model)
    
    # sub model
    submodel = Model(GLPK.Optimizer)
    @variable(submodel, a[1:3]>=0, Int)
    @objective(submodel, Min, 1-sum(dual.(r).*a))
    @constraint(submodel,sum([3,6,7].*a)<=16);
    optimize!(submodel)
    delta = objective_value(submodel)
	
    if delta<-0.001
        A = [[A[i];Int(JuMP.value(a[i]))] for i in 1:3]
        @info objective_value(model),"add",Int.(JuMP.value.(a))
    else
        @info objective_value(model)
        break
    end
end

经过3轮迭代得到最优解
在这里插入图片描述

2.6 python调用cplex代码示例

使用下列代码求出最优解：

from pymprog import *
p = model('example')
y = p.var('y',6,int) # variables
p.minimize(sum(y), 'master')
r1 = 5*y[0]+ y[3]+3*y[4]+y[5]>= 25 
r2 = 2*y[1]+2*y[3]+y[5] >= 20 
r3 = 2*y[2] +y[4]+y[5]>= 15 
p.solve()
for i in range(6):
    print(y[i].primal)

求解结果为： $y_1=1,y_4=3,y_5=1,y_6=14$ ，总计需要19个纸筒。
也可以将最后一轮的 $y$ 取整得到一个可行解： $y_1=2,y_4=3,y_6=15$ ，需要20个纸筒。
这里给出第四轮的限制主问题和子问题的代码供参考：

from pymprog import *
p = model('master problem')
y = p.var('y',6) # variables
p.minimize(sum(y), 'master')
r1 = 5*y[0]+ y[3]+3*y[4]+y[5]>= 25 
r2 = 2*y[1]+2*y[3]+y[5] >= 20 
r3 = 2*y[2] +y[4]+y[5]>= 15 
p.solve()
p.sensitivity()

from pymprog import *
s = model('sub problem')
a = s.var('a',3) # variables
for i in range(3):
    a[i].kind = int
s.minimize(1-(0.2*a[0] + 0.4*a[1] +0.4*a[2]), 'sub')
3*a[0] +6*a[1] +7*a[2]<= 16 # mountain bike limit
s.solve()
print(s.vobj())
for i in range(3):
    print(a[i].primal)

然后下面给出直接求最优解的代码。pymprog用了一天都没有求出来，因此改用cplex进行求解：

import cplex
from cplex.exceptions import CplexError
my_obj = [0.0,0.0,0.0,1.0]*23
my_ub = [5.0,2.0,2.0,1.0]*23
my_lb = [0.0,0.0,0.0,0.0]*23
my_ctype = "IIII"*23
my_colnames = []
my_rhs = [25,20,15]+[0,16]*23
my_sense = "GGG"+"LL"*23
for i in range(23):
    my_colnames = my_colnames + ["x0,"+str(i),"x1,"+str(i),"x2,"+str(i),"y"+str(i)]

def populatebyrow(prob):
    prob.objective.set_sense(prob.objective.sense.minimize)
    prob.variables.add(obj=my_obj, lb=my_lb, ub=my_ub, types=my_ctype,names=my_colnames)
    row0 = [1.0,0.0,0.0,0.0]*23
    row1 = [0.0,1.0,0.0,0.0]*23
    row2 = [0.0,0.0,1.0,0.0]*23    
    rows = [[my_colnames, row0],
            [my_colnames, row1],
            [my_colnames, row2]]
    for i in range(23):
        row3 = [0.0,0.0,0.0,0.0]*23
        row3[i*4:i*4+4]=[1.0,1.0,1.0,-100.0]
        row4 = [0.0,0.0,0.0,0.0]*23
        row4[i*4:i*4+4]=[3.0,6.0,7.0,0.0]
        rows = rows+[[my_colnames, row3],
            [my_colnames, row4]]
    prob.linear_constraints.add(lin_expr=rows, senses=my_sense,rhs=my_rhs)

try:
    my_prob = cplex.Cplex()
    handle = populatebyrow(my_prob)
    my_prob.solve()
    
except CplexError as exc:
    print(exc)

print("Solution status = ", my_prob.solution.status[my_prob.solution.get_status()])
print("Solution value  = ", my_prob.solution.get_objective_value())
x = my_prob.solution.get_values()
for i in range(23):
    if x[i*4+3]>0:
        print(x[i*4:i*4+3])

结果如下：

CPXPARAM_Read_DataCheck                          1
Tried aggregator 1 time.
MIP Presolve modified 23 coefficients.
Reduced MIP has 49 rows, 92 columns, and 230 nonzeros.
Reduced MIP has 23 binaries, 69 generals, 0 SOSs, and 0 indicators.
Presolve time = 0.00 sec. (0.12 ticks)
Found incumbent of value 22.000000 after 0.00 sec. (0.48 ticks)
Probing time = 0.00 sec. (0.01 ticks)
Tried aggregator 1 time.
Reduced MIP has 49 rows, 92 columns, and 230 nonzeros.
Reduced MIP has 23 binaries, 69 generals, 0 SOSs, and 0 indicators.
Presolve time = 0.00 sec. (0.21 ticks)
Probing time = 0.00 sec. (0.01 ticks)
MIP emphasis: balance optimality and feasibility.
MIP search method: dynamic search.
Parallel mode: deterministic, using up to 4 threads.
Root relaxation solution time = 0.00 sec. (0.13 ticks)

        Nodes                                         Cuts/
   Node  Left     Objective  IInf  Best Integer    Best Bound    ItCnt     Gap

*     0+    0                           22.0000        0.0000           100.00%
      0     0        6.6667    38       22.0000        6.6667       30   69.70%
      0     0       10.2791    39       22.0000      Cuts: 62      143   53.28%
      0     0       18.0972    40       22.0000      Cuts: 80      213   17.74%
*     0+    0                           21.0000       18.0972            13.82%
      0     0       18.4688    35       21.0000      Cuts: 38      305   12.05%
      0     0       18.7500    19       21.0000      Cuts: 41      395   10.71%
      0     0       18.7500    12       21.0000      Cuts: 18      444   10.71%
      0     0       18.8667    10       21.0000      Cuts: 12      496   10.16%
*     0+    0                           19.0000       18.8667             0.70%
      0     0        cutoff             19.0000       18.8667      496    0.70%
Elapsed time = 0.18 sec. (7.56 ticks, tree = 0.01 MB, solutions = 3)

Implied bound cuts applied:  22
Mixed integer rounding cuts applied:  35
Zero-half cuts applied:  5
Gomory fractional cuts applied:  7

Root node processing (before b&c):
  Real time             =    0.19 sec. (7.57 ticks)
Parallel b&c, 4 threads:
  Real time             =    0.00 sec. (0.00 ticks)
  Sync time (average)   =    0.00 sec.
  Wait time (average)   =    0.00 sec.
                          ------------
Total (root+branch&cut) =    0.19 sec. (7.57 ticks)
Solution status =  MIP_optimal
Solution value  =  19.0
[1.0, 1.0, 1.0]
[1.0, 1.0, 1.0]
[1.0, 2.0, 0.0]
[1.0, 1.0, 1.0]
[1.0, 1.0, 1.0]
[1.0, 1.0, 1.0]
[1.0, 1.0, 1.0]
[1.0, 1.0, 1.0]
[1.0, 1.0, 1.0]
[1.0, 1.0, 1.0]
[1.0, 2.0, 0.0]
[1.0, 1.0, 1.0]
[1.0, 2.0, 0.0]
[1.0, 1.0, 1.0]
[5.0, 0.0, 0.0]
[1.0, 1.0, 1.0]
[1.0, 1.0, 1.0]
[3.0, 0.0, 1.0]
[1.0, 1.0, 1.0]

最终是19个。

3. 例子2：多货物配送问题

多个货物在有路径容量限制的图上G=(V,E)进行配送。令 $x^l_{i,j}$ 为货物 $l$ 是否在路径 $i, j$ 上进行配送， $q^l$ 为货物的配送量，直观建模为：
在这里插入图片描述
$b^l_i = 1$ if $i$ is the source node of commodity $l$ ,let $b^l_i =−1$ if $i$ is the target node of commodity $l$ , and let $b^l_i = 0$ otherwise。
令P为path集合，我们使用set-partitioning模型重新建模：

在这里插入图片描述
其中 $δ^p_{ij}$ be a constant denoting whether or not path p visits edge $\in E$
模型中多了一个等式约束，并且子问题不能用线性规划简单求解。
首先我们来看列生成的规则，Let $π_{ij}$ ≤0 be the dual of constraint (19) and $σ^l ∈R$ be the dual of constraint (20)，
则reduced cost for a column p 为：
在这里插入图片描述

我们有列生成的限制要求：
在这里插入图片描述
我们可以将距离矩阵从 $c_{ij}$ 替换为 $c_{ij}-\pi_{ij}$ ，上面这个不等式的左边可以理解为寻找每个商品 $l$ 在新的距离矩阵下的最短路，然后看是否小于 $\sigma^l$ ，如果满足则将此路径p添加进限制主问题中。

4. 例子3：bin packing问题

bin的容量为 $W$ ，有 $n$ 个物品，大小为： $L={w_1,...,w_n}$ 。将物品放在bin中，使bin的数量最少。
使用flow模型进行建模：
在这里插入图片描述
如上图，W=5，L={3,2}，下半部分是一条可行路径。目标函数是路径条数最小化，约束则是每个w都要有，并且中间的点需要流量平衡。
Min z
s.t.

IE06

关注

19
点赞
踩
88

收藏

觉得还不错? 一键收藏
10
评论
运筹系列8：Set partitioning问题的列生成方法

1. 应用场景列生成算法是不断添加求解变量的算法，可参考论文。列生成算法常常用于如下的场景：使用set-covering构建的模型，变量非常多，约束相对较少。具体来说，场景如下：有nnn个0-1变量z1...znz1...znz_1...z_n，每个ziziz_i带着很多中间变量xi,jxi,jx_{i,j}用来进行约束，这是个变量很多，约束也很多的模型。我们首先使用set-coverin...
复制链接

扫一扫