線(xiàn)性規(guī)劃技巧: 列生成(Column Generation)

列生成(Column Generation)是一種把線(xiàn)性規(guī)劃問(wèn)題分解為小規(guī)模子問(wèn)題的技巧[1] [2]. 它的原理基于單純形算法. 從一個(gè)基本解(Basic Solution)出發(fā), 主問(wèn)題的系數(shù)矩陣只考慮基本解對(duì)應(yīng)的列, 然后求解子問(wèn)題來(lái)生成主問(wèn)題系數(shù)矩陣的列. 這種迭代求解的方法使得我們可以求解一些具有 指數(shù)多個(gè)變量 的線(xiàn)性規(guī)劃.

下料問(wèn)題 (The cutting stock problem)

假設(shè)我們是一家原材料生產(chǎn)工廠(chǎng)(例如紙張, 布料, 鋼管). 原材料的樣子可參考下圖.

原材料. 來(lái)自網(wǎng)絡(luò)

為了生產(chǎn)成品, 我們需要把原材料切分成更小尺寸的 成品材料(final)(示意圖如下).

切割示例

下面我們描述一個(gè)具體的例子.

原材料長(zhǎng)1000 (單位:厘米). 需要把原材料切分成4種長(zhǎng)度分別為450, 360, 310, 140的成品材料, 且滿(mǎn)足每種成品材料對(duì)應(yīng)的需求量: 97, 610, 395, 211. 我們需要計(jì)算原材料的切分方式和數(shù)量, 目標(biāo)是最小化所需原材料的數(shù)量.

數(shù)學(xué)模型

我們先定義可行的切割方式(Feasible cutting pattern).

  • 設(shè)L是原材料的長(zhǎng)度
  • 設(shè)m是成品材料種類(lèi)的總數(shù)(本例m=4)
  • 用向量s=(s_1, s_2, ..., s_m)^T代表每種成本材料的長(zhǎng)度
  • 用向量a_j=(a_{1,j}, a_{2,j}, ..., a_{m,j})代表每種成品材料的切割數(shù)量, 其中j代表切割方式

可行切割. 切割方式j可行當(dāng)且僅當(dāng)s^Ta_j \leq L.

L=1000. s=(450, 360, 310, 140). 那么a=(0, 2, 0, 2)^T代表把一根原材料切成2段360的成品材料和 2段140的成品材料.

對(duì)于每一件原材料, 算法需要給出可行的切割方式. 目標(biāo)是最小化可行切割方式的總數(shù),即消耗原材料的總數(shù). 按照這種思路我們構(gòu)建一個(gè)數(shù)學(xué)模型.

指標(biāo)

  • j - 可行的切割方式
  • i - 成品材料

參數(shù)

  • d_i - 成品材料i的需求量
  • a_{i,j} - 成品材料i在切割方式j中的數(shù)量

決策變量

  • x_j \in \mathbb{N}^+ -- 選擇切割方式j的原材料數(shù)量

目標(biāo)

  • \min \sum_j x_j

約束

  • 滿(mǎn)足成品材料的需求量: \sum_{j} a_{i,j} x_j \geq d_i

綜上所述, 我們得到如下整數(shù)線(xiàn)性規(guī)劃.
\begin{align} \min ~ & \sum_j x_j \\ \text{ s.t. } & \sum_{j} a_{i,j} x_j \geq d_i, \forall i \\ & x_j \geq 0, \text{integer}, \forall j \end{align}

求解思路

直接求解

模型的輸入需要枚舉所有的可行切割方式, 即需要構(gòu)造如下矩陣:
A = [a_1, a_2, ..., a_n],
其中n代表可行切割方式的總數(shù), a_j代表第j種可行切割方式對(duì)應(yīng)的列向量, j=1, 2, ..., n.

如果成品材料的種類(lèi)會(huì)發(fā)生改變, 這種直接計(jì)算的方式可能變得不可行. 原因是可行切割數(shù)量n隨著成品材料種類(lèi)數(shù)量m的增大呈指數(shù)增長(zhǎng), 這樣一來(lái)可能導(dǎo)致矩陣A非常大, 進(jìn)而導(dǎo)致計(jì)算過(guò)程中的內(nèi)存消耗超出限制.

間接求解

