列生成(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)(例如紙張, 布料, 鋼管). 原材料的樣子可參考下圖.

為了生產(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è)
是原材料的長(zhǎng)度
- 設(shè)
是成品材料種類(lèi)的總數(shù)(本例
)
- 用向量
代表每種成本材料的長(zhǎng)度
- 用向量
代表每種成品材料的切割數(shù)量, 其中
代表切割方式
可行切割. 切割方式可行當(dāng)且僅當(dāng)
.
例 =1000.
. 那么
代表把一根原材料切成2段360的成品材料和 2段140的成品材料.
對(duì)于每一件原材料, 算法需要給出可行的切割方式. 目標(biāo)是最小化可行切割方式的總數(shù),即消耗原材料的總數(shù). 按照這種思路我們構(gòu)建一個(gè)數(shù)學(xué)模型.
指標(biāo)
-
- 可行的切割方式
-
- 成品材料
參數(shù)
-
- 成品材料
的需求量
-
- 成品材料
在切割方式
中的數(shù)量
決策變量
-
-- 選擇切割方式
的原材料數(shù)量
目標(biāo)
約束
- 滿(mǎn)足成品材料的需求量:
綜上所述, 我們得到如下整數(shù)線(xiàn)性規(guī)劃.
求解思路
直接求解
模型的輸入需要枚舉所有的可行切割方式, 即需要構(gòu)造如下矩陣:
其中代表可行切割方式的總數(shù),
代表第
種可行切割方式對(duì)應(yīng)的列向量,
.
如果成品材料的種類(lèi)會(huì)發(fā)生改變, 這種直接計(jì)算的方式可能變得不可行. 原因是可行切割數(shù)量隨著成品材料種類(lèi)數(shù)量
的增大呈指數(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)題只考慮 一部分 可行的切割方式, 例如考慮; 通過(guò)求解子問(wèn)題, 我們生成可行的切割方式
, 然后把它加入主問(wèn)題的輸入矩陣, 例如
; 再次求解主問(wèn)題; 反復(fù)迭代求解直到滿(mǎn)足停止條件. 我們把這種計(jì)算方式稱(chēng)為列生成. 這種方式的好處是把大問(wèn)題拆分成了許多小問(wèn)題, 然后通過(guò)迭代的方式來(lái)求解.
列生成
回顧單純形算法
已知線(xiàn)性規(guī)劃的標(biāo)準(zhǔn)形式如下:
在單純形算法(Simplex algorithm)中, 決策變量被分成兩類(lèi): 基變量和非基變量
. 兩類(lèi)變量對(duì)應(yīng)的矩陣我們記作
和
. 因此, 線(xiàn)性規(guī)劃也可以寫(xiě)成:
根據(jù)基變量的定義, 我們有. 因此目標(biāo)函數(shù)可以寫(xiě)成
上述目標(biāo)函數(shù)關(guān)于求導(dǎo), 我們得到 Shadow Price (或Dual Variable):
對(duì)求導(dǎo)我們得到Reduced Cost (非基變量每增加一個(gè)單位, 目標(biāo)函數(shù)的改變值):
單純形算法從一個(gè)基本解開(kāi)始迭代, 當(dāng)存在
(
的第
個(gè)分量), 則把對(duì)應(yīng)的
加入基(basis)中, 同時(shí)剔除一個(gè)已有的基變量, 因此增加
可以降低目標(biāo)值. 如此迭代下去直到
, 目標(biāo)值無(wú)法再降低, 因此得到最優(yōu)解.
列生成模型
本節(jié)我們考慮把原始問(wèn)題拆分成更小的主問(wèn)題和子問(wèn)題, 然后分別迭代求解. 在主問(wèn)題中, 一開(kāi)始只考慮基本解對(duì)應(yīng)的矩陣, 然后通過(guò)求解子問(wèn)題來(lái)生成主問(wèn)題的列, 從而把它加入到主問(wèn)題中進(jìn)行求解. 根據(jù)上面的公式
, 我們只需要生成 能降低目標(biāo)函數(shù)的列, 即
對(duì)應(yīng)的分量嚴(yán)格小于0. (因此沒(méi)有必要枚舉所有可能性.)
設(shè)代表
種成品材料的長(zhǎng)度. 令
是子問(wèn)題要生成的可行切分方式, 因此需要滿(mǎn)足條件:
. 此外, 注意到主問(wèn)題的目標(biāo)函數(shù)的系數(shù)都是1, 即
, 因此
對(duì)應(yīng)的reducde cost 為
. 綜上所述, 子問(wèn)題滿(mǎn)足兩個(gè)約束:
考慮到第二個(gè)約束是不等式, 我們可以把它轉(zhuǎn)換成優(yōu)化目標(biāo): .
注意
- 當(dāng)目標(biāo)函數(shù)
停止迭代, 我們得到主問(wèn)題的最優(yōu)解.
- 子問(wèn)題的參數(shù)
來(lái)自主問(wèn)題求解過(guò)程的副產(chǎn)品(
).
因此, 我們得到如下子問(wèn)題:
求解過(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ò)程如下所述:
- 把主問(wèn)題松弛成線(xiàn)性規(guī)劃.
- 初始化主問(wèn)題一個(gè)基本解對(duì)應(yīng)的系數(shù)矩陣. 在本例中, 令
為單位矩陣
, 即
, 其中列向量
的第
個(gè)分量為1, 代表把原材料切成1根成品材料
(剩下的扔掉).
- 迭代求解主問(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)
- 把主問(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)
子問(wèn)題(sub problem)
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()