列生成算法原理及实现

发表于 更新于 分类于 阅读次数：

列生成算法(Column Generation)是一种用于高效求解大规模线性优化问题的算法，其理论基础是由Danzig等于1960年提出。本质上而言，列生成算法是单纯形算法的一种形式。列生成算法被应用于求解以下著名的NP-hard优化问题：人员调度问题、切割问题、车辆路径问题等。

线性规划模型

单纯形法

对于标准线性规划问题而言， \[ \max z = \sum^{n}_{j=1} c_{j} x_{j} \] subject to \[ \sum^{n}_{j=1} a_{ij}x_{j} = b_{i} \qquad i=1,2, \cdots, m \]

\[ x_{j} \geq 0 \qquad j=1,2, \cdots, n \]

采用单纯形法求解过程时，变量\(x_{j}\)的检验数为 \[ \sigma_{j} = c_{j} - \sum_{i=1}^{m} c_{i} a_{ij} \] \(\sigma_{j} > 0\)表示非基变量\(x_{j}\)变为基变量会引起目标函数的增加，因此通常选 择\(\sigma_{j}\)最大的非基变量作为进基变量。对于最小化问题，\(\sigma_{j} < 0\)表示非基变量\(x_{j}\)变为基变量会引起目标函数的减小。

单纯形的矩阵表示为：

单纯形计算时，总是选取\(I\)为初始基，对应的基变量为\(X_s\)。迭代若干步后，基变量为 \(X_B\)，\(X_B\)在初始单纯形表中的系数矩阵为\(B\)。将\(B\)在初始单纯形中单独列出，而 \(A\)中去掉\(B\)的若干列后剩下的列组成矩阵\(N\)，这样的初始单纯形表可以列成如下形式 ：

非基变量 基变量
\(C_B\) 基 \(b\) \(X_{B}\) \(X_N\) \(X_S\)
0 \(X_S\) \(b\) \(B\) \(N\) \(I\)
\(c_{j} - z_{j}\) \(C_{B}\) \(C_{N}\) 0

迭代若干步后，基变量变为\(X_{B}\)时，该步的单纯形表中由\(X_{B}\)系数组成的矩阵为 \(I\)。此时，新的单纯形表具有如下形式：

基变量 非基变量
\(C_B\) 基 \(b\) \(X_{B}\) \(X_N\) \(X_S\)
\(C_B\) \(X_s\) \(B^{-1}b\) \(I\) \(B^{-1}N\) \(B^{-1}\)
\(c_{j} - z_{j}\) 0 \(C_{N}-C_{B}B^{-1}N\) \(-C_{B}B^{-1}\)

变量的 检验数的向量表示为 \[ \sigma_{j} = c_{j} - C_{B} B^{-1}P_{j} \] 其中,\(B\)为基变量；\(C_{B}\)为基变量在目标函数中系数组成的向量；\(P_{j}\)为系数矩阵中的第\(j\)列。

对偶问题

原问题的对偶问题为 \[ \min w = \sum_{i=1}^{m} b_{i} y_{i} \] subject to \[ \sum_{i=1}^{n} a_{ij}y_{i} \geq c_{j} \qquad j=1,2,\cdots, n \]

\[ y_{i} \geq 0 \qquad i=1,2, \cdots, m \]

根据原问题和对偶问题的关系，检验数\(\sigma_{j}\)可以写为 \[ \sigma_{j} = c_{j} - C_{B} B^{-1} P_{j} = c_{j} - \sum_{i=1}^{m} a_{ij}y_{i} \]

检验数的经济含义为reduced cost，他表示将该变量由非基变量变为基变量带来的目标函数的改进。对于最大化问题，如果所有变量的reduced cost<= 0，则表示已经无法通过换入换出基变量使目标函数变大；相反对于最小化问题，如果所有变量的reduced cost >=0，则表示已经无法通过换入换出基变量使目标函数值更小。

Cutting Stock Problem

有一堆固定长度为16m的原纸卷，顾客们分别需要25个3m长，20个6m长和18个7m长的纸卷。问怎样切割才能使得浪费最小？

分析：顾客总共需要3种类别，用\(j\)对类别进行索引。原纸卷长度16m，可以切割为5个3m 长的，也可以切割为2个6m长，甚至更多种且法方法。我们将这些不同的切割方法称为 切割方案，即一张16m的原纸卷要切割成几张3m、几张6m和几张7m的，用\(i\)对方案进行 索引。在切割方案\(i\)中类别\(j\)的数量记为\(a_{ij}\)。

定义决策变量\(x_{i}\)为采用第\(i\)种切割方案的原纸卷数，可以建立如下数学规划模型： \[ \min \sum_{i} x_{j} \] subject to \[ \sum_{i} a_{ij} x_{i} \geq d_{j} \qquad \forall j \]

\[ x_{i} \geq 0 \]

我们将此问题称为主问题（Master Problem）。 很难直接对上述模型进行求解，因为无论是目标函数还是约束条件都需要知道所有的切割 方案。按照单纯形法，在列出所有切割方案的情况下，我们才能计算每种切割方案（变量 \(x_{j}\)）的reduced cost。

如果我们能够找到一个初始基，而要对这个解进行改进就是要找出一个进基变量，其 reduced cost应该是最大的。

很明显， \[ \begin{bmatrix} 5 & 0 & 0\\ 0 & 2 & 0\\ 0 & 0 & 2\\ \end{bmatrix} \] 可以是一个初始基，基于初始可行基可将模型写为： \[ \min f = x_1 + x_2 + x_3 \] subject to \[ 5 x_{1} \geq 25 \]