把問(wèn)題分解為主問(wèn)題(Master problem)和子問(wèn)題(Sub problem). 主問(wèn)題只考慮 一部分 可行的切割方式, 例如考慮A'=[a_1, a_2, ..., a_m]; 通過(guò)求解子問(wèn)題, 我們生成可行的切割方式a_{m+1}, 然后把它加入主問(wèn)題的輸入矩陣, 例如A' = [a_1, a_2, ..., a_{m+1}]; 再次求解主問(wèn)題; 反復(fù)迭代求解直到滿(mǎn)足停止條件. 我們把這種計(jì)算方式稱(chēng)為列生成. 這種方式的好處是把大問(wèn)題拆分成了許多小問(wèn)題, 然后通過(guò)迭代的方式來(lái)求解.

列生成

回顧單純形算法

已知線(xiàn)性規(guī)劃的標(biāo)準(zhǔn)形式如下:
\begin{align} \min ~ & c^Tx \\ \text{s.t. } & Ax = b\\ & x \geq 0 \end{align}

在單純形算法(Simplex algorithm)中, 決策變量被分成兩類(lèi): 基變量x_B和非基變量x_N. 兩類(lèi)變量對(duì)應(yīng)的矩陣我們記作BN. 因此, 線(xiàn)性規(guī)劃也可以寫(xiě)成:

\begin{align} \min ~ & c^T_B x_B + c^T_N x_N \\ \text{s.t. } & B x_B + N x_N = b \\ & x_B, x_N \geq 0. \end{align}

根據(jù)基變量的定義, 我們有x_B = B^{-1}(b - Nx_N). 因此目標(biāo)函數(shù)可以寫(xiě)成
c^T_B(B^{-1}b - Nx_N) + c^T_Nx_N = c^T_BB^{-1}b + (c^T_N - c^T_BB^{-1}N)x_N.

上述目標(biāo)函數(shù)關(guān)于b求導(dǎo), 我們得到 Shadow Price (或Dual Variable):
\lambda^T = c^T_B B^{-1}
對(duì)x_N求導(dǎo)我們得到Reduced Cost (非基變量每增加一個(gè)單位, 目標(biāo)函數(shù)的改變值):
\mu^T = c^T_N - c^T_B B^{-1}N = c^T_N - \lambda^T N.

單純形算法從一個(gè)基本解x_B^0開(kāi)始迭代, 當(dāng)存在\mu_j < 0 (\mu的第j個(gè)分量), 則把對(duì)應(yīng)的x_j加入基(basis)中, 同時(shí)剔除一個(gè)已有的基變量, 因此增加x_j可以降低目標(biāo)值. 如此迭代下去直到\mu \geq 0, 目標(biāo)值無(wú)法再降低, 因此得到最優(yōu)解.

列生成模型

本節(jié)我們考慮把原始問(wèn)題拆分成更小的主問(wèn)題和子問(wèn)題, 然后分別迭代求解. 在主問(wèn)題中, 一開(kāi)始只考慮基本解對(duì)應(yīng)的矩陣B, 然后通過(guò)求解子問(wèn)題來(lái)生成主問(wèn)題的列, 從而把它加入到主問(wèn)題中進(jìn)行求解. 根據(jù)上面的公式 \mu^T = c^T_N - \lambda^T N, 我們只需要生成 能降低目標(biāo)函數(shù)的列, 即\mu對(duì)應(yīng)的分量嚴(yán)格小于0. (因此沒(méi)有必要枚舉所有可能性.)

設(shè)s^T=(s_1, s_2, ..., s_m)代表m種成品材料的長(zhǎng)度. 令y是子問(wèn)題要生成的可行切分方式, 因此需要滿(mǎn)足條件: s^T y \leq L. 此外, 注意到主問(wèn)題的目標(biāo)函數(shù)的系數(shù)都是1, 即c_N^T = (1, 1, ... ,1), 因此y對(duì)應(yīng)的reducde cost 為\mu_y = 1 - \lambda^T y < 0. 綜上所述, 子問(wèn)題滿(mǎn)足兩個(gè)約束:

  • s^T y \leq L
  • 1-\lambda^Ty < 0

考慮到第二個(gè)約束是不等式, 我們可以把它轉(zhuǎn)換成優(yōu)化目標(biāo): \max~ \lambda^T y.

注意

  1. 當(dāng)目標(biāo)函數(shù)\lambda^T y \geq 1停止迭代, 我們得到主問(wèn)題的最優(yōu)解.
  2. 子問(wèn)題的參數(shù)\lambda來(lái)自主問(wèn)題求解過(guò)程的副產(chǎn)品(\lambda^T=c_N^TB^{-1}).

因此, 我們得到如下子問(wèn)題:
\begin{align} \max ~ & \sum_{i=1}^m {\lambda}_i y_i \\ \text{s.t. } & \sum_{i=1}^m s_i y_i \leq L\\ & y_i \geq 0, i = 1, 2, ..., m \end{align}

求解過(guò)程

注意到主問(wèn)題和子問(wèn)題都是整數(shù)規(guī)劃, 而且從計(jì)算復(fù)雜性上來(lái)看都是NP-hard問(wèn)題(證明省略, 感興趣的同學(xué)請(qǐng)參考裝箱問(wèn)題(bin packing)和背包問(wèn)題(knapsack)). 當(dāng)問(wèn)題規(guī)模很大時(shí), 求解過(guò)程依然很慢.

實(shí)際的求解過(guò)程如下所述:

  1. 把主問(wèn)題松弛成線(xiàn)性規(guī)劃.
  2. 初始化主問(wèn)題一個(gè)基本解對(duì)應(yīng)的系數(shù)矩陣. 在本例中, 令A為單位矩陣E, 即A = (e_1, e_2, ..., e_m), 其中列向量e_i的第i個(gè)分量為1, 代表把原材料切成1根成品材料i(剩下的扔掉).
  3. 迭代求解主問(wèn)題和子問(wèn)題, 直到得到主問(wèn)題的最優(yōu)解(如下所示) 注意: 此時(shí)的最優(yōu)解可能不是整數(shù)解.
DO
 solve master problem;
 solve sub problem (generate column y);
 add column y to master problem;
WHILE(OPT(sub problem) <= 1)
  1. 把主問(wèn)題的解整數(shù)化(rounding).
1. 把分?jǐn)?shù)解向下取整(round down), 并計(jì)算沒(méi)有被滿(mǎn)足的需求.
2. 把未滿(mǎn)足的成品材料按長(zhǎng)度從大到小排序.
3. 依次滿(mǎn)足上述成品材料(用直觀的方式).

主問(wèn)題 (master problem)
\begin{align} \min ~ & \sum_j x_j \\ \text{s.t. } & \sum_j a_{i,j} x_j \geq d_i, \forall i \\ & x_j \geq 0, \forall j \end{align}
子問(wèn)題(sub problem)
\begin{align} \max ~ & \sum_{i=1}^m {\lambda}_i y_i \\ \text{s.t. } & \sum_{i=1}^m s_i y_i \leq L\\ & y_i \geq 0, i = 1, 2, ..., m \end{align}

Python實(shí)現(xiàn)

主問(wèn)題模型

# master_model.py

from ortools.linear_solver import pywraplp
import numpy as np


class MasterModel(object):

    def __init__(self, A, d):
        """
        :param A: 可行切割矩陣(每列代表一種切割方式)(m*n維矩陣),
                  其中m代表成品材料類(lèi)型總數(shù), n是我們考慮的可行切割方式的數(shù)量
        :param d: 成品材料的需求量(m維向量), m代表成品材料類(lèi)型的總數(shù)
        """
        self._solver = pywraplp.Solver('MasterModel',
                                       pywraplp.Solver.CLP_LINEAR_PROGRAMMING)
        self._A = A
        self._d = d
        self._x = None  # 決策變量
        self._m, self._n = np.array(self._A).shape
        self._constraints = None  # 約束的集合(需要根據(jù)約束得到對(duì)偶變量, 別名shadow price)
        self._solution_x = None  # 計(jì)算結(jié)果
        self._sp = None  # shadow price(m維向量)

    def _init_decision_variables(self):
        self._x = [self._solver.NumVar(0, self._solver.Infinity(), "x[%d]" % j)
                   for j in range(self._n)]

    def _init_constraints(self):
        self._constraints = [None] * self._m
        for i in range(self._m):
            ct = self._solver.Constraint(self._d[i], self._d[i])
            for j in range(self._n):
                ct.SetCoefficient(self._x[j], self._A[i][j])
            self._constraints[i] = ct

    def _init_objective(self):
        obj = self._solver.Objective()
        for j in range(self._n):
            obj.SetCoefficient(self._x[j], 1)
        obj.SetMinimization()

    def solve(self):
        self._init_decision_variables()
        self._init_constraints()
        self._init_objective()
        self._solver.Solve()
        self._solution_x = [self._x[j].solution_value() for j in range(self._n)]
        # shadow price
        self._sp = [self._constraints[i].dual_value() for i in range(self._m)]

    def print_info(self):
        print("[Master problem info]")
        print(" - Objective value =", self.get_objective_value())
        print(" - Shadow price: lambda = ", self._sp)

    def get_sp(self):
        """ 得到shadow price
        """
        return self._sp

    def get_solution(self):
        return self._solution_x

    def get_objective_value(self):
        return sum(self._solution_x)