\[ 2 x_{2} \geq 20 \]

\[ 2 x_{3} \geq 18 \]

将基于初始可行基构建的模型称为RMP模型（Restricted Master Problem） 该模型的对偶问题为： \[ \max w = 25 y_1 + 20 y_2 + 18 y3 \] subject to \[ 5 y_1 \leq 1 \]

\[ 2 y_2 \leq 1 \]

\[ 2 y_3 \leq 1 \]

将该模型称为DMP（Dural Master Problem）。 求解此对偶问题，可以得到\(y^{*}\)，此时寻找检验数最小的换入变量可以通过求解规划 问题： \[ \min c_{j} - \sum_{i} a_{ij}y_{i} = 1 - \sum_{i} a_{ij} y^{*}_{i} \] subject to \[ 3 a_{1j} + 6 a_{2j} + 7 a_{3j} \leq 16 \]

此模型称为Pricing Prolbem，求解此模型可以得到一种新的切割模式（\(a_{ij}\)），因此增加新的决策变量\(x_{j}\)，这也是列生成算法名称的由来。将新生成的变量加入到RMP中，并不断循环直至无法得到新的变量，即求解的最优解。

列生成算法的主要框架如下：

CPLEX实现

借助于CPLEX的docplex可以快速实现列生成算法求解切割问题，代码如下：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
#! /usr/bin/env python
# -*- coding: utf-8 -*-
# vim:fenc=utf-8
#
# Copyright © 2020 Bruce Han <[email protected]>
#
# Distributed under terms of the MIT license.
"""
Column generation with Docplex
"""
from docplex.mp.model import Model


class CutModel():
    """Cut Problem model solved with column generation"""

    def __init__(self, size, amount, width):
        """Define the cut problem."""
        self.size = size
        self.amount = amount
        self.width = width

        self.pricing_prob = None
        self.pattens = [[0 for _ in range(len(size))]
                        for _ in range(len(size))]
        self.rmp, self.vars, self.ctns = self.create_rmp(size, amount, width)

    def create_rmp(self, size, amount, width):
        """Create the restricted master problem"""
        rmp = Model("rmp")
        var = rmp.continuous_var_list(range(len(size)), name="x")

        rmp.minimize(rmp.sum(var))

        for i in range(len(self.size)):
            self.pattens[i][i] = width // size[i]
        ctns = rmp.add_constraints([(width // size[i]) * var[i] >= amount[i]
                                    for i in range(len(size))])

        return rmp, var, ctns

    def create_pricing_problem(self, size, prices, width):
        """Create the prices problem mdl with prices, and return mdl"""
        mdl = Model("pricing")
        mdl.a = mdl.integer_var_list(
            range(len(size)),
            ub=[width // size[i] for i in range(len(size))],
            name="a")
        mdl.minimize(1 - mdl.dot(mdl.a, prices))
        mdl.add_constraint(mdl.dot(mdl.a, size) <= width)

        return mdl

    def update_pricing_problem(self, prices):
        """Update the prices problem mdl with prices"""
        expr = self.pricing_prob.objective_expr
        updated = [(var, -prices[u])
                   for u, var in enumerate(self.pricing_prob.a)]
        expr.set_coefficients(updated)

    def print_sol(self, sol):
        """Print the cut details in solution sol"""
        print("Current need %d" % sol.objective_value)
        for i, p in enumerate(self.pattens):
            if sol[self.vars[i]] > 0:
                print("\t", p, sol[self.vars[i]])

    def solve(self):
        """Solve the prolbem"""
        iter = 0
        while True:
            iter += 1
            self.rmp.export_as_lp("rmp_%d.lp" % iter)

            sol_rmp = self.rmp.solve()
            if not sol_rmp:
                print("Error, can't find solution.")
                break
            else:
                print("The current cut method is:")
                self.print_sol(sol_rmp)

            prices = self.rmp.dual_values(self.ctns)
            if self.pricing_prob:
                self.update_pricing_problem(prices)
            else:
                self.pricing_prob = self.create_pricing_problem(
                    self.size, prices, self.width)

            self.pricing_prob.export_as_lp("pricing_%d.lp" % iter)
            sol = self.pricing_prob.solve()

            new_pattern = None
            if sol:
                if sol.objective_value < 0:
                    print("The retuced cost is %6.3f" % sol.objective_value)
                    new_pattern = [
                        sol[self.pricing_prob.a[i]]
                        for i in range(len(self.size))
                    ]

            if new_pattern:
                print("New patten found", new_pattern)
                self.pattens.append(new_pattern)
                # Add one new variables, i.e., a new patten
                new_var = self.rmp.continuous_var(name="x_%d" % len(self.vars))
                self.vars.append(new_var)
                # Update constraints
                for n, ct in enumerate(self.ctns):
                    ct.lhs += new_pattern[n] * self.vars[-1]
                # Update objective function
                obj_expr = self.rmp.get_objective_expr()
                obj_expr += self.vars[-1]
            else:
                # No new patten, the best scheme found
                print("The best cut method is:")
                self.print_sol(sol_rmp)
                break


def main():
    width = 218
    size = (81, 70, 68)
    amount = (44, 3, 48)

    m = CutModel(size, amount, width)
    m.solve()

    w2 = 16
    s2 = (3, 6, 7)
    a2 = (25, 20, 18)
    m2 = CutModel(s2, a2, w2)
    m2.solve()


if __name__ == "__main__":
    main()