子問(wèn)題模型

# sub_model.py

from ortools.linear_solver import pywraplp


class SubModel(object):

    def __init__(self, sp, s, L):
        """
        :param sp: 主問(wèn)題的shadow price(m維度向量), m代表成品材料的類(lèi)型總數(shù)
        :param s: 成品材料的長(zhǎng)度(m維向量)
        :param L: 原材料的總長(zhǎng)度
        """
        self._solver = pywraplp.Solver('MasterModel',
                                       pywraplp.Solver.CBC_MIXED_INTEGER_PROGRAMMING)
        self._sp = sp
        self._s = s
        self._L = L
        self._m = len(s)
        self._y = None  # 決策變量
        self._solution_y = None  # 計(jì)算結(jié)果

    def _init_decision_variables(self):
        self._y = [self._solver.IntVar(0, self._solver.Infinity(), "y[%d]" % i)
                   for i in range(self._m)]

    def _init_constraints(self):
        ct = self._solver.Constraint(0, self._L)
        for i in range(self._m):
            ct.SetCoefficient(self._y[i], self._s[i])

    def _init_objective(self):
        obj = self._solver.Objective()
        for i in range(self._m):
            obj.SetCoefficient(self._y[i], self._sp[i])
        obj.SetMaximization()

    def solve(self):
        self._init_decision_variables()
        self._init_constraints()
        self._init_objective()
        self._solver.Solve()
        self._solution_y = [self._y[i].solution_value() for i in range(self._m)]

    def print_info(self):
        print("[Sub problem info]")
        print(" - Reduced cost of master problem =", 1 - self.get_objective_value())
        print(" - New column generated:", self._solution_y)

    def get_solution(self):
        return self._solution_y

    def get_objective_value(self):
        return sum(map(lambda x: x[0]*x[1], zip(self._sp, self._solution_y)))

列生成迭代求解的流程

# cg_proc.py

import numpy as np

from master_model import MasterModel
from sub_model import SubModel


class CGProc(object):
    """
    Column Generation Process
    """
    def __init__(self, s, d, L, max_iter=10000):
        """
        
        :param s: 每種成品材料的長(zhǎng)度(m維向量)
        :param d: 每種成品材料的需求量(m維向量)
        :param L: 原材料的長(zhǎng)度
        :param max_iter: 最大循環(huán)次數(shù)
        """
        self._s = s
        for val in s:
            if val > L:
                raise ValueError("final size cannot be larger than raw size!")
        self._d = d  
        self._L = L  
        self._reduced_cost = -1
        self._iter_times = 0
        self._max_iter = max_iter
        self._status = -1  # -1:執(zhí)行錯(cuò)誤; 0:最優(yōu)解; 1: 達(dá)到最大循環(huán)次數(shù)
        self._A = None  # Master problem的輸入(可行切割矩陣)
        self._obj_val = None  # 目標(biāo)函數(shù)值
        self._x = None  # Master problem對(duì)應(yīng)的solution value
        self._solution_x = None  # 切割方式j(luò)需要的原材料數(shù)量
        self._solution_matrix = None  # 切割方式j(luò)

    def _stop_criteria_is_satisfied(self):
        """ 根據(jù)reduced cost判斷是否應(yīng)該停止迭代.
        """
        if self._reduced_cost > -1e-6:
            self._status = 0
            return True
        if self._iter_times >= self._max_iter:
            if self._status == -1:
                self._status = 1
            return True
        return False

    def _init_basic_matrix(self):
        # 生成單位矩陣
        self._A = np.identity(len(self._s))

    def add_column(self, y):
        a = np.array(self._A)
        c = np.array(y).transpose()
        self._A = np.c_[a, c]

    def solve(self):
        print("==== iter %d ====" % self._iter_times)
        self._init_basic_matrix()
        mp = MasterModel(self._A, self._d)
        mp.solve()
        mp.print_info()
        self._iter_times += 1
        while not self._stop_criteria_is_satisfied():
            # 1. 解Sub problem
            print("==== iter %d ====" % self._iter_times)
            sm = SubModel(mp.get_sp(), self._s, self._L)
            sm.solve()

            # 2. 把生成的列加入到Master problem
            self.add_column(sm.get_solution())

            # 3. 解Master problem
            mp = MasterModel(self._A, self._d)
            mp.solve()
            self._x = mp.get_solution()
            self._obj_val = mp.get_objective_value()
            mp.print_info()

            # 4. 更新 reduced cost
            self._reduced_cost = 1 - sm.get_objective_value()
            sm.print_info()
            self._iter_times += 1

        self._get_solution()
        status_str = {-1: "error", 0: "optimal", 1: "attain max iteration"}
        print(">>> Terminated. Status:", status_str[self._status])

    def _get_solution_indices(self):
        # 計(jì)算有效的indices
        return [i for i in range(len(self._x)) if self._x[i] > 0]

    def _get_solution(self):
        """ 剔除冗余的結(jié)果. 這一步不是必須, 僅用來(lái)顯示結(jié)果.
        """
        indices = self._get_solution_indices()
        self._solution_x = [self._x[i] for i in indices]
        self._solution_matrix = np.array([self._A[:, i] for i in indices]).transpose()

    def print_info(self):
        print("==== Column Generation Model Summary ====")
        print(" - Objective: Total number of raws needed =", self._obj_val)
        print(" - Fractional solution: The number of raws needed for each cutting pattern\n", self._solution_x)
        print(" - Solution matrix: Cutting patterns\n", self._solution_matrix)

    def get_solution_x(self):
        return self._solution_x

    def get_solution_matrix(self):
        return self._solution_matrix

把分?jǐn)?shù)解(Fractional solution)取整的算法如下.

# rounding_proc.py

import numpy as np
from copy import deepcopy


class RoundingProc(object):

    def __init__(self, A, x, s, d, L):
        """
        :param A: 可行切割矩陣(每列代表一種切割方式)(m*n維矩陣),
                  其中m代表成品材料類(lèi)型總數(shù), n是我們考慮的可行切割方式的數(shù)量
        :param x: 切割方式需要的原材料數(shù)量(fractional)(n維向量)
        :param s: 成品材料的尺寸(m維向量)
        :param d: 成品材料的尺寸對(duì)應(yīng)的需求量(m維向量)
        :param L: 原材料的長(zhǎng)度
        """
        self._A = np.array(A)
        self._x = x
        self._s = s
        self._d = d
        self._L = L
        self._d0 = None  # 被滿(mǎn)足的需求量
        self._d1 = None  # 未被滿(mǎn)足的需求量
        self._greedy_x = None  # 貪心算法的結(jié)果: 對(duì)應(yīng)滿(mǎn)足d1的部分
        self._greedy_matrix = None  # 貪心算法的結(jié)果: 對(duì)應(yīng)滿(mǎn)足d1的部分
        self._solution_x = None  # 最終結(jié)果: 切割方式對(duì)應(yīng)的數(shù)量
        self._solution_matrix = None  # 最終結(jié)果: 切割方式對(duì)應(yīng)的矩陣

    def _round_down(self):
        """ 把分?jǐn)?shù)解self._x向下取整, 然后計(jì)算未被滿(mǎn)足的需求量.
        """
        self._x = list(map(int, self._x))

        # 計(jì)算被滿(mǎn)足的需求量
        def cal_d0(i):
            return sum(self._A[i] * np.array(self._x))

        self._d0 = [cal_d0(i) for i in range(len(self._d))]

        # 計(jì)算未被滿(mǎn)足的需求量
        self._d1 = (np.array(self._d) - np.array(self._d0)).tolist()

    def _satisfy(self):
        """ 用貪心的方式滿(mǎn)足需求
        :return: 切割方式矩陣
        """
        # 把index按成品材料的長(zhǎng)度從大到小排序
        sorted_indices = list(np.argsort(-np.array(self._s)))
        rows = []
        x = []
        d1 = deepcopy(self._d1)
        while sum(d1) > 0:
            c = self._greedy_cut(d1, sorted_indices)
            if not rows:
                rows.append(c)
                x.append(1)
            else:
                if rows[-1] != c:
                    rows.append(c)
                    x.append(1)
                else:
                    x[-1] += 1
            # 更新需求
            d1 = (np.array(d1) - np.array(c)).tolist()

        return x, np.array(rows).transpose()

    def _greedy_cut(self, d1, sorted_indices):
        """ 用貪心的方式切割1根原材料.
        :param d1: 未被滿(mǎn)足的需求
        :param sorted_indices: 成品材料的長(zhǎng)度從大到小排序?qū)?yīng)的indices
        :return: 切割方式,m維向量(m=成品材料的個(gè)數(shù))
        """
        c = [0] * len(self._d)
        raw_len = self._L
        for i in sorted_indices:
            c[i] = min(raw_len // self._s[i], d1[i])
            raw_len -= c[i] * self._s[i]
        return c

    def solve(self):
        self._round_down()
        self._greedy_x, self._greedy_matrix = self._satisfy()
        self._solution_x = self._x + self._greedy_x
        self._solution_matrix = np.c_[self._A, self._greedy_matrix]

    def print_info(self):
        print("==== Rounding Process Summary ====")
        print("[Round down result]")
        print(" - Integer solution\n", self._x)
        print(" - Unsatisfied demands\n", self._d1)
        print("[Greedy result]")
        print(" - solution (w.r.t. unsatisfied demand)\n", self._greedy_x)
        print(" - Solution matrix(w.r.t. unsatisfied demand)\n", self._greedy_matrix)
        print("[Final result]")
        print(" - Objective value: The total number of raws needed =", sum(self._solution_x))
        print(" - Solution: The number of raws needed for each cutting pattern \n", self._solution_x)
        print(" - Solution matrix: Cutting patterns\n", self._solution_matrix)
        d0 = [sum(self._solution_matrix[i] * np.array(self._solution_x)) for i in range(len(self._d))]
        print(" - Satisfied demands: \n", d0)

主程序

# main.py

from cg_proc import CGProc
from data import s, d, L   # cutting-stock instance
from rouding_proc import RoundingProc


if __name__ == '__main__':
    # 1. 用列生成模型求解原問(wèn)題的LP松弛問(wèn)題
    c = CGProc(s, d, L)
    c.run()
    c.print_info()
    # 2. 得到松弛問(wèn)題的解
    A = c.get_solution_matrix()
    x = c.get_solution_x()
    # 3. 對(duì)分?jǐn)?shù)解向下取整,然后用直觀的方式滿(mǎn)足需求,得到最終的解
    r = RoundingProc(A, x, s, d, L)
    r.solve()
    r.print_info()

完整代碼

參考文獻(xiàn)


  1. P.C. Gilmore and R.E. Gomory. A linear programming approach to the
    cutting stock problem part i.
    Operations Research 9, 849–859, 1961. ?

  2. P.C. Gilmore and R.E. Gomory. A linear programming approach to the
    cutting stock problempart ii.
    Operations Research 11, 863–888, 1963. ?

最后編輯于
?著作權(quán)歸作者所有,轉(zhuǎn)載或內(nèi)容合作請(qǐng)聯(lián)系作者
【社區(qū)內(nèi)容提示】社區(qū)部分內(nèi)容疑似由AI輔助生成,瀏覽時(shí)請(qǐng)結(jié)合常識(shí)與多方信息審慎甄別。
平臺(tái)聲明:文章內(nèi)容(如有圖片或視頻亦包括在內(nèi))由作者上傳并發(fā)布,文章內(nèi)容僅代表作者本人觀點(diǎn),簡(jiǎn)書(shū)系信息發(fā)布平臺(tái),僅提供信息存儲(chǔ)服務(wù)。

相關(guān)閱讀更多精彩內(nèi)容

友情鏈接更多精彩內(nèi)容